New tutorial

soilwater · Apr 11, 2024 · 894ca63 · 894ca63
1 parent b2f78fd
commit 894ca63
Show file tree

Hide file tree

Showing 14 changed files with 3,858 additions and 17 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/_quarto.yml b/_quarto.yml
@@ -105,6 +105,8 @@ book:
         - exercises/newton_law_cooling.ipynb
         - exercises/predator_prey_model.ipynb
         - exercises/reaction_diffusion.ipynb
+        - exercises/green_ampt.ipynb
+
 
     - part: "STOCHASTIC MODELS"
       chapters:
@@ -126,7 +128,7 @@ book:
         - exercises/yield_monitor_zones.ipynb
         - exercises/largest_empty_circle.ipynb
         - exercises/weather_network_coverage.ipynb
-        - exercises/raster_files.ipynb
+        - exercises/raster_and_vector.ipynb
         - exercises/inverse_distance_weighting.ipynb
         - exercises/builtin_spatial_interpolators.ipynb
         - exercises/semivariogram.ipynb

diff --git a/datasets/spatial/konza_watersheds/2024-04-03_Sentinel-2_L2A_B03_B04_B08-16bit.tiff b/datasets/spatial/konza_watersheds/2024-04-03_Sentinel-2_L2A_B03_B04_B08-16bit.tiff
diff --git a/datasets/spatial/konza_watersheds/konza_watersheds.geojson b/datasets/spatial/konza_watersheds/konza_watersheds.geojson
diff --git a/docs/exercises/raster_and_vector.html b/docs/exercises/raster_and_vector.html
diff --git a/docs/exercises/raster_and_vector_files/figure-html/cell-10-output-1.png b/docs/exercises/raster_and_vector_files/figure-html/cell-10-output-1.png
diff --git a/docs/exercises/raster_and_vector_files/figure-html/cell-14-output-1.png b/docs/exercises/raster_and_vector_files/figure-html/cell-14-output-1.png
diff --git a/docs/exercises/raster_and_vector_files/figure-html/cell-17-output-1.png b/docs/exercises/raster_and_vector_files/figure-html/cell-17-output-1.png
diff --git a/docs/exercises/raster_and_vector_files/figure-html/cell-19-output-1.png b/docs/exercises/raster_and_vector_files/figure-html/cell-19-output-1.png
diff --git a/docs/exercises/raster_and_vector_files/figure-html/cell-6-output-1.png b/docs/exercises/raster_and_vector_files/figure-html/cell-6-output-1.png
diff --git a/docs/search.json b/docs/search.json
@@ -2672,5 +2672,68 @@
     "title": "90  Estimate daily vapor pressure",
     "section": "Practice",
     "text": "Practice\n\nCompute the mean absolute error and the mean bias error for the different methods.\nLoad the 5-minute dataset for Ashland Bottoms, KS and compute e_sat, e_act, and vpd with the fitted linear model and the machine learning model. Does it work for another locations without any additional calibration?"
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#read-vector-layer",
+    "href": "exercises/raster_and_vector.html#read-vector-layer",
+    "title": "86  Raster and vector",
+    "section": "Read vector layer",
+    "text": "Read vector layer\n\n# Read file with watersheds\nvector_path = '../datasets/spatial/konza_watersheds/konza_watersheds.geojson'\ngdf = gpd.read_file(vector_path)\ngdf.head(3)\n\n\n\n\n\n\n\n\nCODE_1\nNAME_1\nAREA\nPERIMETER\nACRES\nHECTARES_1\nDATAID_1\nDATACODE_2\ngeometry\n\n\n\n\n0\n0R1A\nR1A\n484815.196184\n2939.146533\n119.800444\n48.481520\nGIS0321\nGIS032\nPOLYGON ((707045.277 4327611.469, 707038.577 4...\n\n\n1\nR20A\nR20A\n263364.016292\n2509.563160\n65.078666\n26.336402\nGIS0322\nGIS032\nPOLYGON ((707080.980 4327453.175, 707080.987 4...\n\n\n2\n002A\n2A\n283017.160279\n2417.724921\n69.935063\n28.301716\nGIS0323\nGIS032\nPOLYGON ((707362.870 4327681.792, 707364.460 4...\n\n\n\n\n\n\n\n\n# Display coordiante reference system of vector layer\nprint(gdf.crs)\n\nEPSG:26914\n\n\nThe first step to perform spatial operations between vector and raster layers is to have both data sources in the same coordinate reference system. Our vector layer has geometries in projected coordinates (meters), but the raster datasets will be in geographic coordinates.\nTo perform operations that involve distance computations is best to have both datasets in projected coordiantes. To perform other operations like clipping layers, then both geographic and projected coordinates should perform similarly. In this example we will convert the vector layer to geographic coordinates.\n\n# Change coordinate reference system from projected to geographic\ngdf.to_crs(epsg=4326, inplace=True)\n\n\n# Visualize vector map of watersheds\ngdf.plot(facecolor='None', edgecolor='k')\nplt.title('Konza Watersheds')\nplt.xlabel('Longitude')\nplt.ylabel('Latitude')\nplt.show()"
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#read-raster-layer",
+    "href": "exercises/raster_and_vector.html#read-raster-layer",
+    "title": "86  Raster and vector",
+    "section": "Read raster layer",
+    "text": "Read raster layer\nThe raster layer is a multi-band geo-referenced image at 10-meter spatial resolution obtained from Sentinel-2 satellite. Image data type is unsigned 16-bit (2^{16}). The bands are as follows:\n\nBand 1: Senteninel B03 (green, 560 nm)\nBand 2: Sentinel B04 (red, 665 nm)\nBand 3: Sentinel B08 (near infrared, 842 nm)\n\n\n# Load raster layer\nraster = xr.open_dataarray('../datasets/spatial/konza_watersheds/2024-04-03_Sentinel-2_L2A_B03_B04_B08-16bit.tiff')\nraster\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n&lt;xarray.DataArray 'band_data' (band: 3, y: 1060, x: 1461)&gt;\n[4645980 values with dtype=float32]\nCoordinates:\n  * band         (band) int64 1 2 3\n  * x            (x) float64 -96.66 -96.66 -96.66 ... -96.53 -96.53 -96.53\n  * y            (y) float64 39.15 39.15 39.15 39.15 ... 39.06 39.06 39.06 39.06\n    spatial_ref  int64 ...\nAttributes:\n    AREA_OR_POINT:           Area\n    TIFFTAG_RESOLUTIONUNIT:  1 (unitless)\n    TIFFTAG_XRESOLUTION:     1\n    TIFFTAG_YRESOLUTION:     1xarray.DataArray'band_data'band: 3y: 1060x: 1461...[4645980 values with dtype=float32]Coordinates: (4)band(band)int641 2 3array([1, 2, 3])x(x)float64-96.66 -96.66 ... -96.53 -96.53array([-96.664536, -96.664446, -96.664356, ..., -96.533314, -96.533224,\n       -96.533134])y(y)float6439.15 39.15 39.15 ... 39.06 39.06array([39.154246, 39.154156, 39.154066, ..., 39.059142, 39.059052, 39.058962])spatial_ref()int64...crs_wkt :GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]semi_major_axis :6378137.0semi_minor_axis :6356752.314245179inverse_flattening :298.257223563reference_ellipsoid_name :WGS 84longitude_of_prime_meridian :0.0prime_meridian_name :Greenwichgeographic_crs_name :WGS 84horizontal_datum_name :World Geodetic System 1984grid_mapping_name :latitude_longitudespatial_ref :GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]GeoTransform :-96.664581 9.000136892538973e-05 0.0 39.154291 0.0 -8.997547169811285e-05[1 values with dtype=int64]Indexes: (3)bandPandasIndexPandasIndex(Index([1, 2, 3], dtype='int64', name='band'))xPandasIndexPandasIndex(Index([-96.66453599931553,  -96.6644459979466, -96.66435599657768,\n       -96.66426599520875, -96.66417599383983,  -96.6640859924709,\n       -96.66399599110198, -96.66390598973305, -96.66381598836414,\n        -96.6637259869952,\n       ...\n        -96.5339440130048, -96.53385401163587, -96.53376401026694,\n       -96.53367400889802, -96.53358400752909, -96.53349400616017,\n       -96.53340400479124, -96.53331400342232, -96.53322400205339,\n       -96.53313400068447],\n      dtype='float64', name='x', length=1461))yPandasIndexPandasIndex(Index([ 39.15424601226415,  39.15415603679246, 39.154066061320755,\n        39.15397608584906,  39.15388611037736,  39.15379613490566,\n       39.153706159433966, 39.153616183962264,  39.15352620849057,\n        39.15343623301887,\n       ...\n       39.059771766981136,  39.05968179150943,  39.05959181603774,\n        39.05950184056604,  39.05941186509434, 39.059321889622645,\n        39.05923191415095,  39.05914193867925,  39.05905196320755,\n        39.05896198773585],\n      dtype='float64', name='y', length=1060))Attributes: (4)AREA_OR_POINT :AreaTIFFTAG_RESOLUTIONUNIT :1 (unitless)TIFFTAG_XRESOLUTION :1TIFFTAG_YRESOLUTION :1\n\n\n\n# Display CRS for raster layer\nprint(raster.rio.crs)\n\nEPSG:4326\n\n\n\n# Create separate variables to easily keep track of each band\nB4 = raster[1] # Red band\nB8 = raster[2] # NIR band\n\n\n# Visualize raster files (1 band each)\nplt.figure(figsize=(12,4))\n\nplt.subplot(1,2,1)\nplt.title('Band 4')\nB4.plot(cmap='Spectral')\n\nplt.subplot(1,2,2)\nplt.title('Band 8')\nB8.plot(cmap='Spectral')\n\nplt.show()"
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#compute-ndvi",
+    "href": "exercises/raster_and_vector.html#compute-ndvi",
+    "title": "86  Raster and vector",
+    "section": "Compute NDVI",
+    "text": "Compute NDVI\nNormalized Difference Vegetation Index (NDVI) is a measure typically used in remote sensing to quantify and assess vegetation health and density. It is calculated using near-infrared (NIR) and red light reflectance values over an area as follows:\n NDVI = \\frac{Red - NIR}{Red+NIR}\nNDVI values range from -1 to 1, where higher values indicate healthier and denser vegetation, typically falling between 0.2 and 0.8 for vegetation. Bare soil and rocks fall have values ranging from 0 to 0.2. Bodies of water like ponds, likes, rivers, and oceans typically have negative values close to -1.\n\n# Calculate NDVI\nndvi = (B8 - B4)/(B8 + B4)"
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#define-custom-colormap",
+    "href": "exercises/raster_and_vector.html#define-custom-colormap",
+    "title": "86  Raster and vector",
+    "section": "Define custom colormap",
+    "text": "Define custom colormap\n\n# Palette of colors NDVI\nhex_palette = ['#CE7E45', '#DF923D', '#F1B555', '#FCD163', '#99B718', '#74A901',\n             '#66A000', '#529400', '#3E8601', '#207401', '#056201', '#004C00', '#023B01',\n             '#012E01', '#011D01', '#011301']\n\n# Use the built-in ListedColormap function to do the conversion\nndvi_cmap = colors.ListedColormap(hex_palette)\nndvi_cmap.set_bad('#FEFEFE')\nndvi_cmap.set_under('#0000FF')\nndvi_cmap\n\nfrom_list  underbad over \n\n\n\n# Add new colormap to list of Matplotlib's colormaps\nplt.colormaps.register(cmap=ndvi_cmap, name='ndvi')\n\n\n# Visualize NDVI for the entire area of the image\nndvi.plot(cmap='ndvi', vmin=0.0, vmax=1.0)\nplt.title('NDVI')\nplt.xlabel('Longitude')\nplt.ylabel('Latitude')\nplt.show()"
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#clip-ndvi-raster-to-watersheds",
+    "href": "exercises/raster_and_vector.html#clip-ndvi-raster-to-watersheds",
+    "title": "86  Raster and vector",
+    "section": "Clip NDVI raster to watersheds",
+    "text": "Clip NDVI raster to watersheds\n\n# Define a function to clip the raster with each polygon and return a numpy array\nclip_raster = lambda polygon, raster: raster.rio.clip([polygon.geometry], \n                                                      crs=raster.rio.crs, \n                                                      all_touched=True)\n\n# Apply the function to each row in the GeoDataFrame to create a new 'clipped_raster' column\ngdf['clipped_raster'] = gdf.apply(lambda row: clip_raster(row, ndvi), axis=1)\n\n# Inspect resulting GeoDataframe\ngdf.head(3)\n\n\n\n\n\n\n\n\nCODE_1\nNAME_1\nAREA\nPERIMETER\nACRES\nHECTARES_1\nDATAID_1\nDATACODE_2\ngeometry\nclipped_raster\n\n\n\n\n0\n0R1A\nR1A\n484815.196184\n2939.146533\n119.800444\n48.481520\nGIS0321\nGIS032\nPOLYGON ((-96.60663 39.07307, -96.60671 39.073...\n[[&lt;xarray.DataArray 'band_data' ()&gt;\\narray(nan...\n\n\n1\nR20A\nR20A\n263364.016292\n2509.563160\n65.078666\n26.336402\nGIS0322\nGIS032\nPOLYGON ((-96.60627 39.07163, -96.60627 39.071...\n[[&lt;xarray.DataArray 'band_data' ()&gt;\\narray(nan...\n\n\n2\n002A\n2A\n283017.160279\n2417.724921\n69.935063\n28.301716\nGIS0323\nGIS032\nPOLYGON ((-96.60294 39.07362, -96.60292 39.073...\n[[&lt;xarray.DataArray 'band_data' ()&gt;\\narray(nan..."
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#visualize-a-specific-clipped-watershed",
+    "href": "exercises/raster_and_vector.html#visualize-a-specific-clipped-watershed",
+    "title": "86  Raster and vector",
+    "section": "Visualize a specific clipped watershed",
+    "text": "Visualize a specific clipped watershed\nWe will choose the K1B watershed to show the NDVI for the entire area.\n\n# Select watershed\nidx = gdf['NAME_1'] == 'K1B'\nrow = gdf[idx].index[0]\nprint(row)\n\n31\n\n\n\n# Create figure of selected watershed\nfig, ax = plt.subplots(figsize=(6, 6))\ngdf.loc[ [row], 'geometry'].boundary.plot(ax=ax, edgecolor='k')\ngdf.loc[row, 'clipped_raster'].plot(ax=ax, cmap='ndvi', \n                                    add_colorbar=True, vmin=0.0, vmax=1.0,\n                                   cbar_kwargs={'label':'NDVI', 'shrink':0.5})\nax.set_title('K1B')\nax.set_xlabel('Longitude')\nax.set_ylabel('Latitude')\n\nplt.show()\n\n\n\n\nNoticed that to plot the boundary of the watershed we had to do: type(gdf.loc[[row],'geometry']) instead of type(gdf.loc[row,'geometry']). This is because we need to access the GeoSeries rather than the Shapely polygon."
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#compute-mean-ndvi-for-each-watershed",
+    "href": "exercises/raster_and_vector.html#compute-mean-ndvi-for-each-watershed",
+    "title": "86  Raster and vector",
+    "section": "Compute mean NDVI for each watershed",
+    "text": "Compute mean NDVI for each watershed\n\n# Create empty list to append watershed NDVI values\nndvi_avg = []\n\n# Iterate over each watershed\nfor k,row in gdf.iterrows():\n    ndvi_avg.append(np.nanmean(row['clipped_raster'].data))\n\n# Add list of NDVI values to GeoDataframe\ngdf.insert(gdf.shape[1]-1, 'ndvi_avg', ndvi_avg)\n\n# Inspect results\ngdf.head(3)\n\n\n\n\n\n\n\n\nCODE_1\nNAME_1\nAREA\nPERIMETER\nACRES\nHECTARES_1\nDATAID_1\nDATACODE_2\ngeometry\nndvi_avg\nclipped_raster\n\n\n\n\n0\n0R1A\nR1A\n484815.196184\n2939.146533\n119.800444\n48.481520\nGIS0321\nGIS032\nPOLYGON ((-96.60663 39.07307, -96.60671 39.073...\n0.187059\n[[&lt;xarray.DataArray 'band_data' ()&gt;\\narray(nan...\n\n\n1\nR20A\nR20A\n263364.016292\n2509.563160\n65.078666\n26.336402\nGIS0322\nGIS032\nPOLYGON ((-96.60627 39.07163, -96.60627 39.071...\n0.196402\n[[&lt;xarray.DataArray 'band_data' ()&gt;\\narray(nan...\n\n\n2\n002A\n2A\n283017.160279\n2417.724921\n69.935063\n28.301716\nGIS0323\nGIS032\nPOLYGON ((-96.60294 39.07362, -96.60292 39.073...\n0.173563\n[[&lt;xarray.DataArray 'band_data' ()&gt;\\narray(nan..."
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#create-static-map",
+    "href": "exercises/raster_and_vector.html#create-static-map",
+    "title": "86  Raster and vector",
+    "section": "Create static map",
+    "text": "Create static map\n\n# Create figure with mean NDVI values for each watershed\ngdf.plot(column='ndvi_avg', edgecolor='k', cmap='ndvi', legend=True, vmin=0.0, vmax=1.0)\nplt.title('Konza Prairie Mean NDVI April 2024')\nplt.xlabel('Longitude')\nplt.ylabel('Latitude')\n#plt.savefig('Konza Mean NDVI APril 2024.jpg', dpi=300)\nplt.show()"
+  },
+  {
+    "objectID": "exercises/raster_and_vector.html#create-interactive-map",
+    "href": "exercises/raster_and_vector.html#create-interactive-map",
+    "title": "86  Raster and vector",
+    "section": "Create interactive map",
+    "text": "Create interactive map\n\n# Create interactive map\ngdf.iloc[:,:-1].explore(column='ndvi_avg', k=10, cmap='ndvi',\n                        style_kwds={'stroke':True,'color':'black','width':1},\n                        legend_kwds={'caption':'NDVI'})\n\n\nMake this Notebook Trusted to load map: File -&gt; Trust Notebook"
   }
 ]
