all_together.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "33881270-367c-45e0-a14c-8c66933eb5a2",
   "metadata": {},
   "source": [
    "# era5land download for test years\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "620a6832-49a7-4def-92a4-774d79b34292",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import cdsapi\n",
    "\n",
    "# Replace 'YOUR_CDS_API_KEY' with your actual API key\n",
    "api_key = '123057:e2f4481e-d1a8-412c-a1b4-fab7dd9b3931'\n",
    "\n",
    "# Create a CDS API client\n",
    "c = cdsapi.Client(key=api_key, url=\"https://cds.climate.copernicus.eu/api/v2\")\n",
    "\n",
    "\n",
    "# Define the variables and other parameters\n",
    "variables = [\"evaporation_from_bare_soil\", \"evaporation_from_vegetation_transpiration\", \"soil_temperature_level_1\", \"soil_temperature_level_2\"]\n",
    "\n",
    "years = [2005, 2007, 2009, 2011, 2013, 2015, 2017, 2019, 2021, 2023]\n",
    "area = [90, -180, -90, 180]  # North, West, South, East bounding box coordinates\n",
    "hours = [\"00:00\", \"01:00\", \"02:00\", \"03:00\", \"04:00\", \"05:00\", \"06:00\", \"07:00\", \"08:00\", \"09:00\", \"10:00\", \"11:00\",\n",
    "         \"12:00\", \"13:00\", \"14:00\", \"15:00\", \"16:00\", \"17:00\", \"18:00\", \"19:00\", \"20:00\", \"21:00\", \"22:00\", \"23:00\"]\n",
    "\n",
    "# Loop over each variable\n",
    "for variable in variables:\n",
    "    # Loop over each year\n",
    "    for year in years:\n",
    "        # Define the request parameters with loops inside\n",
    "        request_params = {\n",
    "            \"product_type\": \"reanalysis\",\n",
    "            \"format\": \"netcdf\",\n",
    "            \"variable\": [variable],\n",
    "            \"year\": [str(year)],\n",
    "            \"month\": [f\"{month:02d}\" for month in range(1, 13)],\n",
    "            \"day\": [f\"{day:02d}\" for day in range(1, 32)],\n",
    "            \"area\": area,\n",
    "            \"time\": hours,\n",
    "        }\n",
    "\n",
    "        # Make the request and download the data\n",
    "        output_file = f\"era5_land_{variable}_{year}.nc\"\n",
    "        c.retrieve(\"reanalysis-era5-land\", request_params, output_file)\n",
    "\n",
    "# Close the CDS API client\n",
    "c.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd7e6f35-dcc6-4080-a689-af932300dd25",
   "metadata": {},
   "source": [
    "# making time series per station for all era5 variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8f0c7a3-7b10-4dbc-80ee-ff00de6fd15b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Что нужно сделать:\n",
    "# По ERA5-land для 26 постов нужно посчитать среднее значение параметра по водосборам. \n",
    "# Всего параметров 23 из списка ниже, отмечено зеленым. \n",
    "# Период расчета – с 1990 года по настоящее время. \n",
    "# dНа выходе нужны осредненные по водосборам суточные ряды по каждому параметру\n",
    "# т.е. таблица вида «site_id - date – mean_value». \n",
    "\n",
    "# read shp files \n",
    "path = '/Users/varyabazilova/Desktop/drivendata_rodeo/data/'\n",
    "gdf = gpd.read_file(path + 'geospatial.gpkg', layer='basins')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1996410e-8679-4767-b2cb-d147c0c46253",
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports \n",
    "\n",
    "import xarray as xr\n",
    "import rioxarray as rio\n",
    "import pandas as pd\n",
    "import numpy as np \n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "# correct longitudes\n",
    "# test['longitude'] = np.arange(-122.4, -104.7 + 0.1, 0.1)\n",
    "\n",
    "\n",
    "# ---------------------------------------------------\n",
    "# ------- functions -----------------\n",
    "# ---------------------------------------------------\n",
    "\n",
    "# stat by shape for the variables that are !stored normally!\n",
    "\n",
    "def stat_by_shape(xr_file, shape, kind='mean'):\n",
    "    '''calculate zonal statistics for each polygon\n",
    "    - write the \"right\" coordinates to the .nc file\n",
    "    - write the same-as-polygons crs\n",
    "    - clip the .nc file with the polygon bounds\n",
    "    - resample hourly data to daily, remove useles variables\n",
    "    - compute statistics\n",
    "    - convert to pandas df, return df\n",
    "    '''\n",
    "    xr_file['longitude'] = np.arange(-122.4, -104.7 + 0.1, 0.1)\n",
    "    xr_file = xr_file.rio.write_crs(shape.crs)\n",
    "    xr_file = xr_file.rio.clip(shape.geometry)\n",
    "    daily_mean = xr_file.resample(time='D').mean().drop_vars('spatial_ref')\n",
    "    \n",
    "    if kind=='mean':\n",
    "        sr_pcpn = daily_mean.mean(dim=['latitude', 'longitude']).to_dataframe()\n",
    "    if kind=='min':\n",
    "        sr_pcpn = daily_mean.min(dim=['latitude', 'longitude']).to_dataframe()\n",
    "    if kind=='max':\n",
    "        sr_pcpn = daily_mean.max(dim=['latitude', 'longitude']).to_dataframe()\n",
    "    if kind=='std':\n",
    "        sr_pcpn = daily_mean.std(dim=['latitude', 'longitude']).to_dataframe()\n",
    "    return sr_pcpn\n",
    "\n",
    "\n",
    "\n",
    "# -------- stat by shape for the variables that are !stored cumulatively! for each forecasting period \n",
    "# (e.g. taking MIN per day for the evaporation)\n",
    "\n",
    "def stat_by_shape_cumul(xr_file, shape, kind='mean'):\n",
    "    '''same as 'stat by shape', \n",
    "    but to resample to dayly - take MIN \n",
    "    (total value per day)\n",
    "    '''\n",
    "    xr_file['longitude'] = np.arange(-122.4, -104.7 + 0.1, 0.1)\n",
    "    xr_file = xr_file.rio.write_crs(shape.crs)\n",
    "    xr_file = xr_file.rio.clip(shape.geometry)\n",
    "    daily_mean = xr_file.resample(time='D').min().drop_vars('spatial_ref') # !!!\n",
    "    \n",
    "    if kind=='mean':\n",
    "        sr_pcpn = daily_mean.mean(dim=['latitude', 'longitude']).to_dataframe()\n",
    "    if kind=='min':\n",
    "        sr_pcpn = daily_mean.min(dim=['latitude', 'longitude']).to_dataframe()\n",
    "    if kind=='max':\n",
    "        sr_pcpn = daily_mean.max(dim=['latitude', 'longitude']).to_dataframe()\n",
    "    if kind=='std':\n",
    "        sr_pcpn = daily_mean.std(dim=['latitude', 'longitude']).to_dataframe()\n",
    "    return sr_pcpn\n",
    "\n",
    "\n",
    "\n",
    "# ----- open all .nc files in a folder, clean the 'expver' dimention\n",
    "\n",
    "def open_and_combine_nc_files(folder_path):\n",
    "    \"\"\"\n",
    "    open all .nc files togetehr\n",
    "    required: folder_path (str)\n",
    "    return: combined dataset with combined 'expver' dim\n",
    "    \"\"\"\n",
    "    # open all data together\n",
    "    dataset = xr.open_mfdataset(os.path.join(folder_path, '*.nc'))\n",
    "\n",
    "    # list of keys (dimensions)\n",
    "    keys_list = list(dataset.dims.keys())\n",
    "    # check if 'expver' is there, if so - combine\n",
    "    if 'expver' in keys_list:\n",
    "        dataset = dataset.sel(expver=1).combine_first(dataset.sel(expver=5))\n",
    "\n",
    "    return dataset\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13c253f2-eb5f-446e-80b7-4103a395415d",
   "metadata": {},
   "source": [
    "## calculations for NORMAL data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdb4a6f1-922b-4318-93b6-aec92dc8a062",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Iterate over the folders and files in the specified path\n",
    "basepath = '/Users/varyabazilova/Desktop/drivendata_rodeo/data/era5_land/'\n",
    "# basepath = '/Users/varyabazilova/Desktop/drivendata_rodeo/data/era5test/'\n",
    "\n",
    "folders = [d for d in os.listdir(basepath) if not d.startswith('.') and os.path.isdir(os.path.join(basepath, d))]\n",
    "\n",
    "for folder in folders:\n",
    "    ''' do the calcuations for all folders\n",
    "    every folder should consist the data with one variable'''\n",
    "    # identify folder\n",
    "    folder_path = os.path.join(basepath, folder)\n",
    "    # open and combine data\n",
    "    dataset = open_and_combine_nc_files(folder_path)\n",
    "    \n",
    "    print(dataset)\n",
    "    \n",
    "    #create a table \n",
    "    table = pd.DataFrame()\n",
    "    # gdf = gdf[:2]\n",
    "    \n",
    "    # iterate over every polygon in a shap file\n",
    "    for site_id in gdf.site_id:\n",
    "        \n",
    "        #take one polygon\n",
    "        shape = gdf[gdf.site_id == site_id]  \n",
    "        \n",
    "        #calculate stats \n",
    "        df1 = stat_by_shape(dataset, shape, kind='mean').rename(columns=lambda x: x + '_mean')        \n",
    "        df2 = stat_by_shape(dataset, shape, kind='min').rename(columns=lambda x: x + '_min')\n",
    "        df3 = stat_by_shape(dataset, shape, kind='max').rename(columns=lambda x: x + '_max')\n",
    "        df5 = stat_by_shape(dataset, shape, kind='std').rename(columns=lambda x: x + '_std')\n",
    "        # concat all togetehr\n",
    "        df = pd.concat([df1, df2, df3, df5], axis=1) \n",
    "        #add site-id of the polygon\n",
    "        df['site_id'] = shape['site_id'].iloc[0]\n",
    "        #append to the table\n",
    "        table = table.append(df)\n",
    "\n",
    "        # Save the table as a CSV file with the folder name\n",
    "    csv_filename = os.path.join(path_out, f'{folder}_table.csv')\n",
    "    table.to_csv(csv_filename, index=True) # important index=True!!! - this will save DATE \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a27eee7-3673-4d71-ace3-3b45c3ef0d7d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f741aa49-3697-4e69-a158-74d05b01a994",
   "metadata": {},
   "source": [
    "# res inflow - hardcoded downloads for ALL years "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f0028ec-9328-4a6c-b68b-61d838664551",
   "metadata": {},
   "outputs": [],
   "source": [
    "#  ----- what is going on here: -------\n",
    "\n",
    "# this is the resercoir inflow data, that is being pulled from rise api by BOR \n",
    "# the API request is hard coded through the web site: \n",
    "# the url is aquired through the time-series quiry \n",
    "# link: https://data.usbr.gov/time-series/search?id=10773.19550301-20231213&it=10773&lo=334\n",
    "# manually search for the name of the reserciir (select \"water\")\n",
    "# choose reservoir inflow, click time-series-queuery -> generate API request \n",
    "# NB! one request-link only has ONE page and is able to store only 2000values\n",
    "# => you need to iterate over PAGES in the link name to get more data\n",
    "\n",
    "\n",
    "# 'site_id' matched with the name of a reservoir name \n",
    "\n",
    "# 0              hungry_horse_reservoir_inflow     Hungry Horse \n",
    "# 2                    pueblo_reservoir_inflow     Pueblo\n",
    "# 8                           boise_r_nr_boise     Lucky Peak\n",
    "# 10              taylor_park_reservoir_inflow     Taylor Park\n",
    "# 12                    ruedi_reservoir_inflow     Ruedi\n",
    "# 13               fontenelle_reservoir_inflow     Fontenelle\n",
    "# 15     san_joaquin_river_millerton_reservoir     Folsom\n",
    "# 17                american_river_folsom_lake     Friant\n",
    "# 23                   boysen_reservoir_inflow     Boysen\n",
    "# 25                    owyhee_r_bl_owyhee_dam     Owyhee\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "cc250e0d-3b1c-478c-8394-da639ce75ee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import pandas as pd\n",
    "\n",
    "# ----- functions ------\n",
    "\n",
    "def process_inflow(json):\n",
    "    '''this function converts json data to pd df,\n",
    "    re-formats time the correct way\n",
    "    add site_id to the df'''\n",
    "    \n",
    "    df = pd.json_normalize(json['data'])\n",
    "    \n",
    "    # df = df[['attributes.dateTime', 'attributes.result']]\n",
    "    df[['date', 'time']] = df['attributes.dateTime'].str.split('T', expand=True)\n",
    "    \n",
    "    return df\n",
    "\n",
    "\n",
    "# --------- getting the data drom the api ----\n",
    "\n",
    "def fetch_data(url, num_pages, site_id):\n",
    "    ''' this fetches the data from the api:\n",
    "    '''\n",
    "    reservoir_data = pd.DataFrame()\n",
    "\n",
    "    for page_num in range(1, num_pages):\n",
    "        # Construct the URL with the current page number\n",
    "        current_url = f'{url}&page={page_num}'\n",
    "        \n",
    "        # Make the request and process the data\n",
    "        r = requests.get(current_url).json()\n",
    "        \n",
    "        #process inflow \n",
    "        df = pd.json_normalize(r['data'])\n",
    "        df = df[['attributes.dateTime', 'attributes.result']]\n",
    "        df[['date', 'time']] = df['attributes.dateTime'].str.split('T', expand=True)\n",
    "        \n",
    "        # df = process_inflow(r)\n",
    "        \n",
    "        # Append the data to the main dataframe\n",
    "        reservoir_data = reservoir_data.append(df)\n",
    "\n",
    "        # print(f'Processing page {page_num}')\n",
    "\n",
    "    # Add the site_id column\n",
    "    reservoir_data['site_id'] = site_id\n",
    "\n",
    "    return reservoir_data\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d4e321d1-e41d-431f-846f-d8ce996de174",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing page 1\n",
      "Processing page 2\n",
      "Processing page 3\n",
      "Processing page 4\n",
      "Processing page 5\n",
      "Processing page 6\n",
      "Processing page 7\n",
      "Processing page 8\n",
      "Processing page 9\n",
      "Processing page 10\n",
      "Processing page 11\n",
      "Processing page 12\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# order: folsom, boysen, taylor park, friant \n",
    "\n",
    "\n",
    "url_folsom = 'https://data.usbr.gov/rise/api/result?itemsPerPage=2000&order%5BdateTime%5D=ASC&itemId=10773&dateTime%5Bafter%5D=19550301&dateTime%5Bstrictly_before%5D=20231212' #folsom\n",
    "url_boysen = 'https://data.usbr.gov/rise/api/result?itemsPerPage=2000&order%5BdateTime%5D=ASC&itemId=197&dateTime%5Bafter%5D=19111001&dateTime%5Bstrictly_before%5D=20230915' #boysen\n",
    "url_TP = 'https://data.usbr.gov/rise/api/result?itemsPerPage=2000&order%5BdateTime%5D=ASC&itemId=794&dateTime%5Bafter%5D=19621001&dateTime%5Bstrictly_before%5D=20231213' #taylor park\n",
    "url_friant = 'https://data.usbr.gov/rise/api/result?itemsPerPage=2000&order%5BdateTime%5D=ASC&itemId=10806&dateTime%5Bafter%5D=19440222&dateTime%5Bstrictly_before%5D=20231212' #friant\n",
    "\n",
    "url_list= [url_folsom, url_boysen, url_TP, url_friant]\n",
    "num_pages_list = [13, 22, 12, 16]\n",
    "site_id_list = ['san_joaquin_river_millerton_reservoir', 'boysen_reservoir_inflow', 'taylor_park_reservoir_inflow', 'american_river_folsom_lake']\n",
    "\n",
    "\n",
    "master_table = pd.DataFrame()\n",
    "\n",
    "# Loop over the lists\n",
    "for url, num_pages, site_id in zip(url_list, num_pages_list, site_id_list):\n",
    "    # Call the fetch_data function\n",
    "    data = fetch_data(url, num_pages, site_id)\n",
    "    data = data[['attributes.result', 'date', 'site_id']]\n",
    "    master_table = master_table.append(data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bf6a0621-9f5e-4a00-b0c0-2eecd2904721",
   "metadata": {},
   "outputs": [],
   "source": [
    "master_table=master_table.rename(columns={'attributes.result': 'res_inflow'})#, 'B': 'Y'}, inplace=True)\n",
    "master_table = master_table.drop_duplicates()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a85f38f6-efa4-416d-be95-5ed1ebfb7b02",
   "metadata": {},
   "source": [
    "# make only test years "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "89f70bc4-5ffa-4ce4-8072-d855e06c46ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "master_table['date'] = pd.to_datetime(master_table['date'])\n",
    "master_table['year'] = master_table['date'].dt.year\n",
    "\n",
    "test_years =  [2005, 2007, 2009, 2011, 2013, 2015, 2017, 2019, 2021, 2023]\n",
    "\n",
    "master_table_test = master_table[master_table['year'].isin(test_years)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5a86ff9-ec0b-47ff-8174-07612d4a37bd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:rodeo] *",
   "language": "python",
   "name": "conda-env-rodeo-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}