add: find nearest grid cell, show selected grid cell in map

415e3f10 · Marco Kulüke · d37a8f0b · 415e3f10
Commit 415e3f10 authored 4 years ago by Marco Kulüke
--- a/notebooks/use-case_count_summer_days.ipynb
+++ b/notebooks/use-case_count_summer_days.ipynb
@@ -11,12 +11,8 @@
   "cell_type": "markdown",
   "metadata": {},
   "source": [
-    "Workflow:\n",
-    "- chose Location either by name (https://pypi.org/project/geopy/) or coordniate (drop-down menu)\n",
-    "- get confirmation via map (https://github.com/python-visualization/folium)\n",
-    "- choose year by drop-down menu\n",
-    "- get time series with 25 threshold shown\n",
-    "- get number of heatdays\n",
+    "TBD:\n",
+    "- choose year by drop-down menu and select correspondig dataset\n",
    "- future:\n",
    "    - compare two places\n",
    "    - choose different index"
@@ -28,8 +24,6 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "# Load required packages\n",
-    "\n",
    "import intake\n",
    "import requests\n",
    "import folium\n",
@@ -37,7 +31,8 @@
    "import matplotlib\n",
    "from IPython.display import display\n",
    "from ipywidgets import widgets\n",
-    "from geopy.geocoders import Nominatim"
+    "from geopy.geocoders import Nominatim\n",
+    "import numpy as np"
   ]
  },
  {
@@ -75,7 +70,7 @@
    "location = geolocator.geocode(place.value)\n",
    "print(location.address)\n",
    "print((location.latitude, location.longitude))\n",
-    "print(location.raw)"
+    "#print(location.raw)"
   ]
  },
  {
@@ -110,15 +105,8 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "x = range(1850, 2020)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "x = range(1850, 2020) # TBD: needs to be historical time frame\n",
+    "\n",
    "year = widgets.Dropdown(\n",
    "    options=x,\n",
    "    description='Select year: ',\n",
@@ -176,7 +164,7 @@
   "cell_type": "markdown",
   "metadata": {},
   "source": [
-    "## Browse in Intake Catalogue"
+    "## Browse Intake Catalogue and open Variable"
   ]
  },
  {
@@ -193,13 +181,20 @@
    "cat = col.search(**query)"
   ]
  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "show query"
+   ]
+  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
-    "cat.df#['path'][0]"
+    "cat.df"
   ]
  },
  {
@@ -208,14 +203,16 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "ds_tasmax = xr.open_dataset(cat.df['path'][0])"
+    "ds_tasmax = xr.open_dataset(cat.df['path'][0]) # TBD: open data set of selected year\n",
+    "\n",
+    "ds_tasmax['tasmax']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
-    "## Open Variable and Plot"
+    "## Compare Model Grid Cell with chosen Location"
   ]
  },
  {
@@ -224,7 +221,18 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "ds_tasmax['tasmax']"
+    "# Find nearest model coordinate\n",
+    "\n",
+    "# First, find the index of the grid point nearest a specific lat/lon.   \n",
+    "abslat = np.abs(ds_tasmax['lat']-location.latitude)\n",
+    "abslon = np.abs(ds_tasmax['lon']-location.longitude)\n",
+    "c = np.maximum(abslon, abslat)\n",
+    "\n",
+    "([xloc], [yloc]) = np.where(c == np.min(c))\n",
+    "\n",
+    "# Now I can use that index location to get the values at the x/y diminsion\n",
+    "\n",
+    "tas_max_hamburg = ds_tasmax['tasmax'][0, yloc, xloc] -273.15"
   ]
  },
  {
@@ -232,13 +240,7 @@
   "execution_count": null,
   "metadata": {},
   "outputs": [],
-   "source": [
-    "lat=ds_tasmax['lat']\n",
-    "lon=ds_tasmax['lon']\n",
-    "\n",
-    "mask_lat = (lat > 40) & (lat < 50)\n",
-    "mask_lon = (lon > 0) & (lon < 10)"
-   ]
+   "source": []
  },
  {
   "cell_type": "code",
@@ -246,16 +248,27 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "mask_lat[0]"
+    "m = folium.Map(location=[location.latitude, location.longitude], zoom_start=8)\n",
+    "tooltip = location.latitude, location.longitude\n",
+    "folium.Marker([ds_tasmax['lat'][yloc], ds_tasmax['lon'][xloc]], tooltip=tooltip, popup='Model Grid Cell Center').add_to(m)\n",
+    "folium.Marker([location.latitude, location.longitude], tooltip=tooltip, popup='Location selected by User').add_to(m)\n",
+    "\n",
+    "# Define rectangle\n",
+    "rect_lat1 = (ds_tasmax['lat'][yloc -1] +ds_tasmax['lat'][yloc])/2\n",
+    "rect_lon1 = (ds_tasmax['lon'][xloc -1] +ds_tasmax['lon'][xloc])/2\n",
+    "rect_lat2 = (ds_tasmax['lat'][yloc +1] +ds_tasmax['lat'][yloc])/2\n",
+    "rect_lon2 = (ds_tasmax['lon'][xloc +1] +ds_tasmax['lon'][xloc])/2\n",
+    "#\n",
+    "\n",
+    "folium.Rectangle(bounds=[[rect_lat1,rect_lon1], [rect_lat2,rect_lon2]], color='#ff7800', fill=True, fill_color='#ffff00', fill_opacity=0.2).add_to(m)\n",
+    "m"
   ]
  },
  {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
   "metadata": {},
-   "outputs": [],
   "source": [
-    "tas_max_hamburg = ds_tasmax['tasmax'][:, mask_lat, mask_lon][:, 0, 0] -273.15"
+    "## Draw Temperature Time Series and count Summer days"
   ]
  },
  {
@@ -264,6 +277,9 @@
   "metadata": {},
   "outputs": [],
   "source": [
+    "# Now I can use that index location to get the values at the lat/lon location\n",
+    "\n",
+    "tas_max_hamburg = ds_tasmax['tasmax'][:, yloc, xloc] -273.15\n",
    "matplotlib.pyplot.plot(tas_max_hamburg)"
   ]
  },

 %% Cell type:markdown id: tags:

 # Count Annual Summer Days for a particular Location

 %% Cell type:markdown id: tags:

