Merge pull request #178 from aleaf/main

Update 11_xarray_mt_rainier...

Author: mnfienen

.github/workflows/test.yaml

@@ -76,3 +76,5 @@ jobs:
       shell: bash -l {0}
       run: |
         pytest tests/test_notebooks.py -v
+      env:
+        GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

notebooks/part0_python_intro/11_xarray_mt_rainier_precip.ipynb

@@ -66,7 +66,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Xarray"
+    "## Xarray basics"
    ]
   },
   {
@@ -163,7 +163,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Datasets\n",
+    "### Datasets\n",
     "A Dataset holds multiple DataArrays, analogous to a DataFrame in pandas housing multiple Series. Like in Pandas, the different DataArrays can share coordinates, and they can be assigned variable names (analogous to column names). The graphic at the beginning of this lesson shows a Dataset with precipitation, temperature, latitude and longitude DataArrays.\n",
     "\n",
     "A Dataset can be created from a DataArray with the `.to_dataset()` method:"
@@ -234,9 +234,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "\n",
-    "\n",
-    "### inputs\n",
+    "## The Mt. Rainer example\n",
+    " \n",
+    "### Inputs\n",
     "* [Daymet](https://daymet.ornl.gov/overview) gridded precipitation estimates for area around Mt. Rainier, 1980-2018. Obtained via the [Daymet Tile Selection Tool](https://daymet.ornl.gov/gridded/)\n",
     "* a [proj string](https://proj.org/usage/quickstart.html) for Daymet could be created from the [dataset description](https://daymet.ornl.gov/overview); however, it was easier to just copy it from [here](http://rpubs.com/tbiggs/GEOG576_EX3_DAYMET)\n",
     "* a 30 m DEM of elevations around Mt. Rainier, in geographic coordinates, created in the `gis_raster_mt_rainier_glaciers.ipynb` exercise"
@@ -329,7 +329,7 @@
    "source": [
     "### Let's add some topography for reference\n",
     "\n",
-    "Load the elevation data with `rasterio`, get the bounding extent, and `rasterio.show` to quickly plot the data with its coordinates (instead of row, column locations)"
+    "Load the elevation data into a Numpy array with `rasterio`, get the bounding extent, and `rasterio.show` to quickly plot the data with its coordinates (instead of row, column locations)"
    ]
   },
   {
@@ -355,7 +355,7 @@
     "* use the `matplotlib` `LightSource.hillshade` method to convert the elevations to a shaded relief raster\n",
     "    * assign `zorder=-1` to ensure that the data plot on the bottom\n",
     "* plot mean precip, but specify `alpha < 1` so that it has some transparency; the colorbar can be controlled with `cbar_kwargs` (arguments that are passed to the `matplotlib` colorbar constructor)\n",
-    "* set the plot limits to the extent of the elevation raster (`*` unpacks list elevations into individual arguments)"
+    "* set the plot limits to the extent of the elevation raster (`*` unpacks list elevations into individual arguments). This puts the elevation data in the correct geographic locations."
    ]
   },
   {
@@ -385,7 +385,7 @@
     "### Getting data at point locations\n",
     "\n",
     "Say we want to get data for the mountain summit, and at the Paradise Visitor Center in Mt. Rainier National Park.\n",
-    "We can easily get the lat, lon coordinates for these locations, but to get data from this dataset, we need to work in the [custom coordinate system for Daymet](https://daymet.ornl.gov/overview). Luckily, we found a PROJ string to define that, which we assigned to the variable `daymet_proj_string`. With the PROJ string, we can use `pyproj` to transform the coordinates to the Daymet CRS."
+    "We can easily get the lat, lon coordinates for these locations, but to get data from this dataset, we need to work in the [custom coordinate system for Daymet](https://daymet.ornl.gov/overview). Luckily, we found a PROJ string to define that, which we assigned to the variable `daymet_proj_string`. With the PROJ string, we can use `pyproj` to transform (reproject) the lat/lon coordinates to the Daymet CRS."
    ]
   },
   {
@@ -426,7 +426,7 @@
     "ax.set_xlim(*elev_extent[:2])\n",
     "for label, (y, x) in coords.items():\n",
     "    ax.scatter(x, y, c='k')\n",
-    "    ax.text(x, y, label, transform=ax.transData)"
+    "    ax.text(x+.003, y+.003, label, transform=ax.transData, c='k')"
    ]
   },
   {
@@ -566,7 +566,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "prcp_subset = ds.prcp.loc[:].where(valid_lat&valid_lon)\n",
+    "prcp_subset = ds.prcp.loc[:].where(valid_lat & valid_lon)\n",
     "prcp_subset = prcp_subset.dropna(dim='x', how='all').dropna(dim='y', how='all')\n",
     "prcp_subset.mean(axis=0).plot()"
    ]
@@ -593,10 +593,10 @@
     "output_folder.mkdir(exist_ok=True)\n",
     "\n",
     "encoding={'prcp': {'zlib': True, # compression algorithm\n",
