Crop Yield Regional Interactive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview \n",
    "\n",
    "This script provides the water satisification requirement index (WSRI) over a given region for a specifc time. It uses the downloaded [Famine Land Data Assimilation System](https://ldas.gsfc.nasa.gov/FLDAS/) (FLDAS)* data which can be acquired through the Crop Yield Data Download script. \n",
    "\n",
    "\\**McNally, A., Arsenault, K., Kumar, S., Shukla, S., Peterson, P., Wang, S., Funk, C., Peters-Lidard, C.D., & Verdin, J. P. (2017). A land data assimilation system for sub-Saharan Africa food and water security applications. Scientific Data, 4, 170012*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User Tips\n",
    "\n",
    "- **Every cell ran when the notebook opened**\n",
    "\n",
    "\n",
    "- **Each input is \"live\"; entering an input then running the cell will delete your input**\n",
    "\n",
    "\n",
    "- **After making an input each cell below that input needs to be rerun to use that new input**\n",
    "\n",
    "\n",
    "- **If you get an error rerun the cells above the cell that errored and it will rectify**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0: The Dependencies\n",
    "\n",
    "This cell imports the python dependencies we use to download the data. It has already run. If you want to know which dependencies this code is using click Hide/Show Code. If you click the run cell button, the buttons for this cell will disappear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#importing libraries\n",
    "import glob\n",
    "from netCDF4 import Dataset\n",
    "import numpy as np\n",
    "import math\n",
    "import pcse\n",
    "import pandas as pd\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact\n",
    "import folium\n",
    "from folium.plugins import HeatMap\n",
    "import datetime\n",
    "import os\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", message=\"Cannot find a last shown plot to update.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Date selection\n",
    "\n",
    "Select the month you want to explore from part one the data download. \n",
    "\n",
    "\n",
    "*If you are using a binder instance of this notebook you may run out of memory if you select too many dates. This is indicated when the notebook disconnects in step 4.* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#choosing data file\n",
    "#grab dates from FLDAS files\n",
    "all_months = []\n",
    "\n",
    "#iterate through files\n",
    "for file in glob.glob(r'data/*.nc4'):\n",
    "    #print(file)\n",
    "    data = Dataset(file, 'r')\n",
    "    #time = data.variables['time'].units\n",
    "    #print(time)\n",
    "    date = file[26:30] + '-' + file[30:32]\n",
    "    #month = time[11:21]\n",
    "    all_months.append(date)\n",
    "    #print(time)\n",
    "\n",
    "#Widget for user to select desired months\n",
    "dates_data = widgets.SelectMultiple(\n",
    "    options = sorted(all_months),\n",
    "    description='Dates',\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "#interact function to store user input\n",
    "dates_selected = []\n",
    "def update(dates):\n",
    "    date_selected = list(dates)\n",
    "    return date_selected\n",
    "    \n",
    "#Retrieve elevevation data and stores in a dictionary of\n",
    "# Key = Long, Lat pair\n",
    "# Value = Elevation\n",
    "elev_dict = {}\n",
    "for file in glob.glob(r'data/*.csv'):\n",
    "    elev = pd.read_csv(file)\n",
    "    for idx, row in elev.iterrows():\n",
    "        elev_dict[(row[\"Longitude\"],row[\"Latitude\"])] = row[\"Elevation\"]\n",
    "\n",
    "dates_selection = interact(update, dates= dates_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Location Selection\n",
    "\n",
    "In the next step we determine the longitude and latitude window in which we want to draw data from and display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "date_selection = dates_selection.__dict__[\"widget\"].children[0].__dict__[\"_trait_values\"][\"value\"]\n",
    "\n",
    "if date_selection: \n",
    "    file_loc = r\"data/FLDAS_NOAH01_C_GL_M.A\" + date_selection[0][0:4]+date_selection[0][5:7] + \".001.nc.SUB.nc4\"\n",
    "    data_loc =Dataset(file_loc, 'r')\n",
    "    lat_array = data_loc.variables['Y'][:]#.compressed() # unmasks numpy_array lat stored in the FLDAS files \n",
    "    long_array = data_loc.variables['X'][:]#.compressed() # unmasks numpy array long stored in FLDAS files\n",
    "    min_lat = lat_array[0]\n",
    "    max_lat = lat_array[-1]\n",
    "    min_long = long_array[0]\n",
    "    max_long = long_array[-1]\n",
    "    long_points = {'globe':[], 'web':[]}\n",
    "    lat_points = {'globe':[], 'web':[]}\n",
    "    for i in range(len(lat_array)):\n",
    "        for j in range(len(long_array)):\n",
    "            if round(lat_array[i],2) not in lat_points['globe']:\n",
    "                lat_points['globe'].append(round(lat_array[i],2))\n",
    "            if round(long_array[j],2) not in long_points['globe']:\n",
    "                long_points['globe'].append(round(long_array[j],2))\n",
    "            #pt = transformer.transform(lat_array[i], long_array[j])\n",
    "            #if pt[0] not in long_points['web']: \n",
    "            #    long_points['web'].append(pt[0])\n",
    "            #if pt[1] not in lat_points['web']:\n",
    "            #    lat_points['web'].append(pt[1])\n",
    "\n",
    "    pts = [(min_lat, min_long), (max_lat, max_long)]\n",
    "    bbox = pts\n",
    "    #for pt in transformer.itransform(pts): \n",
    "    #    bbox.append(pt)       \n",
    "\n",
    "    #Plot is not the best as it has to replot the map each time instead of the just the box\n",
    "    #Need to find a better way\n",
    "    def update(min_latitude, max_latitude, min_longitude, max_longitude):\n",
    "        #Create the base figure\n",
    "        m = folium.Map(location=[(min_latitude + max_latitude) / 2, (min_longitude + max_longitude) / 2], zoom_start=5)\n",
    "        #add the map form the Bokeh map vendor in this case Stamen_Terrain --- see documentation\n",
    "        folium.Polygon(locations=[[(min_latitude, min_longitude),\n",
    "                               (max_latitude, min_longitude),\n",
    "                               (max_latitude, max_longitude),\n",
    "                               (min_latitude, max_longitude)]],\n",
    "                   color='red',\n",
    "                   fill_color='red',\n",
    "                   fill_opacity=0.5).add_to(m)\n",
    "\n",
    "        # Display the map\n",
    "        display(m)\n",
    "\n",
    "\n",
    "\n",
    "    loc_input = interact(update,min_latitude=widgets.FloatSlider(min=lat_points[\"globe\"][0], max=lat_points[\"globe\"][-1], value=lat_points[\"globe\"][0]), \n",
    "                         max_latitude=widgets.FloatSlider(min=lat_points[\"globe\"][0], max=lat_points[\"globe\"][-1], value=lat_points[\"globe\"][-1]),\n",
    "                         min_longitude=widgets.FloatSlider(min=long_points[\"globe\"][0],max=long_points[\"globe\"][-1],value= long_points[\"globe\"][0]),\n",
    "                         max_longitude=widgets.FloatSlider(min=long_points[\"globe\"][0],max=long_points[\"globe\"][-1], value=long_points[\"globe\"][-1]))\n",
    "\n",
    "else: \n",
    "    print (\"Waiting for Inputs\")\n",
    "print (\"  \") #Hides the output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Getting the Elevation\n",
    "\n",
    "This cell retrieves the desired regional map and obtains the elevations for that region. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if date_selection: \n",
    "    #Take user inputs and get new arrays from lat and long\n",
    "    lat_min = round(loc_input.__dict__[\"widget\"].children[0].__dict__[\"_trait_values\"][\"value\"],2)\n",
    "    lat_max = round(loc_input.__dict__[\"widget\"].children[1].__dict__[\"_trait_values\"][\"value\"], 2)\n",
    "    lon_min = round(loc_input.__dict__[\"widget\"].children[2].__dict__[\"_trait_values\"][\"value\"], 2)\n",
    "    lon_max = round(loc_input.__dict__[\"widget\"].children[3].__dict__[\"_trait_values\"][\"value\"], 2)\n",
    "\n",
    "    #Make user input and masked array compatabile and updated array to user inputs\n",
    "    lat=[]\n",
    "    lons =[]\n",
    "    lats_user =[]\n",
    "    lons_user = []\n",
    "    lats = np.round(data.variables['Y'][:],2)\n",
    "    lons = np.round(data.variables['X'][:],2)\n",
    "    lat_min_idx = np.where(lats.data==lat_min)\n",
    "    lat_max_idx = np.where(lats.data==lat_max)\n",
    "    lats_user = lats[lat_min_idx[0][0]:lat_max_idx[0][0]]\n",
    "    lon_min_idx = np.where(lons.data==lon_min)\n",
    "    lon_max_idx = np.where(lons.data==lon_max)\n",
    "    lons_user = lons[lon_min_idx[0][0]:lon_max_idx[0][0]]\n",
    "    #get web projections for use in heatmap\n",
    "    lats_user_globe = lat_points[\"globe\"][lat_min_idx[0][0]:lat_max_idx[0][0]]\n",
    "    lons_user_globe =long_points[\"globe\"][lon_min_idx[0][0]:lon_max_idx[0][0]]\n",
    "\n",
    "\n",
    "    #convert into a grid of lat/long coordinates\n",
    "    lon2, lat2 = np.meshgrid(lons_user, lats_user)\n",
    "    lat_flat = lat2.flatten()\n",
    "    lon_flat = lon2.flatten()\n",
    "    #lon_web,lat_web =np.meshgrid(lons_user_globe, lats_user_globe)\n",
    "    #lon_web_flat = lon_web.flatten()\n",
    "    #lat_web_flat = lat_web.flatten()\n",
    "    area = len(lat_flat)\n",
    "    #reqs = len(lat_flat)*len(lon_flat)\n",
    "    elevation = np.zeros(area)\n",
    "\n",
    "    for x in range(area):\n",
    "       elevation[x] = elev_dict[(lon_flat[x],lat_flat[x])]\n",
    "\n",
    "    print (\"This cell has updated the bounding box based on user input and retrieved the elevation for that region.\")\n",
    "else:\n",
    "    print(\"Waiting for inputs.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Water Requirement ans Satsification Index Calculations\n",
    "\n",
    "The next cell takes the elevation and FLDAS data and uses the [Python Crop Simulation Environment](https://pcse.readthedocs.io/en/stable/), specifically the Penman-Monteith algorithm, to calculate the Water Requirement Satsification Index (WRSI). Each Month is stored in a table (or dataframe). A table for one month shows as an example. It shows the latitude, longitude, WRSI and mercator equivalents on the x an y axis for longitude and latitude.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "if date_selection: \n",
    "    status_dict = {}\n",
    "    print (\"Calculating.....\")\n",
    "    for date in date_selection: \n",
    "        file_loc = r\"data/FLDAS_NOAH01_C_GL_M.A\" + date[0:4]+date[5:7] + \".001.nc.SUB.nc4\"\n",
    "        data =Dataset(file_loc, 'r')\n",
    "        #lat_array = data_loc.variables['Y'][:]#.compressed() # unmasks numpy_array lat stored in the FLDAS files \n",
    "        #long_array = data_loc.variables['X'][:]#.compressed() # unmasks numpy array long stored in FLDAS files    \n",
    "        #defining the date\n",
    "        month_data = datetime.date(int(date[0:4]),int(date[5:8]),1)\n",
    "        #WR = np.arange(area).reshape(len(lats_user),(len(lons_user)))\n",
    "        WR = np.zeros(area)\n",
    "        WRSI = np.zeros(area)\n",
    "\n",
    "        #Get the subset sets of FLDAS downloaded data from user input\n",
    "        air_temp = data.variables['Tair_f_tavg'][0][lat_min_idx[0][0]:lat_max_idx[0][0],lon_min_idx[0][0]:lon_max_idx[0][0]]\n",
    "        humidity = data.variables['Qair_f_tavg'][0][lat_min_idx[0][0]:lat_max_idx[0][0],lon_min_idx[0][0]:lon_max_idx[0][0]]\n",
    "        net_short_radiation = data.variables['SWdown_f_tavg'][0][lat_min_idx[0][0]:lat_max_idx[0][0],lon_min_idx[0][0]:lon_max_idx[0][0]]\n",
    "        net_long_radiation = data.variables['Lwnet_tavg'][0][lat_min_idx[0][0]:lat_max_idx[0][0],lon_min_idx[0][0]:lon_max_idx[0][0]]\n",
    "        wind_speed = data.variables['Wind_f_tavg'][0][lat_min_idx[0][0]:lat_max_idx[0][0],lon_min_idx[0][0]:lon_max_idx[0][0]]\n",
    "        evapotranspiration = data.variables['Evap_tavg'][0][lat_min_idx[0][0]:lat_max_idx[0][0],lon_min_idx[0][0]:lon_max_idx[0][0]]\n",
    "\n",
    "        #get locations in list\n",
    "        locs = tuple(zip(lat_flat, lon_flat))\n",
    "\n",
    "        #grabbing air temp, humidity, short radiation,, long radiation, wind speeds, and evapotranspiration data from the FLDAS dataset\n",
    "        # unit conversions for the PM equation\n",
    "        # converts kg/m^2/s to mm/day\n",
    "        humidity_mm = humidity * 86400\n",
    "\n",
    "        #convert wind speed from 10m height to 2m height\n",
    "        wind_speed_2 = wind_speed * 4.87 / math.log ( 67.8 * 10 - 5.42 ) \n",
    "\n",
    "        #convert air temp from kelvin to celcius\n",
    "        air_temp_C = air_temp -273.15\n",
    "        #flatten all arrays\n",
    "        air_temp_flat = air_temp_C.flatten()\n",
    "        wind_speed_flat = wind_speed_2.flatten()\n",
    "        net_short_radiation_flat =  net_short_radiation.flatten()\n",
    "        evapotranspiration_flat = evapotranspiration.flatten()\n",
    "        evapotranspiration_flat = evapotranspiration_flat * 86400\n",
    "        #Calculate the WR for all co-ordinates in the window\n",
    "        WRSI_dict = {\"latitude\":[], \"longitude\":[], \"WRSI\":[]}# \"mercator_x\":[], \"mercator_y\":[]}\n",
    "        for x in range(area):\n",
    "            WR[x] = pcse.util.penman_monteith(month_data, lat_flat[x], elevation[x], air_temp_flat[x] , air_temp_flat[x] , \n",
    "                                              net_short_radiation_flat[x], 14, wind_speed_flat[x])\n",
    "\n",
    "            WRSI_dict[\"latitude\"].append(lat_flat[x])\n",
    "            WRSI_dict[\"WRSI\"].append(evapotranspiration_flat[x]*100/WR[x])\n",
    "            WRSI_dict[\"longitude\"].append(lon_flat[x])\n",
    "            #WRSI_dict[\"mercator_x\"].append(lon_web_flat[x])\n",
    "            #WRSI_dict[\"mercator_y\"].append(lat_web_flat[x])\n",
    "\n",
    "\n",
    "        status_dict[str(month_data)[0:7]] = pd.DataFrame(WRSI_dict)\n",
    "        print (\"Water Satisifaction Requirement Index calculation for \", str(month_data)[0:7], \"complete.\")\n",
    "    print(\"Calculations Complete\")\n",
    "    #Displaye the first table\n",
    "    keys = list(status_dict.keys())\n",
    "    status_dict[keys[0]]\n",
    "else: \n",
    "    print(\"Waiting for Inputs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Plotting the Regional Data\n",
    "\n",
    "This cell creates a map and plots the WRSI data on it to form a heatmap for the longitude and latitude coordinates. Users can select the month they would like to see. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if date_selection: \n",
    "    bbox2 = [(lat_min, lon_min), (lat_max, lon_max)]\n",
    "    #bbox2 = []\n",
    "    #for pt in transformer.itransform(min_max_pts): \n",
    "    #    bbox2.append(pt)   \n",
    "        \n",
    "    month =widgets.Dropdown(\n",
    "        options=date_selection,\n",
    "        value=date_selection[0],\n",
    "        # rows=10,\n",
    "        description=\"Plot Month\",\n",
    "        disabled=False)\n",
    "\n",
    "    size = widgets.IntText(\n",
    "        value=600,\n",
    "        description='Plot Size',\n",
    "        disabled=False\n",
    "    )\n",
    "\n",
    "    title = widgets.Text(\n",
    "        value='Regional Maximal Crop Yield',\n",
    "        description='Title',\n",
    "        disabled=False\n",
    "    )\n",
    "\n",
    "    #colors = list(RdYlGn[11])   \n",
    "    #colors.reverse()\n",
    "    #mapper = LinearColorMapper(palette=colors, low=status_dict[month.value].WRSI.min(), high=status_dict[month.value].WRSI.max())\n",
    "\n",
    "    #color_bar = ColorBar(color_mapper=mapper, location=(0, 0),\n",
    "    #                 ticker=BasicTicker(desired_num_ticks=len(colors)),\n",
    "    #                 formatter=PrintfTickFormatter(format=\"%d%%\"))\n",
    "\n",
    "    def plot_heatmap(month, size, title):\n",
    "        df1 = status_dict[month]\n",
    "        df = df1.dropna()\n",
    "        m = folium.Map(location=[df['latitude'].mean(), df['longitude'].mean()], zoom_start=5)\n",
    "        \n",
    "        \n",
    "        heat_data = [[row['latitude'], row['longitude'], row['WRSI']] for index, row in df.iterrows()]\n",
    "        #print(heat_data)\n",
    "        HeatMap(heat_data).add_to(m)\n",
    "        \n",
    "        #m.save(title + '.html')\n",
    "        display(m)\n",
    "\n",
    "    interact(plot_heatmap, month=month, size=size, title=title)\n",
    "else: \n",
    "    print(\"Waiting for Inputs\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Save the WRSI Data\n",
    "\n",
    "If you would like to save the WSRI data. Enter the below inputs and each table (dataframe) of location and WSRI will be saved as a .csv (comma seperated value) file (which can be opened with any spreadsheet program or read into a ABM). \n",
    "\n",
    "If you do not want to save the dataframes, do not enter anything in the next cell. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hide_code()\n",
    "run_code()\n",
    "\n",
    "#get filepath of user\n",
    "filepath = widgets.Text(\n",
    "    value='Enter File Path',\n",
    "    placeholder='Type something',\n",
    "    description='Filepath:',\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "def save_files(filepath):\n",
    "    if filepath == 'Enter File Path':\n",
    "        pass\n",
    "    else: \n",
    "        for k,v in status_dict.items(): \n",
    "            #get file name\n",
    "            filename = \"WRSI_\"+ k +\".csv\"\n",
    "            #combine filepath with name\n",
    "            file = os.path.join(filepath,filename)\n",
    "            #save pandas dataframe as .csv\n",
    "            v.to_csv(file)\n",
    "            print(\"{} has been saved.\".format(filename))\n",
    "\n",
    "file_complete = interact(save_files, filepath=filepath)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.9.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}