From 7f9544e8b724f122f4c4804d37f6527af5a70959 Mon Sep 17 00:00:00 2001 From: Qiusheng Wu Date: Thu, 26 Sep 2024 00:42:51 -0400 Subject: [PATCH] Add get_nhd function --- leafmap/basemaps.py | 14 +++++-- leafmap/common.py | 98 ++++++++++++++++++++++++++++++++++++++++----- leafmap/leafmap.py | 2 + 3 files changed, 100 insertions(+), 14 deletions(-) diff --git a/leafmap/basemaps.py b/leafmap/basemaps.py index 12b14019d..9a6f7da79 100644 --- a/leafmap/basemaps.py +++ b/leafmap/basemaps.py @@ -92,7 +92,7 @@ # Custom WMS tile services. WMS_TILES = { "FWS NWI Wetlands": { - "url": "https://www.fws.gov/wetlands/arcgis/services/Wetlands/MapServer/WMSServer?", + "url": "https://fwspublicservices.wim.usgs.gov/wetlandsmapservice/services/Wetlands/MapServer/WMSServer?", "layers": "1", "name": "FWS NWI Wetlands", "attribution": "FWS", @@ -100,7 +100,7 @@ "transparent": True, }, "FWS NWI Wetlands Raster": { - "url": "https://www.fws.gov/wetlands/arcgis/services/Wetlands_Raster/ImageServer/WMSServer?", + "url": "https://fwspublicservices.wim.usgs.gov/wetlandsmapservice/services/WetlandsRaster/ImageServer/WMSServer?", "layers": "0", "name": "FWS NWI Wetlands Raster", "attribution": "FWS", @@ -213,12 +213,20 @@ }, "USGS 3DEP Elevation": { "url": "https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?", - "layers": "33DEPElevation:Hillshade Elevation Tinted", + "layers": "3DEPElevation:Hillshade Elevation Tinted", "name": "USGS 3DEP Elevation", "attribution": "USGS", "format": "image/png", "transparent": True, }, + "USGS 3DEP Elevation Index": { + "url": "https://index.nationalmap.gov/arcgis/services/3DEPElevationIndex/MapServer/WMSServer?", + "layers": "30", + "name": "USGS 3DEP Elevation Index", + "attribution": "USGS", + "format": "image/png", + "transparent": True, + }, "ESA WorldCover 2020": { "url": "https://services.terrascope.be/wms/v2", "layers": "WORLDCOVER_2020_MAP", diff --git a/leafmap/common.py b/leafmap/common.py index 94e9ad289..b1850d44a 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -6056,7 +6056,7 @@ def image_metadata(image, **kwargs): _, client = get_local_tile_layer(image, return_client=True, **kwargs) else: client = image - return client.metadata() + return client.metadata def image_bandcount(image, **kwargs): @@ -9850,14 +9850,10 @@ def get_3dep_dem( try: import py3dep except ImportError: - raise ImportError("py3dep is not installed. Install it with pip install py3dep") + print("py3dep is not installed. Installing py3dep...") + install_package("py3dep") - try: - import geopandas as gpd - except ImportError: - raise ImportError( - "geopandas is not installed. Install it with pip install geopandas" - ) + import geopandas as gpd if output is not None and os.path.exists(output) and not overwrite: print(f"File {output} already exists. Set overwrite=True to overwrite it") @@ -9872,14 +9868,12 @@ def get_3dep_dem( if output is not None: if not output.endswith(".tif"): output += ".tif" - print(output) dem.rio.to_raster(output, **kwargs) if to_cog: image_to_cog(output, output) - else: - return dem + return dem def vector_set_crs(source, output=None, crs="EPSG:4326", **kwargs): @@ -14263,3 +14257,85 @@ def get_overture_data( gdf.to_file(output) return gdf + + +def construct_bbox( + *args: Union[float, Tuple[float, float, float, float]], + buffer: float = 0.001, + crs: str = "EPSG:4326", + return_gdf: bool = False, +) -> Union["Polygon", "gpd.GeoDataFrame"]: + """ + Construct a bounding box (bbox) geometry based on either a centroid point or bbox. + + Args: + *args: Coordinates for the geometry. + - If 2 arguments are provided, it is interpreted as a centroid (x, y) with a buffer. + - If 4 arguments are provided, it is interpreted as a bbox (minx, miny, maxx, maxy). + buffer (float): The buffer distance around the centroid point (default is 0.01 degrees). + crs (str): The coordinate reference system (default is "EPSG:4326"). + return_gdf (bool): Whether to return a GeoDataFrame (default is False). + + Returns: + shapely.geometry.Polygon: The constructed bounding box (Polygon). + """ + from shapely.geometry import box + + if len(args) == 2: + # Case 1: Create a bounding box around the centroid point with a buffer + x, y = args + minx, miny = x - buffer, y - buffer + maxx, maxy = x + buffer, y + buffer + geometry = box(minx, miny, maxx, maxy) + + elif len(args) == 4: + # Case 2: Create a bounding box directly from the given coordinates + geometry = box(args[0], args[1], args[2], args[3]) + + else: + raise ValueError( + "Provide either 2 arguments for centroid (x, y) or 4 arguments for bbox (minx, miny, maxx, maxy)." + ) + + if return_gdf: + return gpd.GeoDataFrame(geometry=[geometry], columns=["geometry"], crs=crs) + else: + return geometry + + +def get_nhd( + geometry, + geo_crs=4326, + xy=True, + buffer=0.001, + dataset="wbd08", + predicate="intersects", + sort_attr: str | None = None, + **kwargs, +): + try: + import pynhd + except ImportError: + print("The pynhd package is required for this function. Installing...") + install_package("pynhd") + + import geopandas as gpd + from pynhd import WaterData + + if isinstance(geometry, (list, tuple)): + crs = f"EPSG:{geo_crs}" + geometry = construct_bbox(*geometry, buffer=buffer, crs=crs, return_gdf=False) + elif isinstance(geometry, gpd.GeoDataFrame): + geometry = geometry.unary_union + elif isinstance(geometry, str): + geometry = gpd.read_file(geometry).unary_union + + water_data = WaterData(dataset) + + try: + gdf = water_data.bygeom(geometry, geo_crs, xy, predicate, sort_attr, **kwargs) + except Exception as e: + print(e) + gdf = None + + return gdf diff --git a/leafmap/leafmap.py b/leafmap/leafmap.py index df80a971e..17ff87a45 100644 --- a/leafmap/leafmap.py +++ b/leafmap/leafmap.py @@ -370,6 +370,8 @@ def add_basemap(self, basemap="HYBRID", show=True, **kwargs) -> None: elif basemap in basemaps and basemaps[basemap].name not in layer_names: self.add(basemap) self.layers[-1].visible = show + for param in kwargs: + setattr(self.layers[-1], param, kwargs[param]) arc_add_layer(basemaps[basemap].url, basemap) elif basemap in basemaps and basemaps[basemap].name in layer_names: print(f"{basemap} has been already added before.")