-    "                     'complevel': 4, # compression level\n",
-    "                     'dtype': 'float32',\n",
-    "                     '_FillValue': -9999,\n",
-    "                     }}\n",
+    "                   'complevel': 4, # compression level\n",
+    "                   'dtype': 'float32',\n",
+    "                   '_FillValue': -9999,\n",
+    "                   }}\n",
     "\n",
     "prcp_subset.to_netcdf(output_folder / 'rainier_prcp_subset.nc',\n",
     "                     encoding=encoding)"
@@ -840,7 +840,7 @@
     "for i, site in enumerate(sites):\n",
     "    df = ds.prcp.sel(x=x[0], \n",
     "                     y=y[0],\n",
-    "                     method='nearest').drop(['lat', 'lon', 'x', 'y']).to_dataframe()\n",
+    "                     method='nearest').drop_vars(['lat', 'lon', 'x', 'y']).to_dataframe()\n",
     "    df.columns = [site]\n",
     "    dfs.append(df)\n",
     "\n",
@@ -937,9 +937,7 @@
    "metadata": {},
    "source": [
     "### Reprojecting the dataset\n",
-    "Here we reproject the dataset into geographic coordinates (Lat/Lon). We can use the spatial reference ID 4269 (an [EPSG code](https://www.linkedin.com/pulse/whats-epsg-how-do-i-use-salina-morrow#:~:text=The%20EPSG%20number%20is%20a,that%20can%20display%20those%20files.)) to specify the destination coordinate reference system (geographic coordinates with the NAD83 datum).\n",
-    "\n",
-    "Note that this dataset already includes `lat`/`lon` coordinates. Because the x and y axes are in a projected coordinate reference system (Lambert Conformal Conic), the `lat`/`lon` coordinates are two dimensional (each point in the dataset has its own `lat`/`lon` coordinate). **`rioxarray`** doesn't like this for some reason-- as of version 0.15, it produces a `ValueError: IndexVariable objects must be 1-dimensional` when we try to reproject. To get around this, we can drop the `lat`/`lon` coordinates before calling `.rio.project()`."
+    "Here we reproject the dataset into geographic coordinates (Lat/Lon). We can use the spatial reference ID 4269 (an [EPSG code](https://www.linkedin.com/pulse/whats-epsg-how-do-i-use-salina-morrow#:~:text=The%20EPSG%20number%20is%20a,that%20can%20display%20those%20files.)) to specify the destination coordinate reference system (geographic coordinates with the NAD83 datum)."
    ]
   },
   {
@@ -948,7 +946,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds_4269 = ds['prcp'].drop_vars(['lat', 'lon']).rio.reproject(4269)"
+    "ds_4269 = ds['prcp'].rio.reproject(4269)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that this dataset already includes `lat`/`lon` coordinates. Because the x and y axes are in a projected coordinate reference system (Lambert Conformal Conic), the `lat`/`lon` coordinates are two dimensional (each point in the dataset has its own `lat`/`lon` coordinate). Earlier versions of **`rioxarray`** didn't like this for some reason-- as of version 0.15, it produced a `ValueError: IndexVariable objects must be 1-dimensional` when calling `rio.reproject()`. Droping the `lat`/`lon` coordinates before calling `.rio.project()` was a way to get around this. For example: ```ds_4269 = ds['prcp'].drop_vars(['lat', 'lon']).rio.reproject(4269)```"
    ]
   },
   {
@@ -968,7 +973,7 @@
     "\n",
     "Write part of the dataset out to a GeoTIFF. \n",
     "\n",
-    "Earlier versions of rioxarray require us to remove the `'grid_mapping'` attribute first (at least in version 0.15, you will get an error otherwise). `grid_mapping` is a variable from the original NetCDF dataset that specifies that coordinate reference information is contained in a variable called `'lambert_conformal_conic'` (we used it above to specify the coordinate reference information for **rioxarray**).\n",
+    "Earlier versions of rioxarray required us to remove the `'grid_mapping'` attribute first (e.g. `del precip_slice.attrs['grid_mapping']`; at least in version 0.15, you would get an error otherwise). `grid_mapping` is a variable from the original NetCDF dataset that specifies that coordinate reference information is contained in a variable called `'lambert_conformal_conic'` (we used it above to specify the coordinate reference information for **rioxarray**).\n",
     "\n",
     "#### Write a single 2D timeslice to a raster"
    ]
@@ -979,10 +984,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "precip_slice = ds.sel(time='2018-12-30')['prcp']\n",
-    "\n",
-    "# remove grid_mapping if needed\n",
-    "#del precip_slice.attrs['grid_mapping']"
+    "precip_slice = ds.sel(time='2018-12-30')['prcp']"
    ]
   },
   {
@@ -1046,7 +1048,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "data = ds_4269.sel(time='2018')#.mean(axis=0)"
+    "data = ds_4269.sel(time='2018')"
    ]
   },
   {
@@ -1083,13 +1085,32 @@
     "gdf.plot(ax=ax)\n",
     "clipped.mean(axis=0).plot(ax=ax)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "pyclass",
    "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.12.11"
   }
  },
  "nbformat": 4,