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 }} 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/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 57440c18d7..bd556709bb 100644 --- a/docs/tutorials.md +++ b/docs/tutorials.md @@ -109,6 +109,8 @@ 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)) +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 6c02488628..a05d16820c 100644 --- a/examples/README.md +++ b/examples/README.md @@ -116,6 +116,8 @@ 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)) +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 4428b87f4e..9df27f63c7 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -14362,3 +14362,423 @@ def get_nhd( gdf = None 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", + 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 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 + # API URL for querying wetlands + url = "https://fwspublicservices.wim.usgs.gov/wetlandsmapservice/rest/services/Wetlands/FeatureServer/0/query" + + # 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", + "WETLAND_TYPE": "WETLAND_TY", + }, + 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 + + +def get_wbd( + geometry: Union["gpd.GeoDataFrame", Dict[str, Any]] = None, + searchText: Optional[str] = None, + 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.") + + 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] + 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 + elif searchText is None: + raise ValueError( + "Invalid geometry input. Must be a GeoDataFrame or a dictionary." + ) + else: + geometry_dict = None + + 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(): + 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() + + 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] + 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=crs, + ) + if outSR != "3857": + gdf = gdf.to_crs(outSR) + + if output is not None: + gdf.to_file(output) + + 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 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 40ffd6f7e3..4a59c4fcb8 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -322,3 +322,5 @@ nav: - notebooks/95_edit_vector.ipynb - notebooks/96_batch_edit_vector.ipynb - notebooks/97_overture_data.ipynb + - notebooks/98_watershed.ipynb + - notebooks/99_wetlands.ipynb