Skip to content

Commit

Permalink
Add the hydrofabric subsetting jupyter notebook (#54)
Browse files Browse the repository at this point in the history 
* Create readme.md for the hydrofab-subsetting folder

* Add the hydrofabric subset jupyter notebook

The jupyter notebook uses the subset.py script. This is similar to the one available at ../subsetting/subset.py with a few changes. The current subset.py retrieves the hydrofabric v1.2 data.
  • Loading branch information
@igarousi
igarousi authored Sep 11, 2023
1 parent 195b44b commit c3bf273
Show file tree
Hide file tree
Showing 3 changed files with 784 additions and 0 deletions.
1 change: 1 addition & 0 deletions hydrofab-subsetting/README.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
This contains subsetting jupyter notebooks developed at CUAHSI.
370 changes: 370 additions & 0 deletions hydrofab-subsetting/hydrofabric-subset-sidecar.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,370 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "cc3872d5-fe42-49a0-aee9-823d58af0e11",
"metadata": {},
"source": [
"# Subset NGEN HydroFabric on S3\n",
"\n",
"**Authors:** \n",
" - Tony Castronova <[email protected]> \n",
" - Irene Garousi-Nejad <[email protected]> \n",
" \n",
"**Last Updated:** 04.21.2023 \n",
"\n",
"**Description**: \n",
"\n",
"The purpose of this Jupyter Notebook is to demonstrate how to extract the National Hydrologic Geospatial (Hydrofabric) at a specific area of interest. Once the Hydrofabric has been extracted, it can be used to execute the [NOAA Next Generation (NextGen) Water Resource Modeling Framework](https://github.com/NOAA-OWP/ngen). The Hydrofabric data are publicly available in the AWS catalog. In this notebook, we use Version 1.2 of the NGen Hydrofabric which is the most recent version available on the Amazon S3 Bucket as of the time of developing this notebook (https://nextgen-hydrofabric.s3.amazonaws.com/index.html#v1.2/).\n",
"\n",
"**Software Requirements**: \n",
"\n",
"The software and operating system versions used to develop this notebook are listed below. To avoid encountering issues related to version conflicts among Python packages, we recommend creating a new environment variable and installing the required packages specifically for this notebook.\n",
"\n",
"Tested on: MacOS Ventura 13.2.1 \n",
"\n",
"> boto3: 1.26.76 \n",
" dask-core: 2023.4.0 \n",
" fiona: 1.9.3 \n",
" fsspec: 2023.4.0 \n",
" geopandas: 0.12.2 \n",
" ipyleaflet: 0.17.2 \n",
" ipywidgets: 7.7.5 \n",
" matplotlib: 3.7.1 \n",
" netcdf4: 1.6.3 \n",
" numpy: 1.24.2 \n",
" pandas: 2.0.0 \n",
" requests: 2.28.2 \n",
" s3fs: 2023.4.0 \n",
" scipy: 1.10.1 \n",
" xarray: 2023.4.1\n",
" \n",
"**Supplementary Code**\n",
"\n",
"This notebook relies on the following external scripts: \n",
"- `subset.py` - A script originally written by Nels Frazier to subset the NGen Hydrofabric"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c354f815-4e05-4a57-86d1-213a1536d801",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import time\n",
"import pandas\n",
"import subset\n",
"import pyproj\n",
"import datetime\n",
"import geopandas\n",
"import ipyleaflet\n",
"import geopandas as gpd\n",
"from pathlib import Path\n",
"from sidecar import Sidecar\n",
"from requests import Request\n",
"from ipywidgets import Layout"
]
},
{
"cell_type": "markdown",
"id": "c9c672fc-9635-460b-bef0-64338759f036",
"metadata": {},
"source": [
"\n",
"Create a map and load the HydroFabric VPU geometries. These geometries have been prepared ahead of time and are stored in a HydroShare [resource](https://www.hydroshare.org/resource/35e8c6023c154b6298fcda280beda849/). HydroShare provides [WMS](https://docs.geoserver.org/latest/en/user/services/wms/index.html) and [WFS](https://docs.geoserver.org/latest/en/user/services/wfs/index.html) capabilities that enable us to easily display them on an interactive map.\n",
"\n",
"Create the map using the following code:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a8c882f2-484a-4d14-9bdc-f0c1f2a3152f",
"metadata": {},
"outputs": [],
"source": [
"defaultLayout=Layout(width='960px', height='940px')\n",
"\n",
"map_center = (41.74614949822607, -111.76617850993877) # logan, UT\n",
"m = ipyleaflet.Map(\n",
" basemap=ipyleaflet.basemap_to_tiles(ipyleaflet.basemaps.OpenStreetMap.Mapnik, layout=defaultLayout),\n",
" center=map_center,\n",
" zoom=9,\n",
" scroll_wheel_zoom=True,\n",
" tap=False\n",
" )\n",
"\n",
"# # add rivers\n",
"# m.add_layer(\n",
"# ipyleaflet.WMSLayer(\n",
"# url='https://hydro.nationalmap.gov/arcgis/services/nhd/MapServer/WMSServer',\n",
"# # url='http://arcgis.cuahsi.org/arcgis/services/NHD/usgs_gages/MapServer/WmsServer',\n",
"# layers='4',\n",
"# transparent=True,\n",
"# format='image/png',\n",
"# min_zoom=8,\n",
"# max_zoom=18,\n",
"# )\n",
"# )\n",
"\n",
"# add USGS Gages\n",
"m.add_layer(\n",
" ipyleaflet.WMSLayer(\n",
" url='http://arcgis.cuahsi.org/arcgis/services/NHD/usgs_gages/MapServer/WmsServer',\n",
" layers='0',\n",
" transparent=True,\n",
" format='image/png',\n",
" min_zoom=8,\n",
" max_zoom=18,\n",
" )\n",
")\n",
"\n",
"# add the CONUS VPU boundaries\n",
"m.add_layer(\n",
" ipyleaflet.WMSLayer(\n",
" url='https://geoserver.hydroshare.org/geoserver/HS-e8ddee6a8a90484fa7a976458e79c0c3/wms?',\n",
" layers='vpu_boundaries',\n",
" format='image/png',\n",
" transparent=True,\n",
" opacity=0.5,\n",
" min_zoom=4,\n",
" max_zoom=8\n",
" )\n",
")\n",
"\n",
"# add the watershed VPU boundaries for each region.\n",
"for vpu in ['16']: \n",
" #['01', '02','03N','03S','03W', '04','05','06','07','08','09','10L','10U','11','12','13','14','15','16','17','18']:\n",
" m.add_layer(\n",
" ipyleaflet.WMSLayer(\n",
" url='https://geoserver.hydroshare.org/geoserver/HS-e8ddee6a8a90484fa7a976458e79c0c3/wms?',\n",
" layers=f'{vpu}_boundaries',\n",
" format='image/png',\n",
" transparent=True,\n",
" opacity=0.5,\n",
" min_zoom=8,\n",
" max_zoom=18\n",
" )\n",
" )\n"
]
},
{
"cell_type": "markdown",
"id": "70d69c43-6c97-4e2d-8a9e-c24f1ee74174",
"metadata": {},
"source": [
"Add an event handler that will enable us to highlight and store geometries that have been clicked on the map. When a geometry is clicked on the map, we'll call the `WFS` endpoint to retrieve the boundary of the shape. We can draw this boundary on the map and store information about the selected area for later use."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "807970de-8af1-43d0-8511-c616e2b7eab5",
"metadata": {},
"outputs": [],
"source": [
"selected_df = None\n",
"\n",
"def handle_map_interaction(**kwargs):\n",
" global selected_df\n",
" \n",
" if kwargs.get('type') == 'click':\n",
" coords = kwargs['coordinates'] \n",
" url = f'https://geoserver.hydroshare.org/geoserver/HS-e8ddee6a8a90484fa7a976458e79c0c3/wfs?' \\\n",
" 'service=wfs&version=2.0.0&' \\\n",
" f'request=getFeature&' \\\n",
" 'srsName=EPSG:4269&' \\\n",
" f'bbox={coords[1]},{coords[0]},{coords[1]},{coords[0]},EPSG:4269&' \\\n",
" f'typeName=vpu_boundaries&' \\\n",
" 'outputFormat=json&' \\\n",
" 'PropertyName=VPU'\n",
" q = Request('GET', url).prepare().url\n",
" df = gpd.read_file(q, format='json')\n",
" \n",
" # exit if a VPU is not found, i.e. a user doesn't click on the layer\n",
" if len(df) == 0: return\n",
" \n",
" VPU = df.VPU.values[0]\n",
" url = f'https://geoserver.hydroshare.org/geoserver/HS-e8ddee6a8a90484fa7a976458e79c0c3/wfs?' \\\n",
" 'service=wfs&version=2.0.0&' \\\n",
" f'request=getFeature&' \\\n",
" 'srsName=EPSG:4269&' \\\n",
" f'bbox={coords[1]},{coords[0]},{coords[1]},{coords[0]},EPSG:4269&' \\\n",
" f'typeName={VPU}_boundaries&' \\\n",
" 'outputFormat=json&'\n",
" q = Request('GET', url).prepare().url\n",
" df = gpd.read_file(q, format='json')\n",
"\n",
" # exit if a VPU is not found, i.e. a user doesn't click on the layer\n",
" if len(df) == 0: return\n",
" \n",
" # save vpu region, convert crs, and save selection for later\n",
" df['VPU'] = VPU \n",
" df = df.to_crs('EPSG:4269')\n",
" selected_df = df\n",
" \n",
" if type(m.layers[-1]) == ipyleaflet.WKTLayer:\n",
" m.remove_layer(m.layers[-1])\n",
" \n",
" # display the watershed boundary on the map\n",
" m.add_layer(ipyleaflet.WKTLayer(wkt_string=df.iloc[0].geometry.wkt))\n",
" \n",
"m.on_interaction(handle_map_interaction)"
]
},
{
"cell_type": "markdown",
"id": "091faa22-689e-4e2f-8e3c-02244e530d2e",
"metadata": {},
"source": [
"Display the map"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c1b6fa44-1657-42c5-9fa4-b550eede388f",
"metadata": {},
"outputs": [],
"source": [
"sc = Sidecar(title='NGEN HydroFabric Map')\n",
"with sc:\n",
" display(m)"
]
},
{
"cell_type": "markdown",
"id": "d1b8729c-fd1e-4206-83b0-1fc0d1fbe0c8",
"metadata": {},
"source": [
"Pass the `ids` of the geometries selected in the map to the hydrofabric subsetting script."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2ac5a477-93df-48d3-b3a6-f9b09d606665",
"metadata": {},
"outputs": [],
"source": [
"selected_df.id"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "28af8c32-011d-4981-a614-1d5f62de424c",
"metadata": {},
"outputs": [],
"source": [
"\n",
"ids = list(selected_df.id.values)\n",
"vpus = list(selected_df.VPU.values)\n",
"\n",
"output_files = []\n",
"for i in range(0, len(ids)): \n",
" print(50*'-'+f'\\nProcessing VPU {vpus[i]}, {ids[i]} ')\n",
" st = time.time()\n",
" # build the hydrofabric_url\n",
" hydrofabric_url = f's3://nextgen-hydrofabric/pre-release/nextgen_{vpus[i]}.gpkg'\n",
" subset.subset_upstream_prerelease(hydrofabric_url, ids[i])\n",
" \n",
" # move these files into a subdir to keep things orderly\n",
" counter = 1\n",
" outpath = ids[i]\n",
" while os.path.exists(outpath):\n",
" outpath = ids[i] + \" (\" + str(counter) + \")\"\n",
" counter += 1\n",
" os.mkdir(outpath)\n",
" for f in [f'{ids[i]}_upstream_subset.gpkg',\n",
" 'catchments.geojson',\n",
" 'crosswalk.json',\n",
" 'flowpath_edge_list.json',\n",
" 'flowpaths.geojson',\n",
" 'nexus.geojson']:\n",
" os.rename(f, f'{outpath}/{f}')\n",
" \n",
" # output_files.append(f'{ids[i]}_upstream_subset.gpkg')\n",
" print(f'Output files located at: {outpath}')\n",
" print(f'Completed in {time.time() - st} seconds\\n'+50*'-') \n",
"\n",
"outdir = Path(outpath)"
]
},
{
"cell_type": "markdown",
"id": "e6410e7b-f225-406d-adb9-c7cefdc07171",
"metadata": {},
"source": [
"### Add the Subset HydroFabric to the Map\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c042621b-5525-41cf-97ee-2c5bfe886921",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# read the shapes\n",
"catchments = geopandas.read_file(f'{outdir/outdir.name}_upstream_subset.gpkg',\n",
" layer='divides')\n",
"rivers = geopandas.read_file(f'{outdir/outdir.name}_upstream_subset.gpkg',\n",
" layer='flowpaths')\n",
"\n",
"# define the target projection as EPSG:4269 - Web Mercator\n",
"# this is the default crs for leaflet\n",
"target_crs = pyproj.Proj('4269')\n",
"\n",
"# transform the shapefile into EPSG:4269\n",
"catchments = catchments.to_crs(target_crs.crs)\n",
"rivers = rivers.to_crs(target_crs.crs)\n",
"\n",
"# add the catchments to the map\n",
"for idx, shape in catchments.iterrows():\n",
" wkt = ipyleaflet.WKTLayer(wkt_string=shape.geometry.wkt)\n",
" wkt.style = {'color': 'green'}\n",
" m.add_layer(wkt)\n",
" \n",
"# add the rivers to the map\n",
"for idx, shape in rivers.iterrows():\n",
" wkt = ipyleaflet.WKTLayer(wkt_string=shape.geometry.wkt)\n",
" wkt.style = {'color': 'blue'}\n",
" m.add_layer(wkt)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a888d53d-47cd-4ff9-b955-87e163bce26c",
"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.10.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Loading

0 comments on commit c3bf273

Please sign in to comment.