-Workflow:
- chose Location either by name (https://pypi.org/project/geopy/) or coordniate (drop-down menu)
- get confirmation via map (https://github.com/python-visualization/folium)
- choose year by drop-down menu
- get time series with 25 threshold shown
- get number of heatdays
+TBD:
+- choose year by drop-down menu and select correspondig dataset
 - future:
    - compare two places
    - choose different index

 %% Cell type:code id: tags:

 ``` python
-# Load required packages
-
 import intake
 import requests
 import folium
 import xarray as xr
 import matplotlib
 from IPython.display import display
 from ipywidgets import widgets
 from geopy.geocoders import Nominatim
+import numpy as np
 ```

 %% Cell type:markdown id: tags:

 ## Choose Location
 If ambiguous the more likely location will be chosen

 %% Cell type:code id: tags:

 ``` python
 place = widgets.Text()
 display(place)
 ```

 %% Cell type:markdown id: tags:

 Find Place by using GeoPy

 %% Cell type:code id: tags:

 ``` python
 geolocator = Nominatim(user_agent="any_agent")
 location = geolocator.geocode(place.value)
 print(location.address)
 print((location.latitude, location.longitude))
-print(location.raw)
+#print(location.raw)
 ```

 %% Cell type:markdown id: tags:

 ## Draw Map

 %% Cell type:code id: tags:

 ``` python
 m = folium.Map(location=[location.latitude, location.longitude])
 tooltip = location.latitude, location.longitude
 folium.Marker([location.latitude, location.longitude], tooltip=tooltip).add_to(m)
 m
 ```

 %% Cell type:markdown id: tags:

 ## Choose Year

 %% Cell type:code id: tags:

 ``` python
-x = range(1850, 2020)
-```
-
-%% Cell type:code id: tags:
+x = range(1850, 2020) # TBD: needs to be historical time frame

-``` python
 year = widgets.Dropdown(
    options=x,
    description='Select year: ',
    disabled=False,
 )
 year
 ```

 %% Cell type:code id: tags:

 ``` python
 year.value
 ```

 %% Cell type:markdown id: tags:

 ## Load Intake Catalogue
 The catalog files are uploaded to the cloud. There is a gsv.gz file and a json file under
 https://swiftbrowser.dkrz.de/public/dkrz_a44962e3ba914c309a7421573a6949a6/intake-esm/
 We load the catalog descriptor with intake:

 %% Cell type:code id: tags:

 ``` python
 col_url = "https://swift.dkrz.de/v1/dkrz_a44962e3ba914c309a7421573a6949a6/intake-esm/mistral-cmip6.json"
 col = intake.open_esm_datastore(col_url)
 ```

 %% Cell type:markdown id: tags:

 Let's see what is in the intake catalog. The underlying data base is given as a panda dataframe which we can access with 'col.df'. col.df.head() shows us the first rows of the table:

 %% Cell type:code id: tags:

 ``` python
 col.df.head()
 ```

 %% Cell type:markdown id: tags:

-## Browse in Intake Catalogue
+## Browse Intake Catalogue and open Variable

 %% Cell type:code id: tags:

 ``` python
 query = dict(
    source_id = 'MPI-ESM1-2-HR',
    variable_id="tasmax",
    table_id="day",
    experiment_id=["historical"])
 cat = col.search(**query)
 ```

+%% Cell type:markdown id: tags:
+
+show query
+
 %% Cell type:code id: tags:

 ``` python
-cat.df#['path'][0]
+cat.df
 ```

 %% Cell type:code id: tags:

 ``` python
-ds_tasmax = xr.open_dataset(cat.df['path'][0])
+ds_tasmax = xr.open_dataset(cat.df['path'][0]) # TBD: open data set of selected year
+
+ds_tasmax['tasmax']
 ```

 %% Cell type:markdown id: tags:

-## Open Variable and Plot
+## Compare Model Grid Cell with chosen Location

 %% Cell type:code id: tags:

 ``` python
-ds_tasmax['tasmax']
-```
+# Find nearest model coordinate

-%% Cell type:code id: tags:
+# First, find the index of the grid point nearest a specific lat/lon.
+abslat = np.abs(ds_tasmax['lat']-location.latitude)
+abslon = np.abs(ds_tasmax['lon']-location.longitude)
+c = np.maximum(abslon, abslat)

-``` python
-lat=ds_tasmax['lat']
-lon=ds_tasmax['lon']
+([xloc], [yloc]) = np.where(c == np.min(c))
+
+# Now I can use that index location to get the values at the x/y diminsion

-mask_lat = (lat > 40) & (lat < 50)
-mask_lon = (lon > 0) & (lon < 10)
+tas_max_hamburg = ds_tasmax['tasmax'][0, yloc, xloc] -273.15
 ```

 %% Cell type:code id: tags:

 ``` python
-mask_lat[0]
 ```

 %% Cell type:code id: tags:

 ``` python
-tas_max_hamburg = ds_tasmax['tasmax'][:, mask_lat, mask_lon][:, 0, 0] -273.15
+m = folium.Map(location=[location.latitude, location.longitude], zoom_start=8)
+tooltip = location.latitude, location.longitude
+folium.Marker([ds_tasmax['lat'][yloc], ds_tasmax['lon'][xloc]], tooltip=tooltip, popup='Model Grid Cell Center').add_to(m)
+folium.Marker([location.latitude, location.longitude], tooltip=tooltip, popup='Location selected by User').add_to(m)
+
+# Define rectangle
+rect_lat1 = (ds_tasmax['lat'][yloc -1] +ds_tasmax['lat'][yloc])/2
+rect_lon1 = (ds_tasmax['lon'][xloc -1] +ds_tasmax['lon'][xloc])/2
+rect_lat2 = (ds_tasmax['lat'][yloc +1] +ds_tasmax['lat'][yloc])/2
+rect_lon2 = (ds_tasmax['lon'][xloc +1] +ds_tasmax['lon'][xloc])/2
+#
+
+folium.Rectangle(bounds=[[rect_lat1,rect_lon1], [rect_lat2,rect_lon2]], color='#ff7800', fill=True, fill_color='#ffff00', fill_opacity=0.2).add_to(m)
+m
 ```

+%% Cell type:markdown id: tags:
+
+## Draw Temperature Time Series and count Summer days
+
 %% Cell type:code id: tags:

 ``` python
+# Now I can use that index location to get the values at the lat/lon location
+
+tas_max_hamburg = ds_tasmax['tasmax'][:, yloc, xloc] -273.15
 matplotlib.pyplot.plot(tas_max_hamburg)
 ```

 %% Cell type:code id: tags:

 ``` python
 tas_max_hamburg[tas_max_hamburg > 25]
 ```