diff --git a/docs/maplibre/overview.md b/docs/maplibre/overview.md index 4cc2ce4cb2..ba35282f86 100644 --- a/docs/maplibre/overview.md +++ b/docs/maplibre/overview.md @@ -681,6 +681,12 @@ Add an external Web Map Service raster layer to the map using addSource's tiles [![](https://i.imgur.com/itFOq8z.png)](https://leafmap.org/maplibre/wms_source) +## Export GeoTIFF from WMS + +Export a GeoTIFF image from a Web Map Service (WMS) within a specified bounding box and scale using the `wms_to_geotiff` function. + +[![](https://github.com/user-attachments/assets/01c24828-4535-4352-8398-894955a6ca90)](https://leafmap.org/maplibre/wms_to_geotiff) + ## Fit to the bounds of a LineString Get the bounds of a LineString. diff --git a/docs/maplibre/wms_to_geotiff.ipynb b/docs/maplibre/wms_to_geotiff.ipynb new file mode 100644 index 0000000000..b14fcd50b9 --- /dev/null +++ b/docs/maplibre/wms_to_geotiff.ipynb @@ -0,0 +1,222 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/leafmap/blob/master/docs/maplibre/wms_to_geotiff.ipynb)\n", + "[![image](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/opengeos/leafmap/HEAD)\n", + "\n", + "**Export GeoTIFF from WMS**\n", + "\n", + "This notebook demonstrates how to export a GeoTIFF image from a Web Map Service (WMS) within a specified bounding box and scale using the `wms_to_geotiff` function.\n", + "\n", + "Uncomment the following line to install [leafmap](https://leafmap.org) if needed.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# %pip install \"leafmap[maplibre]\" owslib rasterio" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import leafmap.common as common\n", + "import leafmap.maplibregl as leafmap" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Explore Available WMS Layers\n", + "\n", + "First, let's explore the available layers from the New Jersey Natural 2015 imagery WMS service.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "url = \"https://img.nj.gov/imagerywms/Natural2015\"\n", + "layers = common.get_wms_layers(url)\n", + "print(f\"Available layers: {layers}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize the WMS Layer on a Map\n", + "\n", + "Let's visualize the WMS layer on an interactive map to help identify the area of interest.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m = leafmap.Map(center=[-74.5447, 40.6892], zoom=10, style=\"positron\")\n", + "m.add_wms_layer(\n", + " url, layers=\"Natural2015\", name=\"NJ Natural 2015\", before_id=m.first_symbol_layer_id\n", + ")\n", + "m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Export GeoTIFF from WMS\n", + "\n", + "Now let's export a GeoTIFF from the WMS service. We'll specify:\n", + "- `url`: The WMS service URL\n", + "- `layers`: The layer(s) to request\n", + "- `bbox`: Bounding box as [minx, miny, maxx, maxy] in EPSG:4326\n", + "- `output`: The output file path\n", + "- `scale`: The resolution in degrees per pixel (for EPSG:4326)\n", + "\n", + "### Example 1: Export using scale parameter\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define bounding box [minx, miny, maxx, maxy] - an area in New Jersey\n", + "bbox = [-74.1834231, 40.6934163, -74.168899, 40.6996788]\n", + "\n", + "# Export GeoTIFF with scale (degrees per pixel)\n", + "output = \"nj_natural_2015_scale.tif\"\n", + "common.wms_to_geotiff(\n", + " url=url,\n", + " layers=\"Natural2015\",\n", + " bbox=bbox,\n", + " output=output,\n", + " scale=0.00001, # ~1 meters per pixel at this latitude\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Example 2: Export using explicit width and height\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Export GeoTIFF with explicit dimensions\n", + "output2 = \"nj_natural_2015_dims.tif\"\n", + "common.wms_to_geotiff(\n", + " url=url,\n", + " layers=\"Natural2015\",\n", + " bbox=bbox,\n", + " output=output2,\n", + " width=1452,\n", + " height=626,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize the Exported GeoTIFF\n", + "\n", + "Let's visualize the exported GeoTIFF on the map to verify the export.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m2 = leafmap.Map(style=\"positron\")\n", + "m2.add_raster(output, layer_name=\"Exported GeoTIFF\")\n", + "m2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Export with Different CRS\n", + "\n", + "You can also export in different coordinate reference systems. Here's an example using EPSG:3857 (Web Mercator).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# You can provide bbox in EPSG:4326 and output in EPSG:3857\n", + "# The function will automatically transform the coordinates\n", + "output3 = \"nj_natural_2015_3857.tif\"\n", + "common.wms_to_geotiff(\n", + " url=url,\n", + " layers=\"Natural2015\",\n", + " bbox=bbox, # bbox in EPSG:4326\n", + " output=output3,\n", + " scale=1, # 1 meters per pixel for projected CRS\n", + " bbox_crs=\"EPSG:4326\", # Input bbox CRS (default)\n", + " output_crs=\"EPSG:3857\", # Output CRS\n", + " quiet=True,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "m3 = leafmap.Map(style=\"positron\")\n", + "m3.add_raster(output3, layer_name=\"Exported GeoTIFF\")\n", + "m3" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "leafmap", + "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.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/leafmap/common.py b/leafmap/common.py index c251e0a69d..bd6e706dd0 100644 --- a/leafmap/common.py +++ b/leafmap/common.py @@ -2769,6 +2769,273 @@ def get_wms_layers(url): return layers +def wms_to_geotiff( + url: str, + layers: str, + bbox: List[float], + output: str, + width: Optional[int] = None, + height: Optional[int] = None, + scale: Optional[float] = None, + bbox_crs: str = "EPSG:4326", + output_crs: Optional[str] = None, + format: str = "image/tiff", + version: str = "1.1.1", + styles: str = "", + transparent: bool = True, + quiet: bool = False, + **kwargs: Any, +) -> str: + """Export a GeoTIFF from a WMS service within a specified bounding box and scale. + + Args: + url (str): The URL of the WMS service. + layers (str): Comma-separated list of WMS layers to request. + bbox (List[float]): Bounding box in the format [minx, miny, maxx, maxy]. + The coordinates should be in the CRS specified by 'bbox_crs'. + output (str): The output file path for the GeoTIFF. + width (int, optional): Width of the output image in pixels. If not provided, + will be calculated from scale. Either width/height or scale must be provided. + height (int, optional): Height of the output image in pixels. If not provided, + will be calculated from scale. Either width/height or scale must be provided. + scale (float, optional): The scale (resolution) in units of the output CRS per pixel. + For EPSG:4326, this is degrees per pixel. For projected CRS like EPSG:3857, + this is meters per pixel. Either width/height or scale must be provided. + bbox_crs (str, optional): The coordinate reference system of the input bbox. + Defaults to "EPSG:4326". + output_crs (str, optional): The coordinate reference system for the output GeoTIFF + and WMS request. If None, uses the same CRS as bbox_crs. Defaults to None. + format (str, optional): The image format to request from WMS. Defaults to "image/tiff". + Some WMS servers may support "image/geotiff" or only "image/png". + version (str, optional): WMS version. Defaults to "1.1.1". + styles (str, optional): WMS styles parameter. Defaults to "". + transparent (bool, optional): Request transparent background. Defaults to True. + quiet (bool, optional): Suppress output messages. Defaults to False. + **kwargs: Additional keyword arguments to pass to the WMS GetMap request. + + Returns: + str: The output file path. + + Raises: + ValueError: If neither width/height nor scale is provided. + ImportError: If required packages are not installed. + Exception: If the WMS request fails. + + Example: + >>> url = "https://img.nj.gov/imagerywms/Natural2015" + >>> layers = "Natural2015" + >>> bbox = [-74.6, 40.5, -74.4, 40.7] # [minx, miny, maxx, maxy] in EPSG:4326 + >>> wms_to_geotiff(url, layers, bbox, "output.tif", scale=0.0001) + >>> # Export with different output CRS + >>> wms_to_geotiff(url, layers, bbox, "output_3857.tif", scale=10, output_crs="EPSG:3857") + """ + import io + + try: + from owslib.wms import WebMapService + except ImportError: + raise ImportError("Please install owslib using 'pip install owslib'.") + + try: + import rasterio + from rasterio.transform import from_bounds + except ImportError: + raise ImportError("Please install rasterio using 'pip install rasterio'.") + + try: + from PIL import Image + except ImportError: + raise ImportError("Please install Pillow using 'pip install Pillow'.") + + # Set output_crs to bbox_crs if not provided + if output_crs is None: + output_crs = bbox_crs + + # Validate inputs + if len(bbox) != 4: + raise ValueError("bbox must be a list of 4 values: [minx, miny, maxx, maxy]") + + # Ensure bbox values are numeric and ordered correctly: minx < maxx, miny < maxy + try: + minx, miny, maxx, maxy = map(float, bbox) + except (TypeError, ValueError) as exc: + raise ValueError( + "bbox coordinates must be numeric values: [minx, miny, maxx, maxy]" + ) from exc + + if not (minx < maxx and miny < maxy): + raise ValueError( + "bbox coordinates must be in order: minx < maxx and miny < maxy" + ) + # Transform bbox if bbox_crs differs from output_crs + if bbox_crs.upper() != output_crs.upper(): + try: + from pyproj import Transformer + except ImportError: + raise ImportError( + "Please install pyproj using 'pip install pyproj' for CRS transformation." + ) + + transformer = Transformer.from_crs(bbox_crs, output_crs, always_xy=True) + minx, miny = transformer.transform(minx, miny) + maxx, maxy = transformer.transform(maxx, maxy) + + if not quiet: + print( + f"Transformed bbox from {bbox_crs} to {output_crs}: [{minx}, {miny}, {maxx}, {maxy}]" + ) + + bbox_width = maxx - minx + bbox_height = maxy - miny + + if width is None and height is None and scale is None: + raise ValueError( + "Either width/height or scale must be provided to determine output size." + ) + + # Calculate width/height from scale if not provided + if scale is not None: + if width is None: + width = int(bbox_width / scale) + if height is None: + height = int(bbox_height / scale) + + # Ensure minimum size + width = max(1, width) + height = max(1, height) + + # Some servers have maximum size limits, cap at reasonable values + max_size = 8192 + if width > max_size or height > max_size: + ratio = min(max_size / width, max_size / height) + width = int(width * ratio) + height = int(height * ratio) + if not quiet: + print(f"Image size capped to {width}x{height} due to server limits.") + + if not quiet: + print(f"Requesting WMS image: {width}x{height} pixels") + print(f"Bounding box ({output_crs}): [{minx}, {miny}, {maxx}, {maxy}]") + print(f"Scale: {scale if scale else 'calculated from dimensions'}") + + # Connect to WMS service + wms = WebMapService(url, version=version) + + # Determine correct coordinate reference system parameter name based on WMS version + crs_param = "srs" + if version is not None: + try: + version_parts = tuple(int(part) for part in str(version).split(".")) + if version_parts >= (1, 3, 0): + crs_param = "crs" + except Exception: + # Fallback: simple string prefix check if parsing fails + if str(version).startswith("1.3"): + crs_param = "crs" + + # Build GetMap request parameters + request_params = { + "layers": [layers] if "," not in layers else layers.split(","), + crs_param: output_crs, + "bbox": (minx, miny, maxx, maxy), + "size": (width, height), + "format": format, + "styles": [styles] if styles else [""], + "transparent": transparent, + } + request_params.update(kwargs) + + # Make the WMS GetMap request + try: + response = wms.getmap(**request_params) + except Exception as e: + raise Exception(f"WMS GetMap request failed: {str(e)}") + + # Read the response + img_data = response.read() + + # Check if the response is a valid image + if b"ServiceException" in img_data or b"Exception" in img_data[:1000]: + raise Exception( + f"WMS service returned an exception: {img_data[:500].decode('utf-8', errors='ignore')}" + ) + + output = check_file_path(output) + + # Determine if response is already a GeoTIFF or needs conversion + is_geotiff = format.lower() in ["image/tiff", "image/geotiff", "image/tif"] + + if is_geotiff: + # Try to save directly as GeoTIFF + try: + # First, check if it's a valid GeoTIFF with georeferencing + with rasterio.open(io.BytesIO(img_data)) as src: + if src.crs is not None: + # Valid GeoTIFF with CRS, write it directly + output_temp = temp_file_path(extension=".tif") + with open(output_temp, "wb") as f: + f.write(img_data) + image_to_cog(output_temp, output, nodata=None) + try: + os.remove(output_temp) + except Exception: + pass + if not quiet: + print(f"GeoTIFF saved to {output}") + return output + except Exception: + pass # Fall through to manual georeferencing + + # Convert image to GeoTIFF with georeferencing + try: + # Open as PIL Image + img = Image.open(io.BytesIO(img_data)) + img_array = np.array(img) + except Exception as e: + raise Exception(f"Failed to process image data: {str(e)}") + + # Handle different image modes + if img_array.ndim == 2: + # Grayscale + count = 1 + img_array = img_array.reshape((1, img_array.shape[0], img_array.shape[1])) + elif img_array.ndim == 3: + # RGB or RGBA + count = 3 + img_array = np.transpose(img_array, (2, 0, 1)) + else: + raise ValueError(f"Unexpected image array dimensions: {img_array.ndim}") + + # Create the transform + transform = from_bounds(minx, miny, maxx, maxy, width, height) + + # Write the GeoTIFF + profile = { + "driver": "GTiff", + "dtype": img_array.dtype, + "width": width, + "height": height, + "count": count, + "crs": output_crs, + "transform": transform, + "compress": "deflate", + } + + output_temp = temp_file_path(extension=".tif") + with rasterio.open(output_temp, "w", **profile) as dst: + dst.write(img_array) + + image_to_cog(output_temp, output, nodata=None) + try: + os.remove(output_temp) + except Exception: + pass + if not quiet: + print(f"GeoTIFF saved to {output}") + + return output + + def streamlit_legend(html, width=None, height=None, scrolling=True): """Streamlit function to display a legend. @@ -4017,7 +4284,9 @@ def image_to_geotiff(image, dst_path, dtype=None, to_cog=True, **kwargs) -> None image_to_cog(dst_path, dst_path) -def image_to_cog(source, dst_path=None, profile="deflate", BIGTIFF=None, **kwargs: Any): +def image_to_cog( + source, dst_path=None, profile="deflate", BIGTIFF=None, nodata=None, **kwargs: Any +): """Converts an image to a COG file. Args: @@ -4025,6 +4294,9 @@ def image_to_cog(source, dst_path=None, profile="deflate", BIGTIFF=None, **kwarg dst_path (str, optional): An output dataset path or or PathLike object. Defaults to None. profile (str, optional): COG profile. More at https://cogeotiff.github.io/rio-cogeo/profile. Defaults to "deflate". BIGTIFF (str, optional): Create a BigTIFF file. Can be "IF_SAFER" or "YES". Defaults to None. + nodata (int, float, or None, optional): NoData value for the output COG. Set to None + to prevent any value from being treated as NoData (avoids 255 being changed to 254 + in uint8 images). Defaults to None. Raises: ImportError: If rio-cogeo is not installed. @@ -4059,7 +4331,7 @@ def image_to_cog(source, dst_path=None, profile="deflate", BIGTIFF=None, **kwarg if BIGTIFF is not None: dst_profile.update({"BIGTIFF": BIGTIFF}) - cog_translate(source, dst_path, dst_profile, **kwargs) + cog_translate(source, dst_path, dst_profile, nodata=nodata, **kwargs) def cog_validate(source, verbose=False) -> Tuple[bool, List[str], List[str]]: diff --git a/mkdocs.yml b/mkdocs.yml index 02ec27058b..4dc1a9005b 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -295,6 +295,7 @@ nav: - maplibre/video_on_a_map.ipynb - maplibre/visualize_population_density.ipynb - maplibre/wms_source.ipynb + - maplibre/wms_to_geotiff.ipynb - maplibre/zoom_to_linestring.ipynb - Notebooks: - notebooks/00_key_features.ipynb