diff --git a/exercises/count_seeds.ipynb b/exercises/count_seeds.ipynb
@@ -211,8 +211,7 @@
    "source": [
     "# Closing to fill in holes or connect small, nearby patches\n",
     "kernel_size = 5\n",
-    "image_binary = binary_closing(image_binary, \n",
-    "                              disk(kernel_size))\n",
+    "image_binary = binary_closing(image_binary, disk(kernel_size))\n",
     "\n",
     "# Let's inspect the structuring element\n",
     "print(disk(kernel_size))\n"
@@ -333,7 +332,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 28,
    "metadata": {},
    "outputs": [
     {
@@ -353,17 +352,12 @@
     "\n",
     "for seed in range(36):\n",
     "    plt.subplot(6,6,seed+1)\n",
-    "    plt.plot(contours[seed][:, 1], contours[seed][:, 0], '-r', linewidth=2)\n",
+    "    plt.plot(contours[seed][:, 1], \n",
+    "             contours[seed][:, 0], \n",
+    "             '-r', linewidth=2)\n",
     "    plt.tight_layout()\n",
     "plt.show()\n"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

diff --git a/exercises/green_ampt.ipynb b/exercises/green_ampt.ipynb
@@ -33,27 +33,34 @@
     "\n",
     "#### 2) Cumulative infiltration at time of ponding\n",
     "\n",
-    "$$ F_p = \\frac{K_{sat}P}{w-K_{sat}}$$\n",
+    "$$ F_p = \\frac{K_{sat}P}{w-K_{sat}} ~~~ , w>K_{sat}$$\n",
     "\n",
     "$F_p$ is the cumulative infiltration when ponding occurs<br>\n",
     "$w$ is the rainfall rate\n",
     "\n",
     "#### 3) Cumulative infiltration under ponded conditions\n",
     "\n",
+    "$$ 0 = t - t_p - \\frac{F-F_p}{K_{sat}} + \\frac{P}{K_{sat}} ln \\Bigg(\\frac{F_p+P}{F+P}  \\Bigg)$$\n",
+    "\n",
+    "where $P$ is the moisture difference suction parameter. We solve this function implicitly for F\n",
+    "\n",
+    "\n",
     "\n",
     "In the code below, we are also using the following variables:\n",
+    "\n",
     "$dt^{'}$ is the partial time interval required for ponding to occur\n",
     "\n",
     "$t_p$ is time during which ponding occurs\n",
     "\n",
-    "$f_t$ is infiltration during time step\n",
+    "$f_t$ is infiltration during time step $t$\n",
     "\n",
-    "$r_t$ is runoff during time step.\n",
+    "$r_t$ is runoff during time step $t$\n",
     "\n",
-    "The model works on the premise that, if rainfall rate is lower than the infiltration capacity ($f_c$), then all the rainfall infiltrates the soil. However, if the rainfall rate is greater than the infiltration capacity, then ponding occurs. Much of the coding needed to solve the Green-Ampt model deals with what happens within the considered time interval ($dt$). In other words, the premise that we stated above, is for the beginning of the time interval. Even in the case where there is no ponding at the beginning of the time interval, it's still possible for ponding to occur within the time interval. To resolve this, tentative values of infiltration capacity ($f_c^{'}$) and cumulative infiltration ($F^{'}$) are computed within the interval assuming that all the water infiltrates. If the rainfall rate is greater than ($f_c^{'}$), then ponding occurs at some point during the time step. Note that if ponding does not occur during the entire interval, computations are simpler and we move onto the next step. Additional calculations using a solver are required when we find that ponding occurs at some point during the time interval.\n",
+    "The model works on the premise that, if rainfall rate is lower than the infiltration capacity ($f_c$), then all the rainfall infiltrates the soil. However, if the rainfall rate is greater than the infiltration capacity of the soil, then ponding occurs. \n",
     "\n",
+    "Much of the coding needed to solve the Green-Ampt model deals with what happens within the time interval ($dt$). In other words, even in the case where there is no ponding at the beginning of the time interval, it's still possible for ponding to occur within the time interval. To resolve this, tentative values of infiltration capacity ($f_c^{'}$) and cumulative infiltration ($F^{'}$) are computed within the interval assuming that all the water infiltrates. If the rainfall rate is greater than ($f_c^{'}$), then ponding occurs at some point during the time step. Note that if ponding does not occur during the entire interval, computations are simpler and we move onto the next step. Additional calculations using a solver are required when we find that ponding occurs at some point during the time interval.\n",
     "\n",
-    "Something important to remember is that rainfall rate can increase or decrease over time, but infiltration capacity decreases as more water infiltrates the soil. Thus, "
+    "Something important to remember is that rainfall rate can increase or decrease over time, but the infiltration capacity of the soil decreases as more water infiltrates the soil. "
    ]
   },
   {