From 8ddadb9cb02ee0d00a3ea154e97d6013dd71246c Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Fri, 27 Sep 2024 14:51:07 -0400 Subject: [PATCH 1/6] Add support for downloading NWI wetlands --- leafmap/common.py | 126 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) diff --git a/leafmap/common.py b/leafmap/common.py index 4428b87f4e..28e25ac75e 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -14362,3 +14362,129 @@ def get_nhd( gdf = None return gdf + + +def get_nwi( + geometry: Dict[str, Any], + inSR: str = "4326", + outSR: str = "3857", + spatialRel: str = "esriSpatialRelIntersects", + return_geometry: bool = True, + outFields: str = "*", + output: Optional[str] = None, + **kwargs: Any, +) -> Union["gpd.GeoDataFrame", "pd.DataFrame", Dict[str, str]]: + """ + Query the NWI (National Wetlands Inventory) API using various geometry types. + https://fwspublicservices.wim.usgs.gov/wetlandsmapservice/rest/services/Wetlands/FeatureServer + + Args: + geometry (dict): The geometry data (e.g., point, polygon, polyline, multipoint, etc.). + inSR (str): The input spatial reference (default is EPSG:4326). + outSR (str): The output spatial reference (default is EPSG:3857). + spatialRel (str): The spatial relationship (default is "esriSpatialRelIntersects"). + return_geometry (bool): Whether to return the geometry (default is True). + outFields (str): The fields to be returned (default is "*"). + output (str): The output file path to save the GeoDataFrame (default is None). + **kwargs: Additional keyword arguments to pass to the API. + + Returns: + gpd.GeoDataFrame: The queried NWI data as a GeoDataFrame. + """ + + import geopandas as gpd + import pandas as pd + from shapely.geometry import Polygon + + def detect_geometry_type(geometry): + """ + Automatically detect the geometry type based on the keys in the geometry dictionary. + + Args: + geometry (dict): The geometry data (e.g., point, polygon, polyline, multipoint, etc.). + + Returns: + str: The detected geometry type (e.g., "esriGeometryPoint", "esriGeometryPolygon"). + """ + if "x" in geometry and "y" in geometry: + return "esriGeometryPoint" + elif ( + "xmin" in geometry + and "ymin" in geometry + and "xmax" in geometry + and "ymax" in geometry + ): + return "esriGeometryEnvelope" + elif "rings" in geometry: + return "esriGeometryPolygon" + elif "paths" in geometry: + return "esriGeometryPolyline" + elif "points" in geometry: + return "esriGeometryMultipoint" + else: + raise ValueError("Unsupported geometry type or invalid geometry structure") + + # Automatically detect the geometry type based on the structure of the geometry + geometry_type = detect_geometry_type(geometry) + + # API URL for querying wetlands + url = "https://fwspublicservices.wim.usgs.gov/wetlandsmapservice/rest/services/Wetlands/FeatureServer/0/query" + + # Convert geometry to a JSON string (required by the API) + geometry_json = json.dumps(geometry) + + # Construct the query parameters + params = { + "geometry": geometry_json, # The geometry as a JSON string + "geometryType": geometry_type, # Geometry type (automatically detected) + "inSR": inSR, # Spatial reference system (default is WGS84) + "spatialRel": spatialRel, # Spatial relationship (default is intersects) + "outFields": outFields, # Which fields to return (default is all fields) + "returnGeometry": str( + return_geometry + ).lower(), # Whether to return the geometry + "f": "json", # Response format + } + + for key, value in kwargs.items(): + params[key] = value + + # Make the GET request + response = requests.get(url, params=params) + + # Check if the request was successful + if response.status_code == 200: + data = response.json() # Return the data as a Python dictionary + else: + return {"error": f"Request failed with status code {response.status_code}"} + + # Extract the features + features = data["features"] + + # Prepare the attribute data and geometries + attributes = [feature["attributes"] for feature in features] + + # Create a DataFrame for attributes + df = pd.DataFrame(attributes) + df.rename( + columns={"Shape__Length": "Shape_Length", "Shape__Area": "Shape_Area"}, + inplace=True, + ) + + if return_geometry: + geometries = [Polygon(feature["geometry"]["rings"][0]) for feature in features] + # Create a GeoDataFrame by combining the attributes and geometries + gdf = gpd.GeoDataFrame( + df, + geometry=geometries, + crs=f"EPSG:{data['spatialReference']['latestWkid']}", + ) + if outSR != "3857": + gdf = gdf.to_crs(outSR) + + if output is not None: + gdf.to_file(output) + + return gdf + else: + return df From d95072e35d1b91d385c510600ac4e8d94148ceff Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Fri, 27 Sep 2024 14:54:04 -0400 Subject: [PATCH 2/6] Disable uv cache --- .github/workflows/docs-build.yml | 2 +- .github/workflows/docs.yml | 2 +- .github/workflows/installation.yml | 2 +- .github/workflows/macos.yml | 2 +- .github/workflows/py313.yml | 2 +- .github/workflows/ubuntu.yml | 2 +- .github/workflows/windows.yml | 2 +- 7 files changed, 7 insertions(+), 7 deletions(-) diff --git a/.github/workflows/docs-build.yml b/.github/workflows/docs-build.yml index 45230464da..5a1ebcbfc8 100644 --- a/.github/workflows/docs-build.yml +++ b/.github/workflows/docs-build.yml @@ -31,7 +31,7 @@ jobs: uses: astral-sh/setup-uv@v3 with: version: "0.4.16" - enable-cache: true + # enable-cache: true - name: Set up Python ${{ matrix.python-version }} run: uv python install ${{ matrix.python-version }} diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 0d8d2f1d69..1cfd686e5d 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -31,7 +31,7 @@ jobs: uses: astral-sh/setup-uv@v3 with: version: "0.4.16" - enable-cache: true + # enable-cache: true - name: Set up Python ${{ matrix.python-version }} run: uv python install ${{ matrix.python-version }} diff --git a/.github/workflows/installation.yml b/.github/workflows/installation.yml index bcd247a99c..6c9bd1c55d 100644 --- a/.github/workflows/installation.yml +++ b/.github/workflows/installation.yml @@ -21,7 +21,7 @@ jobs: uses: astral-sh/setup-uv@v3 with: version: "0.4.16" - enable-cache: true + # enable-cache: true - name: Set up Python ${{ matrix.python-version }} run: uv python install ${{ matrix.python-version }} diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index 56f36e6674..77226270da 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -25,7 +25,7 @@ jobs: uses: astral-sh/setup-uv@v3 with: version: "0.4.16" - enable-cache: true + # enable-cache: true - name: Set up Python ${{ matrix.config.py }} run: uv python install ${{ matrix.config.py }} diff --git a/.github/workflows/py313.yml b/.github/workflows/py313.yml index 1f3c880178..9cec996dcb 100644 --- a/.github/workflows/py313.yml +++ b/.github/workflows/py313.yml @@ -33,7 +33,7 @@ jobs: uses: astral-sh/setup-uv@v3 with: version: "0.4.16" - enable-cache: true + # enable-cache: true - name: Set up Python ${{ matrix.python-version }} run: uv python install ${{ matrix.python-version }} diff --git a/.github/workflows/ubuntu.yml b/.github/workflows/ubuntu.yml index 7ce4c0d9e8..9c5560432e 100644 --- a/.github/workflows/ubuntu.yml +++ b/.github/workflows/ubuntu.yml @@ -33,7 +33,7 @@ jobs: uses: astral-sh/setup-uv@v3 with: version: "0.4.16" - enable-cache: true + # enable-cache: true - name: Set up Python ${{ matrix.python-version }} run: uv python install ${{ matrix.python-version }} diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 8d8b816fd8..d7826a5863 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -30,7 +30,7 @@ jobs: uses: astral-sh/setup-uv@v3 with: version: "0.4.16" - enable-cache: true + # enable-cache: true - name: Set up Python ${{ matrix.python-version }} run: uv python install ${{ matrix.python-version }} From fb6b1dd8d8bf60575bdcbe18d2b413fcc0b1ba2d Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Fri, 27 Sep 2024 22:52:30 -0400 Subject: [PATCH 3/6] Add get_wbd function --- leafmap/common.py | 229 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 215 insertions(+), 14 deletions(-) diff --git a/leafmap/common.py b/leafmap/common.py index 28e25ac75e..cf4dd2c41c 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -14303,7 +14303,7 @@ def construct_bbox( return geometry -def get_nhd( +def get_nhd_wbd( geometry: Union[ "gpd.GeoDataFrame", str, List[float], Tuple[float, float, float, float] ], @@ -14364,6 +14364,53 @@ def get_nhd( return gdf +def _convert_geometry_to_esri_format(geometry): + """ + Convert shapely geometry to ESRI's format for point, polygon, polyline, or multipoint. + + Args: + geometry (shapely.geometry.base.BaseGeometry): The shapely geometry to convert. + + Returns: + dict: The geometry in ESRI format. + """ + from shapely.geometry import Point, Polygon, LineString, MultiPoint + + if isinstance(geometry, Point): + # Convert point to ESRI format + return {"x": geometry.x, "y": geometry.y} + elif isinstance(geometry, Polygon): + # Convert polygon to ESRI format (rings) + return {"rings": [list(geometry.exterior.coords)]} + elif isinstance(geometry, LineString): + # Convert polyline to ESRI format (paths) + return {"paths": [list(geometry.coords)]} + elif isinstance(geometry, MultiPoint): + # Convert multipoint to ESRI format (points) + return {"points": [list(point.coords[0]) for point in geometry.geoms]} + else: + raise ValueError(f"Unsupported geometry type: {type(geometry)}") + + +def _convert_geodataframe_to_esri_format(gdf: "gpd.GeoDataFrame"): + """ + Convert all geometries in a GeoDataFrame to ESRI's format. + + Args: + gdf (geopandas.GeoDataFrame): A GeoDataFrame containing geometries. + + Returns: + list of dict: A list of geometries in ESRI format. + """ + esri_geometries = [] + + for geom in gdf.geometry: + esri_format = _convert_geometry_to_esri_format(geom) + esri_geometries.append(esri_format) + + return esri_geometries + + def get_nwi( geometry: Dict[str, Any], inSR: str = "4326", @@ -14398,13 +14445,7 @@ def get_nwi( def detect_geometry_type(geometry): """ - Automatically detect the geometry type based on the keys in the geometry dictionary. - - Args: - geometry (dict): The geometry data (e.g., point, polygon, polyline, multipoint, etc.). - - Returns: - str: The detected geometry type (e.g., "esriGeometryPoint", "esriGeometryPolygon"). + Automatically detect the geometry type based on the structure of the geometry dictionary. """ if "x" in geometry and "y" in geometry: return "esriGeometryPoint" @@ -14422,17 +14463,30 @@ def detect_geometry_type(geometry): elif "points" in geometry: return "esriGeometryMultipoint" else: - raise ValueError("Unsupported geometry type or invalid geometry structure") + raise ValueError("Unsupported geometry type or invalid geometry structure.") - # Automatically detect the geometry type based on the structure of the geometry - geometry_type = detect_geometry_type(geometry) + # Convert GeoDataFrame to a dictionary if needed + if isinstance(geometry, gpd.GeoDataFrame): + geometry_dict = _convert_geodataframe_to_esri_format(geometry)[0] + geometry_type = detect_geometry_type(geometry_dict) + elif isinstance(geometry, dict): + geometry_type = detect_geometry_type(geometry) + geometry_dict = geometry + elif isinstance(geometry, str): + geometry_dict = geometry + else: + raise ValueError( + "Invalid geometry input. Must be a GeoDataFrame or a dictionary." + ) + # Convert geometry to a JSON string (required by the API) + if isinstance(geometry_dict, dict): + geometry_json = json.dumps(geometry_dict) + else: + geometry_json = geometry_dict # API URL for querying wetlands url = "https://fwspublicservices.wim.usgs.gov/wetlandsmapservice/rest/services/Wetlands/FeatureServer/0/query" - # Convert geometry to a JSON string (required by the API) - geometry_json = json.dumps(geometry) - # Construct the query parameters params = { "geometry": geometry_json, # The geometry as a JSON string @@ -14488,3 +14542,150 @@ def detect_geometry_type(geometry): return gdf else: return df + + +def get_wbd( + geometry: Union["gpd.GeoDataFrame", Dict[str, Any]], + inSR: str = "4326", + outSR: str = "3857", + digit: int = 8, + spatialRel: str = "esriSpatialRelIntersects", + return_geometry: bool = True, + outFields: str = "*", + output: Optional[str] = None, + **kwargs: Any, +) -> Union["gpd.GeoDataFrame", "pd.DataFrame", Dict[str, str]]: + """ + Query the WBD (Watershed Boundary Dataset) API using various geometry types or a GeoDataFrame. + https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer + + Args: + geometry (Union[gpd.GeoDataFrame, Dict]): The geometry data (GeoDataFrame or geometry dict). + inSR (str): The input spatial reference (default is EPSG:4326). + outSR (str): The output spatial reference (default is EPSG:3857). + digit (int): The digit code for the WBD layer (default is 8). + spatialRel (str): The spatial relationship (default is "esriSpatialRelIntersects"). + return_geometry (bool): Whether to return the geometry (default is True). + outFields (str): The fields to be returned (default is "*"). + output (Optional[str]): The output file path to save the GeoDataFrame (default is None). + **kwargs: Additional keyword arguments to pass to the API. + + Returns: + gpd.GeoDataFrame or pd.DataFrame: The queried WBD data as a GeoDataFrame or DataFrame. + """ + + import geopandas as gpd + import pandas as pd + from shapely.geometry import Polygon + + def detect_geometry_type(geometry): + """ + Automatically detect the geometry type based on the structure of the geometry dictionary. + """ + if "x" in geometry and "y" in geometry: + return "esriGeometryPoint" + elif ( + "xmin" in geometry + and "ymin" in geometry + and "xmax" in geometry + and "ymax" in geometry + ): + return "esriGeometryEnvelope" + elif "rings" in geometry: + return "esriGeometryPolygon" + elif "paths" in geometry: + return "esriGeometryPolyline" + elif "points" in geometry: + return "esriGeometryMultipoint" + else: + raise ValueError("Unsupported geometry type or invalid geometry structure.") + + # Convert GeoDataFrame to a dictionary if needed + if isinstance(geometry, gpd.GeoDataFrame): + geometry_dict = _convert_geodataframe_to_esri_format(geometry)[0] + geometry_type = detect_geometry_type(geometry_dict) + elif isinstance(geometry, dict): + geometry_type = detect_geometry_type(geometry) + geometry_dict = geometry + elif isinstance(geometry, str): + geometry_dict = geometry + else: + raise ValueError( + "Invalid geometry input. Must be a GeoDataFrame or a dictionary." + ) + + # Convert geometry to a JSON string (required by the API) + if isinstance(geometry_dict, dict): + geometry_json = json.dumps(geometry_dict) + else: + geometry_json = geometry_dict + + allowed_digit_values = [2, 4, 6, 8, 10, 12, 14, 16] + if digit not in allowed_digit_values: + raise ValueError( + f"Invalid digit value. Allowed values are {allowed_digit_values}" + ) + + layer = allowed_digit_values.index(digit) + 1 + + # API URL for querying the WBD + url = f"https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/{layer}/query" + + # Construct the query parameters + params = { + "geometry": geometry_json, + "geometryType": geometry_type, + "inSR": inSR, + "spatialRel": spatialRel, + "outFields": outFields, + "returnGeometry": str(return_geometry).lower(), + "f": "json", + } + + # Add additional keyword arguments + for key, value in kwargs.items(): + params[key] = value + + # Make the GET request + response = requests.get(url, params=params) + + if response.status_code != 200: + return {"error": f"Request failed with status code {response.status_code}"} + + data = response.json() + + # Extract features from the API response + features = data.get("features", []) + + # Prepare attribute data and geometries + attributes = [feature["attributes"] for feature in features] + df = pd.DataFrame(attributes) + df.rename( + columns={"Shape__Length": "Shape_Length", "Shape__Area": "Shape_Area"}, + inplace=True, + ) + + # Handle geometries + if return_geometry: + geometries = [ + ( + Polygon(feature["geometry"]["rings"][0]) + if "rings" in feature["geometry"] + else None + ) + for feature in features + ] + gdf = gpd.GeoDataFrame( + df, + geometry=geometries, + crs=f"EPSG:{data['spatialReference']['latestWkid']}", + ) + if outSR != "3857": + gdf = gdf.to_crs(outSR) + + if output is not None: + gdf.to_file(output) + + return gdf + else: + return df From e8504fd6efcd04e3cd0dcf69183c8087088d9fa5 Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Fri, 27 Sep 2024 23:43:27 -0400 Subject: [PATCH 4/6] Add get_nwi_by_huc8 function --- leafmap/common.py | 69 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/leafmap/common.py b/leafmap/common.py index cf4dd2c41c..0063084f5f 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -14689,3 +14689,72 @@ def detect_geometry_type(geometry): return gdf else: return df + + +def get_nwi_by_huc8( + huc8: Optional[str] = None, + geometry: Optional[Union["gpd.GeoDataFrame", str]] = None, + out_dir: Optional[str] = None, + quiet: bool = True, + layer: str = "Wetlands", + **kwargs, +) -> "gpd.GeoDataFrame": + """ + Fetches National Wetlands Inventory (NWI) data by HUC8 code. + + Args: + huc8 (Optional[str]): The HUC8 code to query the NWI data. It must be a + string of length 8. + geometry (Optional[Union[gpd.GeoDataFrame, str]]): The geometry to derive + the HUC8 code. It can be a GeoDataFrame or a file path. + out_dir (Optional[str]): The directory to save the downloaded data. + Defaults to a temporary directory. + quiet (bool): Whether to suppress download progress messages. Defaults to True. + layer (str): The layer to fetch from the NWI data. It can be one of the following: + Wetlands, Watershed, Riparian_Project_Metadata, Wetlands_Historic_Map_Info. + Defaults to "Wetlands". + **kwargs: Additional keyword arguments to pass to the download_file function. + + Returns: + gpd.GeoDataFrame: The fetched NWI data as a GeoDataFrame. + + Raises: + ValueError: If the HUC8 code is invalid or the layer is not allowed. + """ + import tempfile + import geopandas as gpd + + if geometry is not None: + wbd = get_wbd(geometry, return_geometry=False) + huc8 = wbd["huc8"].values[0] + + if isinstance(huc8, str) and len(huc8) == 8: + pass + else: + raise ValueError("Invalid HUC8 code. It must be a string of length 8.") + + if out_dir is None: + out_dir = tempfile.gettempdir() + + allowed_layers = [ + "Wetlands", + "Watershed", + "Riparian_Project_Metadata", + "Wetlands_Historic_Map_Info", + "Wetlands_Project_Metadata", + ] + if layer not in allowed_layers: + raise ValueError(f"Invalid layer. Allowed values are {allowed_layers}") + + url = f"https://documentst.ecosphere.fws.gov/wetlands/downloads/watershed/HU8_{huc8}_Watershed.zip" + + filename = os.path.join(out_dir, f"HU8_{huc8}_Watershed.zip") + + download_file(url, filename, quiet=quiet, **kwargs) + + data_dir = os.path.join(out_dir, f"HU8_{huc8}_Watershed") + + filepath = os.path.join(data_dir, f"HU8_{huc8}_{layer}.shp") + + gdf = gpd.read_file(filepath) + return gdf From 567b6eefad433ff6eaf0ecb85c51030540242d9b Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Sat, 28 Sep 2024 00:49:13 -0400 Subject: [PATCH 5/6] Add nhd notebook example --- docs/notebooks/98_watershed.ipynb | 124 ++++++++++++++++++++++++++++++ docs/tutorials.md | 1 + examples/README.md | 1 + leafmap/common.py | 84 ++++++++++++-------- mkdocs.yml | 1 + 5 files changed, 179 insertions(+), 32 deletions(-) create mode 100644 docs/notebooks/98_watershed.ipynb diff --git a/docs/notebooks/98_watershed.ipynb b/docs/notebooks/98_watershed.ipynb new file mode 100644 index 0000000000..c03ff7a508 --- /dev/null +++ b/docs/notebooks/98_watershed.ipynb @@ -0,0 +1,124 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "[![image](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://demo.leafmap.org/lab/index.html?path=notebooks/98_watershed.ipynb)\n", + "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/leafmap/blob/master/docs/notebooks/98_watershed.ipynb)\n", + "[![image](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/opengeos/leafmap/HEAD)\n", + "\n", + "**Retrieving watershed boundaries from the National Hydrography Dataset (NHD)**\n", + "\n", + "Uncomment the following line to install [leafmap](https://leafmap.org) if needed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "# %pip install -U leafmap" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import leafmap" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "point_geometry = {\"x\": -114.452985, \"y\": 36.13323}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = leafmap.get_wbd(point_geometry, digit=8, return_geometry=True)\n", + "gdf.explore()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "bbox_geometry = {\"xmin\": -115.0954, \"ymin\": 36.1000, \"xmax\": -114.6658, \"ymax\": 36.1986}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = leafmap.get_wbd(bbox_geometry, digit=8, return_geometry=True)\n", + "gdf.explore()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "huc_id = \"10160002\"\n", + "gdf = leafmap.get_wbd(searchText=huc_id, digit=8)\n", + "gdf.explore()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = leafmap.get_wbd(searchText=huc_id, digit=10)\n", + "gdf.explore()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "geo", + "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.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorials.md b/docs/tutorials.md index 57440c18d7..016df026a8 100644 --- a/docs/tutorials.md +++ b/docs/tutorials.md @@ -109,6 +109,7 @@ 95. Editing vector data interactively ([notebook](https://leafmap.org/notebooks/95_edit_vector)) 96. Batch editing vector data attributes interactively ([notebook](https://leafmap.org/notebooks/96_batch_edit_vector)) 97. Downloading Overture Maps data ([notebook](https://leafmap.org/notebooks/97_overture_data)) +98. Retrieving watershed boundaries from the National Hydrography Dataset (NHD) ([notebook](https://leafmap.org/notebooks/98_watershed)) ## Demo diff --git a/examples/README.md b/examples/README.md index 6c02488628..73c33d9b8b 100644 --- a/examples/README.md +++ b/examples/README.md @@ -116,6 +116,7 @@ 95. Editing vector data interactively ([notebook](https://leafmap.org/notebooks/95_edit_vector)) 96. Batch editing vector data attributes interactively ([notebook](https://leafmap.org/notebooks/96_batch_edit_vector)) 97. Downloading Overture Maps data ([notebook](https://leafmap.org/notebooks/97_overture_data)) +98. Retrieving watershed boundaries from the National Hydrography Dataset (NHD) ([notebook](https://leafmap.org/notebooks/98_watershed)) ## Demo diff --git a/leafmap/common.py b/leafmap/common.py index 0063084f5f..6140383f44 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -14303,7 +14303,7 @@ def construct_bbox( return geometry -def get_nhd_wbd( +def get_nhd( geometry: Union[ "gpd.GeoDataFrame", str, List[float], Tuple[float, float, float, float] ], @@ -14545,7 +14545,8 @@ def detect_geometry_type(geometry): def get_wbd( - geometry: Union["gpd.GeoDataFrame", Dict[str, Any]], + geometry: Union["gpd.GeoDataFrame", Dict[str, Any]] = None, + searchText: Optional[str] = None, inSR: str = "4326", outSR: str = "3857", digit: int = 8, @@ -14600,6 +14601,14 @@ def detect_geometry_type(geometry): else: raise ValueError("Unsupported geometry type or invalid geometry structure.") + allowed_digit_values = [2, 4, 6, 8, 10, 12, 14, 16] + if digit not in allowed_digit_values: + raise ValueError( + f"Invalid digit value. Allowed values are {allowed_digit_values}" + ) + + layer = allowed_digit_values.index(digit) + 1 + # Convert GeoDataFrame to a dictionary if needed if isinstance(geometry, gpd.GeoDataFrame): geometry_dict = _convert_geodataframe_to_esri_format(geometry)[0] @@ -14609,38 +14618,44 @@ def detect_geometry_type(geometry): geometry_dict = geometry elif isinstance(geometry, str): geometry_dict = geometry - else: + elif searchText is None: raise ValueError( "Invalid geometry input. Must be a GeoDataFrame or a dictionary." ) - - # Convert geometry to a JSON string (required by the API) - if isinstance(geometry_dict, dict): - geometry_json = json.dumps(geometry_dict) else: - geometry_json = geometry_dict - - allowed_digit_values = [2, 4, 6, 8, 10, 12, 14, 16] - if digit not in allowed_digit_values: - raise ValueError( - f"Invalid digit value. Allowed values are {allowed_digit_values}" - ) + geometry_dict = None - layer = allowed_digit_values.index(digit) + 1 - - # API URL for querying the WBD - url = f"https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/{layer}/query" - - # Construct the query parameters - params = { - "geometry": geometry_json, - "geometryType": geometry_type, - "inSR": inSR, - "spatialRel": spatialRel, - "outFields": outFields, - "returnGeometry": str(return_geometry).lower(), - "f": "json", - } + if geometry_dict is not None: + # Convert geometry to a JSON string (required by the API) + if isinstance(geometry_dict, dict): + geometry_json = json.dumps(geometry_dict) + else: + geometry_json = geometry_dict + + # Construct the query parameters + params = { + "geometry": geometry_json, + "geometryType": geometry_type, + "inSR": inSR, + "spatialRel": spatialRel, + "outFields": outFields, + "returnGeometry": str(return_geometry).lower(), + "f": "json", + } + # API URL for querying the WBD + url = f"https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/{layer}/query" + else: + # Construct the query parameters + params = { + "searchText": searchText, + "contains": "true", + "layers": str(layer), + "inSR": inSR, + "outFields": outFields, + "returnGeometry": str(return_geometry).lower(), + "f": "json", + } + url = f"https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/find" # Add additional keyword arguments for key, value in kwargs.items(): @@ -14654,8 +14669,13 @@ def detect_geometry_type(geometry): data = response.json() - # Extract features from the API response - features = data.get("features", []) + if geometry_dict is not None: + # Extract features from the API response + features = data.get("features", []) + crs = f"EPSG:{data['spatialReference']['latestWkid']}" + else: + features = data.get("results", []) + crs = f"EPSG:{data['results'][0]['geometry']['spatialReference']['latestWkid']}" # Prepare attribute data and geometries attributes = [feature["attributes"] for feature in features] @@ -14678,7 +14698,7 @@ def detect_geometry_type(geometry): gdf = gpd.GeoDataFrame( df, geometry=geometries, - crs=f"EPSG:{data['spatialReference']['latestWkid']}", + crs=crs, ) if outSR != "3857": gdf = gdf.to_crs(outSR) diff --git a/mkdocs.yml b/mkdocs.yml index 40ffd6f7e3..b29bc5e6f8 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -322,3 +322,4 @@ nav: - notebooks/95_edit_vector.ipynb - notebooks/96_batch_edit_vector.ipynb - notebooks/97_overture_data.ipynb + - notebooks/98_watershed.ipynb From 632270e64e01673b74fed64ac24dd89e2edc7446 Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Sat, 28 Sep 2024 16:24:10 -0400 Subject: [PATCH 6/6] Add add_nwi method --- docs/notebooks/99_wetlands.ipynb | 241 +++++++++++++++++++++++++++++++ docs/tutorials.md | 1 + examples/README.md | 1 + leafmap/common.py | 6 +- leafmap/foliumap.py | 59 ++++++++ leafmap/leafmap.py | 59 ++++++++ mkdocs.yml | 1 + 7 files changed, 367 insertions(+), 1 deletion(-) create mode 100644 docs/notebooks/99_wetlands.ipynb diff --git a/docs/notebooks/99_wetlands.ipynb b/docs/notebooks/99_wetlands.ipynb new file mode 100644 index 0000000000..45b98fb47f --- /dev/null +++ b/docs/notebooks/99_wetlands.ipynb @@ -0,0 +1,241 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "[![image](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://demo.leafmap.org/lab/index.html?path=notebooks/99_wetlands.ipynb)\n", + "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/leafmap/blob/master/docs/notebooks/99_wetlands.ipynb)\n", + "[![image](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/opengeos/leafmap/HEAD)\n", + "\n", + "**Retrieving wetland boundaries from the National Wetlands Inventory (NWI)**\n", + "\n", + "The [National Wetlands Inventory (NWI)](https://www.fws.gov/program/national-wetlands-inventory/wetlands-data) is a comprehensive geospatial inventory of U.S. wetlands. It is a publicly available resource that provides detailed information on the abundance, characteristics, and distribution of wetlands in the United States. The NWI dataset is maintained by the U.S. Fish and Wildlife Service (USFWS) and is used by a wide range of stakeholders, including federal, state, and local agencies, researchers, and conservation organizations.\n", + "\n", + "The notebook demonstrates how to retrieve wetland boundaries from the NWI using the `leafmap` Python package. The NWI dataset is available as a web service, which allows users to access the data programmatically. The `leafmap` package provides a simple interface for querying the NWI web service and visualizing the wetland boundaries on an interactive map.\n", + "\n", + "Uncomment the following line to install [leafmap](https://leafmap.org) if needed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "# %pip install -U leafmap" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import leafmap" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map(center=[47.223940, -99.885386], zoom=14)\n", + "m.add_basemap(\"HYBRID\")\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Using point geometry" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "point_geometry = {\"x\": -99.907986, \"y\": 47.216359}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = leafmap.get_nwi(point_geometry)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map(center=[47.223940, -99.885386], zoom=14)\n", + "m.add_basemap(\"HYBRID\")\n", + "m.add_nwi(gdf, add_legend=True)\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Using bounding box" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "bbox_geometry = {\"xmin\": -99.9023, \"ymin\": 47.211, \"xmax\": -99.8556, \"ymax\": 47.2325}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = leafmap.get_nwi(bbox_geometry)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map()\n", + "m.add_basemap(\"HYBRID\")\n", + "m.add_nwi(gdf, layer_name=\"Wetlands\", zoom_to_layer=True)\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## Using GeoDataFrame" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "bbox = [-99.8653, 47.3952, -99.7498, 47.4441]\n", + "bbox_geometry = leafmap.bbox_to_gdf(bbox)\n", + "bbox_geometry" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = leafmap.get_nwi(bbox_geometry)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map()\n", + "m.add_basemap(\"HYBRID\")\n", + "m.add_nwi(gdf, layer_name=\"Wetlands\", zoom_to_layer=True)\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "## Using HUC8 watershed boundary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "huc8 = \"11120104\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = leafmap.get_nwi_by_huc8(huc8, quiet=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map()\n", + "m.add_basemap(\"HYBRID\")\n", + "m.add_nwi(gdf, layer_name=\"Wetlands\", zoom_to_layer=True)\n", + "m" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "geo", + "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.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorials.md b/docs/tutorials.md index 016df026a8..bd556709bb 100644 --- a/docs/tutorials.md +++ b/docs/tutorials.md @@ -110,6 +110,7 @@ 96. Batch editing vector data attributes interactively ([notebook](https://leafmap.org/notebooks/96_batch_edit_vector)) 97. Downloading Overture Maps data ([notebook](https://leafmap.org/notebooks/97_overture_data)) 98. Retrieving watershed boundaries from the National Hydrography Dataset (NHD) ([notebook](https://leafmap.org/notebooks/98_watershed)) +99. Retrieving wetland boundaries from the National Wetlands Inventory (NWI) ([notebook](https://leafmap.org/notebooks/99_wetlands)) ## Demo diff --git a/examples/README.md b/examples/README.md index 73c33d9b8b..a05d16820c 100644 --- a/examples/README.md +++ b/examples/README.md @@ -117,6 +117,7 @@ 96. Batch editing vector data attributes interactively ([notebook](https://leafmap.org/notebooks/96_batch_edit_vector)) 97. Downloading Overture Maps data ([notebook](https://leafmap.org/notebooks/97_overture_data)) 98. Retrieving watershed boundaries from the National Hydrography Dataset (NHD) ([notebook](https://leafmap.org/notebooks/98_watershed)) +99. Retrieving wetland boundaries from the National Wetlands Inventory (NWI) ([notebook](https://leafmap.org/notebooks/99_wetlands)) ## Demo diff --git a/leafmap/common.py b/leafmap/common.py index 6140383f44..9df27f63c7 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -14521,7 +14521,11 @@ def detect_geometry_type(geometry): # Create a DataFrame for attributes df = pd.DataFrame(attributes) df.rename( - columns={"Shape__Length": "Shape_Length", "Shape__Area": "Shape_Area"}, + columns={ + "Shape__Length": "Shape_Length", + "Shape__Area": "Shape_Area", + "WETLAND_TYPE": "WETLAND_TY", + }, inplace=True, ) diff --git a/leafmap/foliumap.py b/leafmap/foliumap.py index f33e986bea..000e5117f2 100644 --- a/leafmap/foliumap.py +++ b/leafmap/foliumap.py @@ -3118,6 +3118,65 @@ def oam_search( else: print("No images found.") + def add_nwi( + self, + data: Union[str, "gpd.GeoDataFrame"], + add_legend: bool = True, + style_callback: Optional[Callable[[dict], dict]] = None, + layer_name: str = "Wetlands", + **kwargs, + ) -> None: + """ + Adds National Wetlands Inventory (NWI) data to the map. + + Args: + data (Union[str, gpd.GeoDataFrame]): The NWI data to add. It can be a file path or a GeoDataFrame. + add_legend (bool): Whether to add a legend to the map. Defaults to True. + style_callback (Optional[Callable[[dict], dict]]): A callback function to style the features. Defaults to None. + layer_name (str): The name of the layer to add. Defaults to "Wetlands". + **kwargs: Additional keyword arguments to pass to the add_vector or add_gdf method. + + Returns: + None + """ + + nwi = { + "Freshwater Forested/Shrub Wetland": "#008837", + "Freshwater Emergent Wetland": "#7fc31c", + "Freshwater Pond": "#688cc0", + "Estuarine and Marine Wetland": "#66c2a5", + "Riverine": "#0190bf", + "Lake": "#13007c", + "Estuarine and Marine Deepwater": "#007c88", + "Other": "#b28656", + } + + def nwi_color(feature): + return { + "color": "black", + "fillColor": ( + nwi[feature["properties"]["WETLAND_TY"]] + if feature["properties"]["WETLAND_TY"] in nwi + else "gray" + ), + "fillOpacity": 0.6, + "weight": 1, + } + + if style_callback is None: + style_callback = nwi_color + + if isinstance(data, str): + self.add_vector( + data, style_callback=style_callback, layer_name=layer_name, **kwargs + ) + else: + self.add_gdf( + data, style_callback=style_callback, layer_name=layer_name, **kwargs + ) + if add_legend: + self.add_legend(title="Wetland Type", builtin_legend="NWI") + def remove_labels(self, **kwargs): """Removes a layer from the map.""" print("The folium plotting backend does not support removing labels.") diff --git a/leafmap/leafmap.py b/leafmap/leafmap.py index 17ff87a452..e8651c3d8e 100644 --- a/leafmap/leafmap.py +++ b/leafmap/leafmap.py @@ -5570,6 +5570,65 @@ def batch_edit_lines( **kwargs, ) + def add_nwi( + self, + data: Union[str, "gpd.GeoDataFrame"], + add_legend: bool = True, + style_callback: Optional[Callable[[dict], dict]] = None, + layer_name: str = "Wetlands", + **kwargs, + ) -> None: + """ + Adds National Wetlands Inventory (NWI) data to the map. + + Args: + data (Union[str, gpd.GeoDataFrame]): The NWI data to add. It can be a file path or a GeoDataFrame. + add_legend (bool): Whether to add a legend to the map. Defaults to True. + style_callback (Optional[Callable[[dict], dict]]): A callback function to style the features. Defaults to None. + layer_name (str): The name of the layer to add. Defaults to "Wetlands". + **kwargs: Additional keyword arguments to pass to the add_vector or add_gdf method. + + Returns: + None + """ + + nwi = { + "Freshwater Forested/Shrub Wetland": "#008837", + "Freshwater Emergent Wetland": "#7fc31c", + "Freshwater Pond": "#688cc0", + "Estuarine and Marine Wetland": "#66c2a5", + "Riverine": "#0190bf", + "Lake": "#13007c", + "Estuarine and Marine Deepwater": "#007c88", + "Other": "#b28656", + } + + def nwi_color(feature): + return { + "color": "black", + "fillColor": ( + nwi[feature["properties"]["WETLAND_TY"]] + if feature["properties"]["WETLAND_TY"] in nwi + else "gray" + ), + "fillOpacity": 0.6, + "weight": 1, + } + + if style_callback is None: + style_callback = nwi_color + + if isinstance(data, str): + self.add_vector( + data, style_callback=style_callback, layer_name=layer_name, **kwargs + ) + else: + self.add_gdf( + data, style_callback=style_callback, layer_name=layer_name, **kwargs + ) + if add_legend: + self.add_legend(title="Wetland Type", builtin_legend="NWI") + # The functions below are outside the Map class. diff --git a/mkdocs.yml b/mkdocs.yml index b29bc5e6f8..4a59c4fcb8 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -323,3 +323,4 @@ nav: - notebooks/96_batch_edit_vector.ipynb - notebooks/97_overture_data.ipynb - notebooks/98_watershed.ipynb + - notebooks/99_wetlands.ipynb