From 6f320a58a27f14f979770d842a20ef5e8d80babc Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 10:31:34 -0500 Subject: [PATCH 01/11] Adds NRW DTM provider for geographic data processing Introduces NRWProvider class for handling North Rhine-Westphalia DGM1 data. Integrates data download, merge, reproject, and convert functionalities. Enhances terrain data capabilities with community-driven features. --- maps4fs/__init__.py | 1 + maps4fs/generator/dtm/nrw.py | 363 +++++++++++++++++++++++++++++++++++ 2 files changed, 364 insertions(+) create mode 100644 maps4fs/generator/dtm/nrw.py diff --git a/maps4fs/__init__.py b/maps4fs/__init__.py index d0d8e756..a5707d66 100644 --- a/maps4fs/__init__.py +++ b/maps4fs/__init__.py @@ -2,6 +2,7 @@ from maps4fs.generator.dtm.dtm import DTMProvider from maps4fs.generator.dtm.srtm import SRTM30Provider from maps4fs.generator.dtm.usgs import USGSProvider +from maps4fs.generator.dtm.nrw import NRWProvider from maps4fs.generator.game import Game from maps4fs.generator.map import Map from maps4fs.generator.settings import ( diff --git a/maps4fs/generator/dtm/nrw.py b/maps4fs/generator/dtm/nrw.py new file mode 100644 index 00000000..35ad4561 --- /dev/null +++ b/maps4fs/generator/dtm/nrw.py @@ -0,0 +1,363 @@ +"""This module contains provider of USGS data.""" + +import os +from datetime import datetime +from zipfile import ZipFile + +import numpy as np +from pyproj import Transformer +import rasterio +import requests +from rasterio.enums import Resampling +from rasterio.merge import merge +from rasterio.warp import calculate_default_transform, reproject +from rasterio.windows import from_bounds +from owslib.wcs import WebCoverageService +from owslib.util import Authentication + +from maps4fs.generator.dtm.dtm import DTMProvider, DTMProviderSettings + + +class NRWProviderSettings(DTMProviderSettings): + """Settings for the USGS provider.""" + + max_local_elevation: int = 255 + # wcs_url: str = '' + + +class NRWProvider(DTMProvider): + """Generic provider of WCS sources.""" + + _code = "NRW" + _name = "North Rhine-Westphalia DGM1" + _region = "DE" + _icon = "🇩🇪" + _resolution = 1 + _data: np.ndarray | None = None + _settings = NRWProviderSettings + _author = "[kbrandwijk](https://github.com/kbrandwijk)" + _is_community = True + _instructions = ( + "ℹ️ Set the max local elevation to approx the local max elevation for your area in" + " meters. This will allow you to use heightScale 255 in GE with minimal tweaking." + " Setting this value too low can cause a flat map!" + ) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + timestamp = datetime.now().strftime("%Y-%m-%d_%H-%M-%S") + self.shared_tiff_path = os.path.join(self._tile_directory, "shared") + os.makedirs(self.shared_tiff_path, exist_ok=True) + self.output_path = os.path.join(self._tile_directory, f"timestamp_{timestamp}") + os.makedirs(self.output_path, exist_ok=True) + + def download_tif_files(self, urls: list[str]) -> list[str]: + """Download GeoTIFF files from the given URLs. + + Arguments: + urls (list): List of URLs to download GeoTIFF files from. + + Returns: + list: List of paths to the downloaded GeoTIFF files. + """ + tif_files = [] + for url in urls: + file_name = os.path.basename(url) + self.logger.debug("Retrieving TIFF: %s", file_name) + file_path = os.path.join(self.shared_tiff_path, file_name) + if not os.path.exists(file_path): + try: + # Send a GET request to the file URL + response = requests.get(url, stream=True) # pylint: disable=W3101 + response.raise_for_status() # Raise an error for HTTP status codes 4xx/5xx + + # Write the content of the response to the file + with open(file_path, "wb") as file: + for chunk in response.iter_content(chunk_size=8192): # Download in chunks + file.write(chunk) + self.logger.info("File downloaded successfully: %s", file_path) + if file_name.endswith('.zip'): + with ZipFile(file_path, "r") as f_in: + f_in.extract(file_name.replace('.zip', '.img'), self.shared_tiff_path) + tif_files.append(file_path.replace('.zip', '.img')) + else: + tif_files.append(file_path) + except requests.exceptions.RequestException as e: + self.logger.error("Failed to download file: %s", e) + else: + self.logger.debug("File already exists: %s", file_name) + if file_name.endswith('.zip'): + if not os.path.exists(file_path.replace('.zip', '.img')): + with ZipFile(file_path, "r") as f_in: + f_in.extract(file_name.replace('.zip', '.img'), self.shared_tiff_path) + tif_files.append(file_path.replace('.zip', '.img')) + else: + tif_files.append(file_path) + + return tif_files + + def merge_geotiff(self, input_files: list[str], output_file: str) -> None: + """Merge multiple GeoTIFF files into a single GeoTIFF file. + + Arguments: + input_files (list): List of input GeoTIFF files to merge. + output_file (str): Path to save the merged GeoTIFF file. + """ + # Open all input GeoTIFF files as datasets + self.logger.debug("Merging tiff files...") + datasets = [rasterio.open(file) for file in input_files] + + # Merge datasets + mosaic, out_transform = merge(datasets, nodata=0) + + # Get metadata from the first file and update it for the output + out_meta = datasets[0].meta.copy() + out_meta.update( + { + "driver": "GTiff", + "height": mosaic.shape[1], + "width": mosaic.shape[2], + "transform": out_transform, + "count": mosaic.shape[0], # Number of bands + } + ) + + # Write merged GeoTIFF to the output file + with rasterio.open(output_file, "w", **out_meta) as dest: + dest.write(mosaic) + + self.logger.debug("GeoTIFF images merged successfully into %s", output_file) + + def reproject_geotiff(self, input_tiff: str, output_tiff: str, target_crs: str) -> None: + """Reproject a GeoTIFF file to a new coordinate reference system (CRS). + + Arguments: + input_tiff (str): Path to the input GeoTIFF file. + output_tiff (str): Path to save the reprojected GeoTIFF file. + target_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). + """ + # Open the source GeoTIFF + self.logger.debug("Reprojecting GeoTIFF to %s CRS...", target_crs) + with rasterio.open(input_tiff) as src: + # Get the transform, width, and height of the target CRS + transform, width, height = calculate_default_transform( + src.crs, target_crs, src.width, src.height, *src.bounds + ) + + # Update the metadata for the target GeoTIFF + kwargs = src.meta.copy() + kwargs.update( + {"crs": target_crs, "transform": transform, "width": width, "height": height} + ) + + # Open the destination GeoTIFF file and reproject + with rasterio.open(output_tiff, "w", **kwargs) as dst: + for i in range(1, src.count + 1): # Iterate over all raster bands + reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs, + dst_transform=transform, + dst_crs=target_crs, + resampling=Resampling.nearest, # Choose resampling method + ) + self.logger.debug("Reprojected GeoTIFF saved to %s", output_tiff) + + def extract_roi(self, input_tiff: str) -> np.ndarray: # pylint: disable=W0237 + """ + Crop a GeoTIFF based on given geographic bounding box and save to a new file. + + Arguments: + input_tiff (str): Path to the input GeoTIFF file. + + Returns: + np.ndarray: Numpy array of the cropped GeoTIFF. + """ + self.logger.debug("Extracting ROI...") + # Open the input GeoTIFF + with rasterio.open(input_tiff) as src: + + # Create a rasterio window from the bounding box + (north, south, east, west) = self.get_bbox() + window = from_bounds(west, south, east, north, transform=src.transform) + + data = src.read(1, window=window) + self.logger.debug("Extracted ROI") + return data + + # pylint: disable=R0914, R0917, R0913 + def convert_geotiff_to_geotiff( + self, + input_tiff: str, + output_tiff: str, + min_height: float, + max_height: float, + target_crs: str, + ) -> None: + """ + Convert a GeoTIFF to a scaled GeoTIFF with UInt16 values using a specific coordinate + system and output size. + + Arguments: + input_tiff (str): Path to the input GeoTIFF file. + output_tiff (str): Path to save the output GeoTIFF file. + min_height (float): Minimum terrain height (input range). + max_height (float): Maximum terrain height (input range). + target_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). + """ + # Open the input GeoTIFF file + self.logger.debug("Converting to uint16") + with rasterio.open(input_tiff) as src: + # Ensure the input CRS matches the target CRS (reprojection may be required) + if str(src.crs) != str(target_crs): + raise ValueError( + f"The GeoTIFF CRS is {src.crs}, but the target CRS is {target_crs}. " + "Reprojection may be required." + ) + + # Read the data from the first band + data = src.read(1) # Assuming the input GeoTIFF has only a single band + + # Identify the input file's NoData value + input_nodata = src.nodata + if input_nodata is None: + input_nodata = -999999.0 # Default fallback if no NoData value is defined + nodata_value = 0 + # Replace NoData values (e.g., -999999.0) with the new NoData value + # (e.g., 65535 for UInt16) + data[data == input_nodata] = nodata_value + + # Scale the data to the 0–65535 range (UInt16), avoiding NoData areas + scaled_data = np.clip( + (data - min_height) * (65535 / (max_height - min_height)), 0, 65535 + ).astype(np.uint16) + scaled_data[data == nodata_value] = ( + nodata_value # Preserve NoData value in the scaled array + ) + + # Compute the proper transform to ensure consistency + # Get the original transform, width, and height + transform = src.transform + width = src.width + height = src.height + left, bottom, right, top = src.bounds + + # Adjust the transform matrix to make sure bounds and transform align correctly + transform = rasterio.transform.from_bounds(left, bottom, right, top, width, height) + + # Prepare metadata for the output GeoTIFF + metadata = src.meta.copy() + metadata.update( + { + "dtype": rasterio.uint16, # Update dtype for uint16 + "crs": target_crs, # Update CRS if needed + "nodata": nodata_value, # Set the new NoData value + "transform": transform, # Use the updated, consistent transform + } + ) + + # Write the scaled data to the output GeoTIFF + with rasterio.open(output_tiff, "w", **metadata) as dst: + dst.write(scaled_data, 1) # Write the first band + + self.logger.debug( + "GeoTIFF successfully converted and saved to %s, with nodata value: %s.", + output_tiff, + nodata_value, + ) + + def generate_data(self) -> np.ndarray: + """Generate data from the NRW provider. + + Returns: + np.ndarray: Numpy array of the data. + """ + wcs = WebCoverageService('https://www.wcs.nrw.de/geobasis/wcs_nw_dgm', auth=Authentication(verify=False), timeout=600) + (north, south, east, west) = self.get_bbox() + + # Because the target CRS is rotated compared to 4326, we need to transform all corners and get the bounding box that covers the rotated area + transformer = Transformer.from_crs("epsg:4326", "epsg:25832") + bottom_left_x, bottom_left_y = transformer.transform(xx=south, yy=west) + top_left_x, top_left_y = transformer.transform(xx=north, yy=west) + top_right_x, top_right_y = transformer.transform(xx=north, yy=east) + bottom_right_x, bottom_right_y = transformer.transform(xx=south, yy=east) + + # get the bounding box that covers the rotated area + # coordinate system is reversed, so this looks weird but it is correct + west = min(bottom_left_y, bottom_right_y) + east = max(top_left_y, top_right_y) + south = min(bottom_left_x, top_left_x) + north = max(bottom_right_x, top_right_x) + + print(west, east, south, north) + + # Generate coordinates for the tile corners + x_coords = np.arange(west, east, 1000) + y_coords = np.arange(south, north, 1000) + + # Append the edge coordinates because arange does not include the end value + x_coords = np.append(x_coords, east) + y_coords = np.append(y_coords, north) + + print(x_coords) + print(y_coords) + + # Use meshgrid to create grid combinations + x_min, y_min = np.meshgrid(x_coords[:-1], y_coords[:-1], indexing="ij") + x_max, y_max = np.meshgrid(x_coords[1:], y_coords[1:], indexing="ij") + + # Combine into a single array of bounding boxes + tiles = np.stack([x_min.ravel(), y_min.ravel(), x_max.ravel(), y_max.ravel()], axis=1) + + print(tiles) + + all_tif_files = [] + + for tile in tiles: + file_name = '_'.join(map(str, tile)) + '.tif' + file_path = os.path.join(self.shared_tiff_path, file_name) + + if not os.path.exists(file_path): + output = wcs.getCoverage( + identifier=['nw_dgm'], + subsets=[('y', str(tile[0]), str(tile[2])), ('x', str(tile[1]), str(tile[3]))], + format='image/tiff' + ) + with open(file_path, 'wb') as f: + f.write(output.read()) + + all_tif_files.append(file_path) + + self.merge_geotiff(all_tif_files, os.path.join(self.output_path, "merged.tif")) + + self.reproject_geotiff( + os.path.join(self.output_path, "merged.tif"), + os.path.join(self.output_path, "reprojected.tif"), + "EPSG:4326", + ) + self.convert_geotiff_to_geotiff( + os.path.join(self.output_path, "reprojected.tif"), + os.path.join(self.output_path, "translated.tif"), + min_height=0, + max_height=self.user_settings.max_local_elevation, # type: ignore + target_crs="EPSG:4326", + ) + return self.extract_roi(os.path.join(self.output_path, "translated.tif")) + + def get_numpy(self) -> np.ndarray: + """Get numpy array of the tile. + + Returns: + np.ndarray: Numpy array of the tile. + """ + if not self.user_settings: + raise ValueError("user_settings is 'none'") + if self.user_settings.max_local_elevation <= 0: # type: ignore + raise ValueError( + "Entered 'max_local_elevation' value is unable to be used. " + "Use a value greater than 0." + ) + if not self._data: + self._data = self.generate_data() + return self._data From d336a81c20341be38d62468a9e19d9c5f0d3dab5 Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 11:32:21 -0500 Subject: [PATCH 02/11] fix type stuff --- maps4fs/generator/dtm/nrw.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/maps4fs/generator/dtm/nrw.py b/maps4fs/generator/dtm/nrw.py index 35ad4561..a0f31e69 100644 --- a/maps4fs/generator/dtm/nrw.py +++ b/maps4fs/generator/dtm/nrw.py @@ -274,7 +274,7 @@ def generate_data(self) -> np.ndarray: np.ndarray: Numpy array of the data. """ wcs = WebCoverageService('https://www.wcs.nrw.de/geobasis/wcs_nw_dgm', auth=Authentication(verify=False), timeout=600) - (north, south, east, west) = self.get_bbox() + north, south, east, west = self.get_bbox() # Because the target CRS is rotated compared to 4326, we need to transform all corners and get the bounding box that covers the rotated area transformer = Transformer.from_crs("epsg:4326", "epsg:25832") @@ -297,8 +297,8 @@ def generate_data(self) -> np.ndarray: y_coords = np.arange(south, north, 1000) # Append the edge coordinates because arange does not include the end value - x_coords = np.append(x_coords, east) - y_coords = np.append(y_coords, north) + x_coords = np.append(x_coords, east).astype(x_coords.dtype) + y_coords = np.append(y_coords, north).astype(y_coords.dtype) print(x_coords) print(y_coords) From 5f83e02f9173c631d46b7c37c400e349c008baed Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 12:47:58 -0500 Subject: [PATCH 03/11] some code cleanup --- maps4fs/generator/dtm/nrw.py | 69 ++++++++++++++++++++++++------------ 1 file changed, 47 insertions(+), 22 deletions(-) diff --git a/maps4fs/generator/dtm/nrw.py b/maps4fs/generator/dtm/nrw.py index a0f31e69..66272757 100644 --- a/maps4fs/generator/dtm/nrw.py +++ b/maps4fs/generator/dtm/nrw.py @@ -267,53 +267,63 @@ def convert_geotiff_to_geotiff( nodata_value, ) - def generate_data(self) -> np.ndarray: - """Generate data from the NRW provider. + def transform_bbox(self, bbox: tuple[float, float, float, float], to_crs: str) -> tuple[float, float, float, float]: + """Transform the bounding box to a different coordinate reference system (CRS). + + Arguments: + bbox (tuple): Bounding box to transform (north, south, east, west). + to_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). Returns: - np.ndarray: Numpy array of the data. + tuple: Transformed bounding box (north, south, east, west). """ - wcs = WebCoverageService('https://www.wcs.nrw.de/geobasis/wcs_nw_dgm', auth=Authentication(verify=False), timeout=600) - north, south, east, west = self.get_bbox() - - # Because the target CRS is rotated compared to 4326, we need to transform all corners and get the bounding box that covers the rotated area - transformer = Transformer.from_crs("epsg:4326", "epsg:25832") + transformer = Transformer.from_crs("epsg:4326", to_crs) + north, south, east, west = bbox bottom_left_x, bottom_left_y = transformer.transform(xx=south, yy=west) top_left_x, top_left_y = transformer.transform(xx=north, yy=west) top_right_x, top_right_y = transformer.transform(xx=north, yy=east) bottom_right_x, bottom_right_y = transformer.transform(xx=south, yy=east) - # get the bounding box that covers the rotated area - # coordinate system is reversed, so this looks weird but it is correct west = min(bottom_left_y, bottom_right_y) east = max(top_left_y, top_right_y) south = min(bottom_left_x, top_left_x) north = max(bottom_right_x, top_right_x) - print(west, east, south, north) + return north, south, east, west - # Generate coordinates for the tile corners - x_coords = np.arange(west, east, 1000) - y_coords = np.arange(south, north, 1000) + def tile_bbox(self, bbox: tuple[float, float, float, float], tile_size: int) -> list[tuple[float, float, float, float]]: + """Tile the bounding box into smaller bounding boxes of a specified size. - # Append the edge coordinates because arange does not include the end value + Arguments: + bbox (tuple): Bounding box to tile (north, south, east, west). + tile_size (int): Size of the tiles in meters. + + Returns: + list: List of smaller bounding boxes (north, south, east, west). + """ + north, south, east, west = bbox + x_coords = np.arange(west, east, tile_size) + y_coords = np.arange(south, north, tile_size) x_coords = np.append(x_coords, east).astype(x_coords.dtype) y_coords = np.append(y_coords, north).astype(y_coords.dtype) - print(x_coords) - print(y_coords) - - # Use meshgrid to create grid combinations x_min, y_min = np.meshgrid(x_coords[:-1], y_coords[:-1], indexing="ij") x_max, y_max = np.meshgrid(x_coords[1:], y_coords[1:], indexing="ij") - # Combine into a single array of bounding boxes tiles = np.stack([x_min.ravel(), y_min.ravel(), x_max.ravel(), y_max.ravel()], axis=1) - print(tiles) + return tiles - all_tif_files = [] + def download_tiles(self, tiles: list[tuple[float, float, float, float]]) -> list[str]: + """Download tiles from the NRW provider. + + Arguments: + tiles (list): List of tiles to download. + Returns: + list: List of paths to the downloaded GeoTIFF files. + """ + all_tif_files = [] for tile in tiles: file_name = '_'.join(map(str, tile)) + '.tif' file_path = os.path.join(self.shared_tiff_path, file_name) @@ -328,6 +338,21 @@ def generate_data(self) -> np.ndarray: f.write(output.read()) all_tif_files.append(file_path) + return all_tif_files + + def generate_data(self) -> np.ndarray: + """Generate data from the NRW provider. + + Returns: + np.ndarray: Numpy array of the data. + """ + north, south, east, west = self.get_bbox() + + north, south, east, west = self.transform_bbox((north, south, east, west), "epsg:25832") + + tiles = self.tile_bbox((north, south, east, west), 1000) + + all_tif_files = self.download_tiles(tiles) self.merge_geotiff(all_tif_files, os.path.join(self.output_path, "merged.tif")) From 1925fd037964ab934124b8a0663db4cad5f64e9f Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 17:27:30 -0500 Subject: [PATCH 04/11] fix providers using settings with the same name --- webui/generator.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/webui/generator.py b/webui/generator.py index 251d8d9c..d75436cb 100644 --- a/webui/generator.py +++ b/webui/generator.py @@ -173,7 +173,7 @@ def get_settings(self): field_name = self.snake_to_human(raw_field_name) disabled = self.is_disabled_on_public(raw_field_name) st.write(getattr(Settings, raw_field_name.upper())) - widget = self._create_widget(field_name, raw_field_name, field_value, disabled) + widget = self._create_widget('main', field_name, raw_field_name, field_value, disabled) category[raw_field_name] = widget @@ -182,11 +182,12 @@ def get_settings(self): self.settings = settings def _create_widget( - self, field_name: str, raw_field_name: str, value: int | bool | str, disabled: bool = False - ) -> int | bool: + self, prefix: str, field_name: str, raw_field_name: str, value: int | bool | str | tuple | dict, disabled: bool = False + ) -> int | bool | str: """Create a widget for the given field. Arguments: + prefix (str): The prefix for the key, used to make sure the key is unique across different contexts (e.g. providers). field_name (str): The field name. raw_field_name (str): The raw field name. value (int | bool): The value of the field. @@ -195,22 +196,23 @@ def _create_widget( Returns: int | bool: The widget for the field. """ + key = f"{prefix}_{raw_field_name}" if disabled: st.warning(Messages.SETTING_DISABLED_ON_PUBLIC.format(setting=field_name)) if type(value) is int: return st.number_input( - label=field_name, value=value, min_value=0, key=raw_field_name, disabled=disabled + label=field_name, value=value, min_value=0, key=key, disabled=disabled ) elif type(value) is bool: - return st.checkbox(label=field_name, value=value, key=raw_field_name, disabled=disabled) + return st.checkbox(label=field_name, value=value, key=key, disabled=disabled) elif type(value) is tuple: - return st.selectbox(label=field_name, options=value) + return st.selectbox(label=field_name, key=key, options=value) elif type(value) is dict: return st.selectbox( label=field_name, options=value, format_func=value.get, - key=raw_field_name, + key=key, disabled=disabled, ) else: @@ -252,7 +254,7 @@ def provider_info(self) -> None: settings_json = provider_settings.model_dump() for raw_field_name, value in settings_json.items(): field_name = self.snake_to_human(raw_field_name) - widget = self._create_widget(field_name, raw_field_name, value) + widget = self._create_widget(provider_code, field_name, raw_field_name, value) settings[raw_field_name] = widget From 9531deb65b834d1f630ddacb12b7d84f02037259 Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 22:29:59 -0500 Subject: [PATCH 05/11] refactor DTM providers --- maps4fs/generator/dtm/dtm.py | 323 ++++++++++++++++++++++++++++----- maps4fs/generator/dtm/nrw.py | 329 ++++------------------------------ maps4fs/generator/dtm/srtm.py | 220 ++++++----------------- maps4fs/generator/dtm/usgs.py | 216 +--------------------- webui/generator.py | 8 +- 5 files changed, 377 insertions(+), 719 deletions(-) diff --git a/maps4fs/generator/dtm/dtm.py b/maps4fs/generator/dtm/dtm.py index e1113d9f..22b0cd7e 100644 --- a/maps4fs/generator/dtm/dtm.py +++ b/maps4fs/generator/dtm/dtm.py @@ -11,7 +11,10 @@ import numpy as np import osmnx as ox # type: ignore import rasterio # type: ignore -import requests +from rasterio.warp import calculate_default_transform, reproject +from rasterio.enums import Resampling +from rasterio.merge import merge + from pydantic import BaseModel from maps4fs.logger import Logger @@ -23,6 +26,9 @@ class DTMProviderSettings(BaseModel): """Base class for DTM provider settings models.""" + easy_mode: bool = True + power_factor: int = 0 + # pylint: disable=too-many-public-methods class DTMProvider(ABC): @@ -44,6 +50,20 @@ class DTMProvider(ABC): _instructions: str | None = None + _base_instructions = ( + "ℹ️ Using **Easy mode** is recommended, as it automatically adjusts the values in the image, " + "so the terrain elevation in Giants Editor will match real world elevation in meters. \n" + "ℹ️ If the terrain height difference in the real world is bigger than 255 meters, " + "the [Height scale](https://github.com/iwatkot/maps4fs/blob/main/docs/dem.md#height-scale)" + " parameter in the **map.i3d** file will be automatically adjusted. \n" + "⚡ If the **Easy mode** option is disabled, you will probably get completely flat " + "terrain, unless you adjust the DEM Multiplier Setting or the Height scale parameter in " + "the Giants Editor. \n" + "💡 You can use the **Power factor** setting to make the difference between heights " + "bigger. Be extremely careful with this setting, and use only low values, otherwise your " + "terrain may be completely broken. \n" + ) + # pylint: disable=R0913, R0917 def __init__( self, @@ -152,7 +172,7 @@ def is_community(cls) -> bool: return cls._is_community @classmethod - def settings(cls) -> Type[DTMProviderSettings] | None: + def settings(cls) -> Type[DTMProviderSettings]: """Settings model of the provider. Returns: @@ -169,6 +189,15 @@ def instructions(cls) -> str | None: """ return cls._instructions + @classmethod + def base_instructions(cls) -> str | None: + """Instructions for using any provider. + + Returns: + str: Instructions for using any provider. + """ + return cls._base_instructions + @property def user_settings(self) -> DTMProviderSettings | None: """User settings of the provider. @@ -212,63 +241,116 @@ def get_provider_descriptions(cls) -> dict[str, str]: """ providers = {} for provider in cls.__subclasses__(): - if not provider.is_base(): - providers[provider._code] = provider.description() # pylint: disable=W0212 + providers[provider._code] = provider.description() # pylint: disable=W0212 return providers # type: ignore - def download_tile(self, output_path: str, **kwargs) -> bool: - """Download a tile from the provider. - - Arguments: - output_path (str): Path to save the downloaded tile. + @abstractmethod + def download_tiles(self) -> list[str]: + """Download tiles from the provider. Returns: - bool: True if the tile was downloaded successfully, False otherwise. + list: List of paths to the downloaded tiles. """ - url = self.formatted_url(**kwargs) - response = requests.get(url, stream=True, timeout=10) - if response.status_code == 200: - with open(output_path, "wb") as file: - for chunk in response.iter_content(chunk_size=1024): - file.write(chunk) - return True - return False - - def get_or_download_tile(self, output_path: str, **kwargs) -> str | None: - """Get or download a tile from the provider. + raise NotImplementedError - Arguments: - output_path (str): Path to save the downloaded tile. + def get_numpy(self) -> np.ndarray: + """Get numpy array of the tile. + Resulting array must be 16 bit (signed or unsigned) integer and it should be already + windowed to the bounding box of ROI. It also must have only one channel. Returns: - str: Path to the downloaded tile or None if the tile not exists and was - not downloaded. + np.ndarray: Numpy array of the tile. """ - if not os.path.exists(output_path): - if not self.download_tile(output_path, **kwargs): - return None - return output_path + # download tiles using DTM provider implementation + tiles = self.download_tiles() + self.logger.debug(f"Downloaded {len(tiles)} DEM tiles") + + # merge tiles if necessary + if len(tiles) > 1: + self.logger.debug("Multiple tiles downloaded. Merging tiles") + tile, _ = self.merge_geotiff(tiles) + else: + tile = tiles[0] + + # determine CRS of the resulting tile and reproject if necessary + with rasterio.open(tile) as src: + crs = src.crs + if crs != "EPSG:4326": + self.logger.debug(f"Reprojecting GeoTIFF from {crs} to EPSG:4326...") + tile = self.reproject_geotiff(tile) + + # extract region of interest from the tile + data = self.extract_roi(tile) + + # process elevation data to be compatible with the game + data = self.process_elevation(data) - def get_tile_parameters(self, *args, **kwargs) -> dict: - """Get parameters for the tile, that will be used to format the URL. - Must be implemented in subclasses. + return data + + def process_elevation(self, data: np.ndarray) -> np.ndarray: + """Process elevation data. + + Arguments: + data (np.ndarray): Elevation data. Returns: - dict: Tile parameters to format the URL. + np.ndarray: Processed elevation data. """ - raise NotImplementedError + self.data_info = {} + self.add_numpy_params(data, "original") - @abstractmethod - def get_numpy(self) -> np.ndarray: - """Get numpy array of the tile. - Resulting array must be 16 bit (signed or unsigned) integer and it should be already - windowed to the bounding box of ROI. It also must have only one channel. + data = self.signed_to_unsigned(data) + self.add_numpy_params(data, "grounded") + + original_deviation = int(self.data_info["original_deviation"]) + in_game_maximum_height = 65535 // 255 + if original_deviation > in_game_maximum_height: + suggested_height_scale_multiplier = ( + original_deviation / in_game_maximum_height # type: ignore + ) + suggested_height_scale_value = int(255 * suggested_height_scale_multiplier) + else: + suggested_height_scale_multiplier = 1 + suggested_height_scale_value = 255 + + self.data_info["suggested_height_scale_multiplier"] = suggested_height_scale_multiplier + self.data_info["suggested_height_scale_value"] = suggested_height_scale_value + + self.map.shared_settings.height_scale_multiplier = ( # type: ignore + suggested_height_scale_multiplier + ) + self.map.shared_settings.height_scale_value = suggested_height_scale_value # type: ignore + + if self.user_settings.easy_mode: # type: ignore + try: + data = self.normalize_dem(data) + self.add_numpy_params(data, "normalized") + + normalized_deviation = self.data_info["normalized_deviation"] + z_scaling_factor = normalized_deviation / original_deviation # type: ignore + self.data_info["z_scaling_factor"] = z_scaling_factor + + self.map.shared_settings.mesh_z_scaling_factor = z_scaling_factor # type: ignore + self.map.shared_settings.change_height_scale = True # type: ignore + + except Exception as e: # pylint: disable=W0718 + self.logger.error( + "Failed to normalize DEM data. Error: %s. Using original data.", e + ) + + return data + + def info_sequence(self) -> dict[str, int | str | float] | None: + """Returns the information sequence for the component. Must be implemented in the child + class. If the component does not have an information sequence, an empty dictionary must be + returned. Returns: - np.ndarray: Numpy array of the tile. + dict[str, int | str | float] | None: Information sequence for the component. """ - raise NotImplementedError + return self.data_info + # region helpers def get_bbox(self) -> tuple[float, float, float, float]: """Get bounding box of the tile based on the center point and size. @@ -281,6 +363,83 @@ def get_bbox(self) -> tuple[float, float, float, float]: bbox = north, south, east, west return bbox + def reproject_geotiff(self, input_tiff: str) -> str: + """Reproject a GeoTIFF file to a new coordinate reference system (CRS). + + Arguments: + input_tiff (str): Path to the input GeoTIFF file. + target_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). + + Returns: + str: Path to the reprojected GeoTIFF file. + """ + output_tiff = os.path.join(self._tile_directory, "merged.tif") + + # Open the source GeoTIFF + self.logger.debug("Reprojecting GeoTIFF to %s CRS...", "EPSG:4326") + with rasterio.open(input_tiff) as src: + # Get the transform, width, and height of the target CRS + transform, width, height = calculate_default_transform( + src.crs, "EPSG:4326", src.width, src.height, *src.bounds + ) + + # Update the metadata for the target GeoTIFF + kwargs = src.meta.copy() + kwargs.update( + {"crs": "EPSG:4326", "transform": transform, "width": width, "height": height} + ) + + # Open the destination GeoTIFF file and reproject + with rasterio.open(output_tiff, "w", **kwargs) as dst: + for i in range(1, src.count + 1): # Iterate over all raster bands + reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs, + dst_transform=transform, + dst_crs="EPSG:4326", + resampling=Resampling.nearest, # Choose resampling method + ) + + self.logger.debug("Reprojected GeoTIFF saved to %s", output_tiff) + return output_tiff + + def merge_geotiff(self, input_files: list[str]) -> tuple[str, str]: + """Merge multiple GeoTIFF files into a single GeoTIFF file. + + Arguments: + input_files (list): List of input GeoTIFF files to merge. + output_file (str): Path to save the merged GeoTIFF file. + """ + output_file = os.path.join(self._tile_directory, "merged.tif") + # Open all input GeoTIFF files as datasets + self.logger.debug("Merging tiff files...") + datasets = [rasterio.open(file) for file in input_files] + + # Merge datasets + crs = datasets[0].crs + mosaic, out_transform = merge(datasets, nodata=0) + + # Get metadata from the first file and update it for the output + out_meta = datasets[0].meta.copy() + out_meta.update( + { + "driver": "GTiff", + "height": mosaic.shape[1], + "width": mosaic.shape[2], + "transform": out_transform, + "count": mosaic.shape[0], # Number of bands + } + ) + + # Write merged GeoTIFF to the output file + with rasterio.open(output_file, "w", **out_meta) as dest: + dest.write(mosaic) + + self.logger.debug("GeoTIFF images merged successfully into %s", output_file) + return output_file, crs + def extract_roi(self, tile_path: str) -> np.ndarray: """Extract region of interest (ROI) from the GeoTIFF file. @@ -304,18 +463,86 @@ def extract_roi(self, tile_path: str) -> np.ndarray: window.width, window.height, ) - data = src.read(1, window=window) + data = src.read(1, window=window, masked=True) if not data.size > 0: raise ValueError("No data in the tile.") return data - def info_sequence(self) -> dict[str, int | str | float] | None: - """Returns the information sequence for the component. Must be implemented in the child - class. If the component does not have an information sequence, an empty dictionary must be - returned. + def normalize_dem(self, data: np.ndarray) -> np.ndarray: + """Normalize DEM data to 16-bit unsigned integer using max height from settings. + + Arguments: + data (np.ndarray): DEM data from SRTM file after cropping. Returns: - dict[str, int | str | float] | None: Information sequence for the component. + np.ndarray: Normalized DEM data. """ - return self.data_info + maximum_height = int(data.max()) + minimum_height = int(data.min()) + deviation = maximum_height - minimum_height + self.logger.debug( + "Maximum height: %s. Minimum height: %s. Deviation: %s.", + maximum_height, + minimum_height, + deviation, + ) + self.logger.debug("Number of unique values in original DEM data: %s.", np.unique(data).size) + + adjusted_maximum_height = maximum_height * 255 + adjusted_maximum_height = min(adjusted_maximum_height, 65535) + scaling_factor = adjusted_maximum_height / maximum_height + self.logger.debug( + "Adjusted maximum height: %s. Scaling factor: %s.", + adjusted_maximum_height, + scaling_factor, + ) + + if self.user_settings.power_factor: # type: ignore + power_factor = 1 + self.user_settings.power_factor / 10 # type: ignore + self.logger.debug( + "Applying power factor: %s to the DEM data.", + power_factor, + ) + data = np.power(data, power_factor).astype(np.uint16) + + normalized_data = np.round(data * scaling_factor).astype(np.uint16) + self.logger.debug( + "Normalized data maximum height: %s. Minimum height: %s. Number of unique values: %s.", + normalized_data.max(), + normalized_data.min(), + np.unique(normalized_data).size, + ) + return normalized_data + + def signed_to_unsigned(self, data: np.ndarray, add_one: bool = True) -> np.ndarray: + """Convert signed 16-bit integer to unsigned 16-bit integer. + + Arguments: + data (np.ndarray): DEM data from SRTM file after cropping. + + Returns: + np.ndarray: Unsigned DEM data. + """ + data = data - data.min() + if add_one: + data = data + 1 + return data.astype(np.uint16) + + def add_numpy_params( + self, + data: np.ndarray, + prefix: str, + ) -> None: + """Add numpy array parameters to the data_info dictionary. + + Arguments: + data (np.ndarray): Numpy array of the tile. + prefix (str): Prefix for the parameters. + """ + self.data_info[f"{prefix}_minimum_height"] = int(data.min()) # type: ignore + self.data_info[f"{prefix}_maximum_height"] = int(data.max()) # type: ignore + self.data_info[f"{prefix}_deviation"] = int(data.max() - data.min()) # type: ignore + self.data_info[f"{prefix}_unique_values"] = int(np.unique(data).size) # type: ignore + + # endregion diff --git a/maps4fs/generator/dtm/nrw.py b/maps4fs/generator/dtm/nrw.py index 66272757..f9bfcf4d 100644 --- a/maps4fs/generator/dtm/nrw.py +++ b/maps4fs/generator/dtm/nrw.py @@ -2,18 +2,11 @@ import os from datetime import datetime -from zipfile import ZipFile import numpy as np -from pyproj import Transformer -import rasterio -import requests -from rasterio.enums import Resampling -from rasterio.merge import merge -from rasterio.warp import calculate_default_transform, reproject -from rasterio.windows import from_bounds from owslib.wcs import WebCoverageService from owslib.util import Authentication +from pyproj import Transformer from maps4fs.generator.dtm.dtm import DTMProvider, DTMProviderSettings @@ -21,9 +14,6 @@ class NRWProviderSettings(DTMProviderSettings): """Settings for the USGS provider.""" - max_local_elevation: int = 255 - # wcs_url: str = '' - class NRWProvider(DTMProvider): """Generic provider of WCS sources.""" @@ -37,259 +27,21 @@ class NRWProvider(DTMProvider): _settings = NRWProviderSettings _author = "[kbrandwijk](https://github.com/kbrandwijk)" _is_community = True - _instructions = ( - "ℹ️ Set the max local elevation to approx the local max elevation for your area in" - " meters. This will allow you to use heightScale 255 in GE with minimal tweaking." - " Setting this value too low can cause a flat map!" - ) + _instructions = None def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) - timestamp = datetime.now().strftime("%Y-%m-%d_%H-%M-%S") self.shared_tiff_path = os.path.join(self._tile_directory, "shared") os.makedirs(self.shared_tiff_path, exist_ok=True) - self.output_path = os.path.join(self._tile_directory, f"timestamp_{timestamp}") - os.makedirs(self.output_path, exist_ok=True) - - def download_tif_files(self, urls: list[str]) -> list[str]: - """Download GeoTIFF files from the given URLs. - - Arguments: - urls (list): List of URLs to download GeoTIFF files from. - - Returns: - list: List of paths to the downloaded GeoTIFF files. - """ - tif_files = [] - for url in urls: - file_name = os.path.basename(url) - self.logger.debug("Retrieving TIFF: %s", file_name) - file_path = os.path.join(self.shared_tiff_path, file_name) - if not os.path.exists(file_path): - try: - # Send a GET request to the file URL - response = requests.get(url, stream=True) # pylint: disable=W3101 - response.raise_for_status() # Raise an error for HTTP status codes 4xx/5xx - - # Write the content of the response to the file - with open(file_path, "wb") as file: - for chunk in response.iter_content(chunk_size=8192): # Download in chunks - file.write(chunk) - self.logger.info("File downloaded successfully: %s", file_path) - if file_name.endswith('.zip'): - with ZipFile(file_path, "r") as f_in: - f_in.extract(file_name.replace('.zip', '.img'), self.shared_tiff_path) - tif_files.append(file_path.replace('.zip', '.img')) - else: - tif_files.append(file_path) - except requests.exceptions.RequestException as e: - self.logger.error("Failed to download file: %s", e) - else: - self.logger.debug("File already exists: %s", file_name) - if file_name.endswith('.zip'): - if not os.path.exists(file_path.replace('.zip', '.img')): - with ZipFile(file_path, "r") as f_in: - f_in.extract(file_name.replace('.zip', '.img'), self.shared_tiff_path) - tif_files.append(file_path.replace('.zip', '.img')) - else: - tif_files.append(file_path) - - return tif_files - - def merge_geotiff(self, input_files: list[str], output_file: str) -> None: - """Merge multiple GeoTIFF files into a single GeoTIFF file. - - Arguments: - input_files (list): List of input GeoTIFF files to merge. - output_file (str): Path to save the merged GeoTIFF file. - """ - # Open all input GeoTIFF files as datasets - self.logger.debug("Merging tiff files...") - datasets = [rasterio.open(file) for file in input_files] - - # Merge datasets - mosaic, out_transform = merge(datasets, nodata=0) - - # Get metadata from the first file and update it for the output - out_meta = datasets[0].meta.copy() - out_meta.update( - { - "driver": "GTiff", - "height": mosaic.shape[1], - "width": mosaic.shape[2], - "transform": out_transform, - "count": mosaic.shape[0], # Number of bands - } - ) - - # Write merged GeoTIFF to the output file - with rasterio.open(output_file, "w", **out_meta) as dest: - dest.write(mosaic) - - self.logger.debug("GeoTIFF images merged successfully into %s", output_file) - - def reproject_geotiff(self, input_tiff: str, output_tiff: str, target_crs: str) -> None: - """Reproject a GeoTIFF file to a new coordinate reference system (CRS). - - Arguments: - input_tiff (str): Path to the input GeoTIFF file. - output_tiff (str): Path to save the reprojected GeoTIFF file. - target_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). - """ - # Open the source GeoTIFF - self.logger.debug("Reprojecting GeoTIFF to %s CRS...", target_crs) - with rasterio.open(input_tiff) as src: - # Get the transform, width, and height of the target CRS - transform, width, height = calculate_default_transform( - src.crs, target_crs, src.width, src.height, *src.bounds - ) - - # Update the metadata for the target GeoTIFF - kwargs = src.meta.copy() - kwargs.update( - {"crs": target_crs, "transform": transform, "width": width, "height": height} - ) - - # Open the destination GeoTIFF file and reproject - with rasterio.open(output_tiff, "w", **kwargs) as dst: - for i in range(1, src.count + 1): # Iterate over all raster bands - reproject( - source=rasterio.band(src, i), - destination=rasterio.band(dst, i), - src_transform=src.transform, - src_crs=src.crs, - dst_transform=transform, - dst_crs=target_crs, - resampling=Resampling.nearest, # Choose resampling method - ) - self.logger.debug("Reprojected GeoTIFF saved to %s", output_tiff) - - def extract_roi(self, input_tiff: str) -> np.ndarray: # pylint: disable=W0237 - """ - Crop a GeoTIFF based on given geographic bounding box and save to a new file. - - Arguments: - input_tiff (str): Path to the input GeoTIFF file. - - Returns: - np.ndarray: Numpy array of the cropped GeoTIFF. - """ - self.logger.debug("Extracting ROI...") - # Open the input GeoTIFF - with rasterio.open(input_tiff) as src: - - # Create a rasterio window from the bounding box - (north, south, east, west) = self.get_bbox() - window = from_bounds(west, south, east, north, transform=src.transform) - - data = src.read(1, window=window) - self.logger.debug("Extracted ROI") - return data - - # pylint: disable=R0914, R0917, R0913 - def convert_geotiff_to_geotiff( - self, - input_tiff: str, - output_tiff: str, - min_height: float, - max_height: float, - target_crs: str, - ) -> None: - """ - Convert a GeoTIFF to a scaled GeoTIFF with UInt16 values using a specific coordinate - system and output size. - - Arguments: - input_tiff (str): Path to the input GeoTIFF file. - output_tiff (str): Path to save the output GeoTIFF file. - min_height (float): Minimum terrain height (input range). - max_height (float): Maximum terrain height (input range). - target_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). - """ - # Open the input GeoTIFF file - self.logger.debug("Converting to uint16") - with rasterio.open(input_tiff) as src: - # Ensure the input CRS matches the target CRS (reprojection may be required) - if str(src.crs) != str(target_crs): - raise ValueError( - f"The GeoTIFF CRS is {src.crs}, but the target CRS is {target_crs}. " - "Reprojection may be required." - ) - - # Read the data from the first band - data = src.read(1) # Assuming the input GeoTIFF has only a single band - - # Identify the input file's NoData value - input_nodata = src.nodata - if input_nodata is None: - input_nodata = -999999.0 # Default fallback if no NoData value is defined - nodata_value = 0 - # Replace NoData values (e.g., -999999.0) with the new NoData value - # (e.g., 65535 for UInt16) - data[data == input_nodata] = nodata_value - - # Scale the data to the 0–65535 range (UInt16), avoiding NoData areas - scaled_data = np.clip( - (data - min_height) * (65535 / (max_height - min_height)), 0, 65535 - ).astype(np.uint16) - scaled_data[data == nodata_value] = ( - nodata_value # Preserve NoData value in the scaled array - ) - - # Compute the proper transform to ensure consistency - # Get the original transform, width, and height - transform = src.transform - width = src.width - height = src.height - left, bottom, right, top = src.bounds - - # Adjust the transform matrix to make sure bounds and transform align correctly - transform = rasterio.transform.from_bounds(left, bottom, right, top, width, height) - - # Prepare metadata for the output GeoTIFF - metadata = src.meta.copy() - metadata.update( - { - "dtype": rasterio.uint16, # Update dtype for uint16 - "crs": target_crs, # Update CRS if needed - "nodata": nodata_value, # Set the new NoData value - "transform": transform, # Use the updated, consistent transform - } - ) - # Write the scaled data to the output GeoTIFF - with rasterio.open(output_tiff, "w", **metadata) as dst: - dst.write(scaled_data, 1) # Write the first band + def download_tiles(self) -> list[str]: + bbox = self.get_bbox() + bbox = self.transform_bbox(bbox, "EPSG:25832") - self.logger.debug( - "GeoTIFF successfully converted and saved to %s, with nodata value: %s.", - output_tiff, - nodata_value, - ) + tiles = self.tile_bbox(bbox, 1000) - def transform_bbox(self, bbox: tuple[float, float, float, float], to_crs: str) -> tuple[float, float, float, float]: - """Transform the bounding box to a different coordinate reference system (CRS). - - Arguments: - bbox (tuple): Bounding box to transform (north, south, east, west). - to_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). - - Returns: - tuple: Transformed bounding box (north, south, east, west). - """ - transformer = Transformer.from_crs("epsg:4326", to_crs) - north, south, east, west = bbox - bottom_left_x, bottom_left_y = transformer.transform(xx=south, yy=west) - top_left_x, top_left_y = transformer.transform(xx=north, yy=west) - top_right_x, top_right_y = transformer.transform(xx=north, yy=east) - bottom_right_x, bottom_right_y = transformer.transform(xx=south, yy=east) - - west = min(bottom_left_y, bottom_right_y) - east = max(top_left_y, top_right_y) - south = min(bottom_left_x, top_left_x) - north = max(bottom_right_x, top_right_x) - - return north, south, east, west + all_tif_files = self.download_all_tiles(tiles) + return all_tif_files def tile_bbox(self, bbox: tuple[float, float, float, float], tile_size: int) -> list[tuple[float, float, float, float]]: """Tile the bounding box into smaller bounding boxes of a specified size. @@ -314,7 +66,7 @@ def tile_bbox(self, bbox: tuple[float, float, float, float], tile_size: int) -> return tiles - def download_tiles(self, tiles: list[tuple[float, float, float, float]]) -> list[str]: + def download_all_tiles(self, tiles: list[tuple[float, float, float, float]]) -> list[str]: """Download tiles from the NRW provider. Arguments: @@ -324,6 +76,10 @@ def download_tiles(self, tiles: list[tuple[float, float, float, float]]) -> list list: List of paths to the downloaded GeoTIFF files. """ all_tif_files = [] + wcs = WebCoverageService( + 'https://www.wcs.nrw.de/geobasis/wcs_nw_dgm', + auth=Authentication(verify=False), + timeout=600) for tile in tiles: file_name = '_'.join(map(str, tile)) + '.tif' file_path = os.path.join(self.shared_tiff_path, file_name) @@ -340,49 +96,26 @@ def download_tiles(self, tiles: list[tuple[float, float, float, float]]) -> list all_tif_files.append(file_path) return all_tif_files - def generate_data(self) -> np.ndarray: - """Generate data from the NRW provider. + def transform_bbox(self, bbox: tuple[float, float, float, float], to_crs: str) -> tuple[float, float, float, float]: + """Transform the bounding box to a different coordinate reference system (CRS). + + Arguments: + bbox (tuple): Bounding box to transform (north, south, east, west). + to_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). Returns: - np.ndarray: Numpy array of the data. + tuple: Transformed bounding box (north, south, east, west). """ - north, south, east, west = self.get_bbox() - - north, south, east, west = self.transform_bbox((north, south, east, west), "epsg:25832") - - tiles = self.tile_bbox((north, south, east, west), 1000) - - all_tif_files = self.download_tiles(tiles) - - self.merge_geotiff(all_tif_files, os.path.join(self.output_path, "merged.tif")) - - self.reproject_geotiff( - os.path.join(self.output_path, "merged.tif"), - os.path.join(self.output_path, "reprojected.tif"), - "EPSG:4326", - ) - self.convert_geotiff_to_geotiff( - os.path.join(self.output_path, "reprojected.tif"), - os.path.join(self.output_path, "translated.tif"), - min_height=0, - max_height=self.user_settings.max_local_elevation, # type: ignore - target_crs="EPSG:4326", - ) - return self.extract_roi(os.path.join(self.output_path, "translated.tif")) + transformer = Transformer.from_crs("epsg:4326", to_crs) + north, south, east, west = bbox + bottom_left_x, bottom_left_y = transformer.transform(xx=south, yy=west) + top_left_x, top_left_y = transformer.transform(xx=north, yy=west) + top_right_x, top_right_y = transformer.transform(xx=north, yy=east) + bottom_right_x, bottom_right_y = transformer.transform(xx=south, yy=east) - def get_numpy(self) -> np.ndarray: - """Get numpy array of the tile. + west = min(bottom_left_y, bottom_right_y) + east = max(top_left_y, top_right_y) + south = min(bottom_left_x, top_left_x) + north = max(bottom_right_x, top_right_x) - Returns: - np.ndarray: Numpy array of the tile. - """ - if not self.user_settings: - raise ValueError("user_settings is 'none'") - if self.user_settings.max_local_elevation <= 0: # type: ignore - raise ValueError( - "Entered 'max_local_elevation' value is unable to be used. " - "Use a value greater than 0." - ) - if not self._data: - self._data = self.generate_data() - return self._data + return north, south, east, west diff --git a/maps4fs/generator/dtm/srtm.py b/maps4fs/generator/dtm/srtm.py index 5ac286d9..f322b070 100644 --- a/maps4fs/generator/dtm/srtm.py +++ b/maps4fs/generator/dtm/srtm.py @@ -7,7 +7,7 @@ import os import shutil -import numpy as np +import requests from maps4fs.generator.dtm.dtm import DTMProvider, DTMProviderSettings @@ -15,9 +15,6 @@ class SRTM30ProviderSettings(DTMProviderSettings): """Settings for SRTM 30m provider.""" - easy_mode: bool = True - power_factor: int = 0 - class SRTM30Provider(DTMProvider): """Provider of Shuttle Radar Topography Mission (SRTM) 30m data.""" @@ -32,22 +29,6 @@ class SRTM30Provider(DTMProvider): _author = "[iwatkot](https://github.com/iwatkot)" - _instructions = ( - "ℹ️ If you don't know how to work with DEM data, it is recommended to use the " - "**Easy mode** option. It will automatically change the values in the image, so the " - "terrain will be visible in the Giants Editor. If you're an experienced modder, it's " - "recommended to disable this option and work with the DEM data in a usual way. \n" - "ℹ️ If the terrain height difference in the real world is bigger than 255 meters, " - "the [Height scale](https://github.com/iwatkot/maps4fs/blob/main/docs/dem.md#height-scale)" - " parameter in the **map.i3d** file will be changed automatically. \n" - "⚡ If the **Easy mode** option is disabled, you will probably get completely flat " - "terrain, unless you adjust the DEM Multiplier Setting or the Height scale parameter in " - "the Giants Editor. \n" - "💡 You can use the **Power factor** setting to make the difference between heights " - "bigger. Be extremely careful with this setting, and use only low values, otherwise your " - "terrain may be completely broken. \n" - ) - _settings = SRTM30ProviderSettings def __init__(self, *args, **kwargs): @@ -58,169 +39,86 @@ def __init__(self, *args, **kwargs): os.makedirs(self.gz_directory, exist_ok=True) self.data_info: dict[str, int | str | float] | None = None # type: ignore - def get_tile_parameters(self, *args, **kwargs) -> dict[str, str]: - """Returns latitude band and tile name for SRTM tile from coordinates. + def download_tiles(self): + """Download SRTM tiles.""" + north, south, east, west = self.get_bbox() - Arguments: - lat (float): Latitude. - lon (float): Longitude. + tiles = [] + # Look at each corner of the bbox in case the bbox spans across multiple tiles + for pair in [(north, east), (south, west), (south, east), (north, west)]: + tile_parameters = self.get_tile_parameters(*pair) + tile_name = tile_parameters["tile_name"] + decompressed_tile_path = os.path.join(self.hgt_directory, f"{tile_name}.hgt") - Returns: - dict: Tile parameters. - """ - lat, lon = args + if not os.path.isfile(decompressed_tile_path): + compressed_tile_path = os.path.join(self.gz_directory, f"{tile_name}.hgt.gz") + if not self.get_or_download_tile(compressed_tile_path, **tile_parameters): + raise FileNotFoundError(f"Tile {tile_name} not found.") - tile_latitude = math.floor(lat) - tile_longitude = math.floor(lon) + with gzip.open(compressed_tile_path, "rb") as f_in: + with open(decompressed_tile_path, "wb") as f_out: + shutil.copyfileobj(f_in, f_out) + tiles.append(decompressed_tile_path) - latitude_band = f"N{abs(tile_latitude):02d}" if lat >= 0 else f"S{abs(tile_latitude):02d}" - if lon < 0: - tile_name = f"{latitude_band}W{abs(tile_longitude):03d}" - else: - tile_name = f"{latitude_band}E{abs(tile_longitude):03d}" + return tiles - self.logger.debug( - "Detected tile name: %s for coordinates: lat %s, lon %s.", tile_name, lat, lon - ) - return {"latitude_band": latitude_band, "tile_name": tile_name} + # region provider specific helpers + def download_tile(self, output_path: str, **kwargs) -> bool: + """Download a tile from the provider. - def get_numpy(self) -> np.ndarray: - """Get numpy array of the tile. + Arguments: + output_path (str): Path to save the downloaded tile. Returns: - np.ndarray: Numpy array of the tile. - """ - tile_parameters = self.get_tile_parameters(*self.coordinates) - tile_name = tile_parameters["tile_name"] - decompressed_tile_path = os.path.join(self.hgt_directory, f"{tile_name}.hgt") - - if not os.path.isfile(decompressed_tile_path): - compressed_tile_path = os.path.join(self.gz_directory, f"{tile_name}.hgt.gz") - if not self.get_or_download_tile(compressed_tile_path, **tile_parameters): - raise FileNotFoundError(f"Tile {tile_name} not found.") - - with gzip.open(compressed_tile_path, "rb") as f_in: - with open(decompressed_tile_path, "wb") as f_out: - shutil.copyfileobj(f_in, f_out) - - data = self.extract_roi(decompressed_tile_path) - - self.data_info = {} - self.add_numpy_params(data, "original") - - data = self.signed_to_unsigned(data) - self.add_numpy_params(data, "grounded") - - original_deviation = int(self.data_info["original_deviation"]) - in_game_maximum_height = 65535 // 255 - if original_deviation > in_game_maximum_height: - suggested_height_scale_multiplier = ( - original_deviation / in_game_maximum_height # type: ignore - ) - suggested_height_scale_value = int(255 * suggested_height_scale_multiplier) - else: - suggested_height_scale_multiplier = 1 - suggested_height_scale_value = 255 - - self.data_info["suggested_height_scale_multiplier"] = suggested_height_scale_multiplier - self.data_info["suggested_height_scale_value"] = suggested_height_scale_value - - self.map.shared_settings.height_scale_multiplier = ( # type: ignore - suggested_height_scale_multiplier - ) - self.map.shared_settings.height_scale_value = suggested_height_scale_value # type: ignore - - if self.user_settings.easy_mode: # type: ignore - try: - data = self.normalize_dem(data) - self.add_numpy_params(data, "normalized") - - normalized_deviation = self.data_info["normalized_deviation"] - z_scaling_factor = normalized_deviation / original_deviation # type: ignore - self.data_info["z_scaling_factor"] = z_scaling_factor - - self.map.shared_settings.mesh_z_scaling_factor = z_scaling_factor # type: ignore - self.map.shared_settings.change_height_scale = True # type: ignore - - except Exception as e: # pylint: disable=W0718 - self.logger.error( - "Failed to normalize DEM data. Error: %s. Using original data.", e - ) - - return data - - def add_numpy_params( - self, - data: np.ndarray, - prefix: str, - ) -> None: - """Add numpy array parameters to the data_info dictionary. - - Arguments: - data (np.ndarray): Numpy array of the tile. - prefix (str): Prefix for the parameters. + bool: True if the tile was downloaded successfully, False otherwise. """ - self.data_info[f"{prefix}_minimum_height"] = int(data.min()) # type: ignore - self.data_info[f"{prefix}_maximum_height"] = int(data.max()) # type: ignore - self.data_info[f"{prefix}_deviation"] = int(data.max() - data.min()) # type: ignore - self.data_info[f"{prefix}_unique_values"] = int(np.unique(data).size) # type: ignore - - def signed_to_unsigned(self, data: np.ndarray, add_one: bool = True) -> np.ndarray: - """Convert signed 16-bit integer to unsigned 16-bit integer. + url = self.formatted_url(**kwargs) + response = requests.get(url, stream=True, timeout=10) + if response.status_code == 200: + with open(output_path, "wb") as file: + for chunk in response.iter_content(chunk_size=1024): + file.write(chunk) + return True + return False + + def get_or_download_tile(self, output_path: str, **kwargs) -> str | None: + """Get or download a tile from the provider. Arguments: - data (np.ndarray): DEM data from SRTM file after cropping. + output_path (str): Path to save the downloaded tile. Returns: - np.ndarray: Unsigned DEM data. + str: Path to the downloaded tile or None if the tile not exists and was + not downloaded. """ - data = data - data.min() - if add_one: - data = data + 1 - return data.astype(np.uint16) + if not os.path.exists(output_path): + if not self.download_tile(output_path, **kwargs): + return None + return output_path - def normalize_dem(self, data: np.ndarray) -> np.ndarray: - """Normalize DEM data to 16-bit unsigned integer using max height from settings. + def get_tile_parameters(self, *args, **kwargs) -> dict[str, str]: + """Returns latitude band and tile name for SRTM tile from coordinates. Arguments: - data (np.ndarray): DEM data from SRTM file after cropping. + lat (float): Latitude. + lon (float): Longitude. Returns: - np.ndarray: Normalized DEM data. + dict: Tile parameters. """ - maximum_height = int(data.max()) - minimum_height = int(data.min()) - deviation = maximum_height - minimum_height - self.logger.debug( - "Maximum height: %s. Minimum height: %s. Deviation: %s.", - maximum_height, - minimum_height, - deviation, - ) - self.logger.debug("Number of unique values in original DEM data: %s.", np.unique(data).size) + lat, lon = args - adjusted_maximum_height = maximum_height * 255 - adjusted_maximum_height = min(adjusted_maximum_height, 65535) - scaling_factor = adjusted_maximum_height / maximum_height - self.logger.debug( - "Adjusted maximum height: %s. Scaling factor: %s.", - adjusted_maximum_height, - scaling_factor, - ) + tile_latitude = math.floor(lat) + tile_longitude = math.floor(lon) - if self.user_settings.power_factor: # type: ignore - power_factor = 1 + self.user_settings.power_factor / 10 # type: ignore - self.logger.debug( - "Applying power factor: %s to the DEM data.", - power_factor, - ) - data = np.power(data, power_factor).astype(np.uint16) + latitude_band = f"N{abs(tile_latitude):02d}" if lat >= 0 else f"S{abs(tile_latitude):02d}" + if lon < 0: + tile_name = f"{latitude_band}W{abs(tile_longitude):03d}" + else: + tile_name = f"{latitude_band}E{abs(tile_longitude):03d}" - normalized_data = np.round(data * scaling_factor).astype(np.uint16) self.logger.debug( - "Normalized data maximum height: %s. Minimum height: %s. Number of unique values: %s.", - normalized_data.max(), - normalized_data.min(), - np.unique(normalized_data).size, + "Detected tile name: %s for coordinates: lat %s, lon %s.", tile_name, lat, lon ) - return normalized_data + return {"latitude_band": latitude_band, "tile_name": tile_name} + # endregion diff --git a/maps4fs/generator/dtm/usgs.py b/maps4fs/generator/dtm/usgs.py index 2741e94a..df2449d1 100644 --- a/maps4fs/generator/dtm/usgs.py +++ b/maps4fs/generator/dtm/usgs.py @@ -54,6 +54,11 @@ class USGSProvider(DTMProvider): ) + def download_tiles(self): + download_urls = self.get_download_urls() + all_tif_files = self.download_tif_files(download_urls) + return all_tif_files + def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) timestamp = datetime.now().strftime("%Y-%m-%d_%H-%M-%S") @@ -138,214 +143,3 @@ def download_tif_files(self, urls: list[str]) -> list[str]: tif_files.append(file_path) return tif_files - - def merge_geotiff(self, input_files: list[str], output_file: str) -> None: - """Merge multiple GeoTIFF files into a single GeoTIFF file. - - Arguments: - input_files (list): List of input GeoTIFF files to merge. - output_file (str): Path to save the merged GeoTIFF file. - """ - # Open all input GeoTIFF files as datasets - self.logger.debug("Merging tiff files...") - datasets = [rasterio.open(file) for file in input_files] - - # Merge datasets - mosaic, out_transform = merge(datasets, nodata=0) - - # Get metadata from the first file and update it for the output - out_meta = datasets[0].meta.copy() - out_meta.update( - { - "driver": "GTiff", - "height": mosaic.shape[1], - "width": mosaic.shape[2], - "transform": out_transform, - "count": mosaic.shape[0], # Number of bands - } - ) - - # Write merged GeoTIFF to the output file - with rasterio.open(output_file, "w", **out_meta) as dest: - dest.write(mosaic) - - self.logger.debug("GeoTIFF images merged successfully into %s", output_file) - - def reproject_geotiff(self, input_tiff: str, output_tiff: str, target_crs: str) -> None: - """Reproject a GeoTIFF file to a new coordinate reference system (CRS). - - Arguments: - input_tiff (str): Path to the input GeoTIFF file. - output_tiff (str): Path to save the reprojected GeoTIFF file. - target_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). - """ - # Open the source GeoTIFF - self.logger.debug("Reprojecting GeoTIFF to %s CRS...", target_crs) - with rasterio.open(input_tiff) as src: - # Get the transform, width, and height of the target CRS - transform, width, height = calculate_default_transform( - src.crs, target_crs, src.width, src.height, *src.bounds - ) - - # Update the metadata for the target GeoTIFF - kwargs = src.meta.copy() - kwargs.update( - {"crs": target_crs, "transform": transform, "width": width, "height": height} - ) - - # Open the destination GeoTIFF file and reproject - with rasterio.open(output_tiff, "w", **kwargs) as dst: - for i in range(1, src.count + 1): # Iterate over all raster bands - reproject( - source=rasterio.band(src, i), - destination=rasterio.band(dst, i), - src_transform=src.transform, - src_crs=src.crs, - dst_transform=transform, - dst_crs=target_crs, - resampling=Resampling.nearest, # Choose resampling method - ) - self.logger.debug("Reprojected GeoTIFF saved to %s", output_tiff) - - def extract_roi(self, input_tiff: str) -> np.ndarray: # pylint: disable=W0237 - """ - Crop a GeoTIFF based on given geographic bounding box and save to a new file. - - Arguments: - input_tiff (str): Path to the input GeoTIFF file. - - Returns: - np.ndarray: Numpy array of the cropped GeoTIFF. - """ - self.logger.debug("Extracting ROI...") - # Open the input GeoTIFF - with rasterio.open(input_tiff) as src: - - # Create a rasterio window from the bounding box - (north, south, east, west) = self.get_bbox() - window = from_bounds(west, south, east, north, transform=src.transform) - - data = src.read(1, window=window) - self.logger.debug("Extracted ROI") - return data - - # pylint: disable=R0914, R0917, R0913 - def convert_geotiff_to_geotiff( - self, - input_tiff: str, - output_tiff: str, - min_height: float, - max_height: float, - target_crs: str, - ) -> None: - """ - Convert a GeoTIFF to a scaled GeoTIFF with UInt16 values using a specific coordinate - system and output size. - - Arguments: - input_tiff (str): Path to the input GeoTIFF file. - output_tiff (str): Path to save the output GeoTIFF file. - min_height (float): Minimum terrain height (input range). - max_height (float): Maximum terrain height (input range). - target_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84). - """ - # Open the input GeoTIFF file - self.logger.debug("Converting to uint16") - with rasterio.open(input_tiff) as src: - # Ensure the input CRS matches the target CRS (reprojection may be required) - if str(src.crs) != str(target_crs): - raise ValueError( - f"The GeoTIFF CRS is {src.crs}, but the target CRS is {target_crs}. " - "Reprojection may be required." - ) - - # Read the data from the first band - data = src.read(1) # Assuming the input GeoTIFF has only a single band - - # Identify the input file's NoData value - input_nodata = src.nodata - if input_nodata is None: - input_nodata = -999999.0 # Default fallback if no NoData value is defined - nodata_value = 0 - # Replace NoData values (e.g., -999999.0) with the new NoData value - # (e.g., 65535 for UInt16) - data[data == input_nodata] = nodata_value - - # Scale the data to the 0–65535 range (UInt16), avoiding NoData areas - scaled_data = np.clip( - (data - min_height) * (65535 / (max_height - min_height)), 0, 65535 - ).astype(np.uint16) - scaled_data[data == nodata_value] = ( - nodata_value # Preserve NoData value in the scaled array - ) - - # Compute the proper transform to ensure consistency - # Get the original transform, width, and height - transform = src.transform - width = src.width - height = src.height - left, bottom, right, top = src.bounds - - # Adjust the transform matrix to make sure bounds and transform align correctly - transform = rasterio.transform.from_bounds(left, bottom, right, top, width, height) - - # Prepare metadata for the output GeoTIFF - metadata = src.meta.copy() - metadata.update( - { - "dtype": rasterio.uint16, # Update dtype for uint16 - "crs": target_crs, # Update CRS if needed - "nodata": nodata_value, # Set the new NoData value - "transform": transform, # Use the updated, consistent transform - } - ) - - # Write the scaled data to the output GeoTIFF - with rasterio.open(output_tiff, "w", **metadata) as dst: - dst.write(scaled_data, 1) # Write the first band - - self.logger.debug( - "GeoTIFF successfully converted and saved to %s, with nodata value: %s.", - output_tiff, - nodata_value, - ) - - def generate_data(self) -> np.ndarray: - """Generate data from the USGS 1m provider. - - Returns: - np.ndarray: Numpy array of the data. - """ - download_urls = self.get_download_urls() - all_tif_files = self.download_tif_files(download_urls) - self.merge_geotiff(all_tif_files, os.path.join(self.output_path, "merged.tif")) - self.reproject_geotiff( - os.path.join(self.output_path, "merged.tif"), - os.path.join(self.output_path, "reprojected.tif"), - "EPSG:4326", - ) - self.convert_geotiff_to_geotiff( - os.path.join(self.output_path, "reprojected.tif"), - os.path.join(self.output_path, "translated.tif"), - min_height=0, - max_height=self.user_settings.max_local_elevation, # type: ignore - target_crs="EPSG:4326", - ) - return self.extract_roi(os.path.join(self.output_path, "translated.tif")) - - def get_numpy(self) -> np.ndarray: - """Get numpy array of the tile. - - Returns: - np.ndarray: Numpy array of the tile. - """ - if not self.user_settings: - raise ValueError("user_settings is 'none'") - if self.user_settings.max_local_elevation <= 0: # type: ignore - raise ValueError( - "Entered 'max_local_elevation' value is unable to be used. " - "Use a value greater than 0." - ) - if not self._data: - self._data = self.generate_data() - return self._data diff --git a/webui/generator.py b/webui/generator.py index d75436cb..9dd2342a 100644 --- a/webui/generator.py +++ b/webui/generator.py @@ -233,6 +233,9 @@ def provider_info(self) -> None: provider_code = self.dtm_provider_code provider = mfs.DTMProvider.get_provider_by_code(provider_code) + if provider is None: + return + self.provider_settings = None self.provider_info_container.empty() sleep(0.1) @@ -245,6 +248,9 @@ def provider_info(self) -> None: if provider.contributors() is not None: st.write(f"Contributors: {provider.contributors()}") + if provider.base_instructions() is not None: + st.write(provider.base_instructions()) + if provider.instructions() is not None: st.write(provider.instructions()) @@ -336,7 +342,7 @@ def add_left_widgets(self) -> None: on_change=self.provider_info, ) self.provider_settings = None - self.provider_info_container = st.empty() + self.provider_info_container: st.delta_generator.DeltaGenerator = st.empty() self.provider_info() # Rotation input. From 9889c29caa96ebe425917f3916164918fc69dd3a Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 22:58:38 -0500 Subject: [PATCH 06/11] clean up USGS, fix type error --- maps4fs/generator/dtm/dtm.py | 2 +- maps4fs/generator/dtm/usgs.py | 12 +----------- 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/maps4fs/generator/dtm/dtm.py b/maps4fs/generator/dtm/dtm.py index 22b0cd7e..1bbffbd9 100644 --- a/maps4fs/generator/dtm/dtm.py +++ b/maps4fs/generator/dtm/dtm.py @@ -172,7 +172,7 @@ def is_community(cls) -> bool: return cls._is_community @classmethod - def settings(cls) -> Type[DTMProviderSettings]: + def settings(cls) -> Type[DTMProviderSettings] | None: """Settings model of the provider. Returns: diff --git a/maps4fs/generator/dtm/usgs.py b/maps4fs/generator/dtm/usgs.py index df2449d1..f98befad 100644 --- a/maps4fs/generator/dtm/usgs.py +++ b/maps4fs/generator/dtm/usgs.py @@ -5,12 +5,7 @@ from zipfile import ZipFile import numpy as np -import rasterio import requests -from rasterio.enums import Resampling -from rasterio.merge import merge -from rasterio.warp import calculate_default_transform, reproject -from rasterio.windows import from_bounds from maps4fs.generator.dtm.dtm import DTMProvider, DTMProviderSettings @@ -18,7 +13,6 @@ class USGSProviderSettings(DTMProviderSettings): """Settings for the USGS provider.""" - max_local_elevation: int = 255 dataset: tuple | str = ( 'Digital Elevation Model (DEM) 1 meter', 'Alaska IFSAR 5 meter DEM', @@ -43,11 +37,7 @@ class USGSProvider(DTMProvider): _author = "[ZenJakey](https://github.com/ZenJakey)" _contributors = "[kbrandwijk](https://github.com/kbrandwijk)" _is_community = True - _instructions = ( - "ℹ️ Set the max local elevation to approx the local max elevation for your area in" - " meters. This will allow you to use heightScale 255 in GE with minimal tweaking." - " Setting this value too low can cause a flat map!" - ) + _instructions = None _url = ( "https://tnmaccess.nationalmap.gov/api/v1/products?prodFormats=GeoTIFF,IMG" From 72f6f137cf5492a2ebc101134338aa0bbca876bb Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 23:01:47 -0500 Subject: [PATCH 07/11] add missing requirement, clean up NRW --- maps4fs/generator/dtm/nrw.py | 1 - requirements.txt | 3 ++- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/maps4fs/generator/dtm/nrw.py b/maps4fs/generator/dtm/nrw.py index f9bfcf4d..25d5a5c9 100644 --- a/maps4fs/generator/dtm/nrw.py +++ b/maps4fs/generator/dtm/nrw.py @@ -1,7 +1,6 @@ """This module contains provider of USGS data.""" import os -from datetime import datetime import numpy as np from owslib.wcs import WebCoverageService diff --git a/requirements.txt b/requirements.txt index 4cfb6c5c..c815c319 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,4 +11,5 @@ streamlit-stl==0.0.5 pympler pydantic maps4fs -pygmdl \ No newline at end of file +pygmdl +OWSLib \ No newline at end of file From 9c12bbabb59249df3155fc63e957bcd7b8c01c62 Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 23:08:46 -0500 Subject: [PATCH 08/11] linting --- maps4fs/generator/dtm/dtm.py | 7 ++++--- maps4fs/generator/dtm/nrw.py | 11 +++++++++-- maps4fs/generator/dtm/srtm.py | 2 +- 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/maps4fs/generator/dtm/dtm.py b/maps4fs/generator/dtm/dtm.py index 1bbffbd9..be53428e 100644 --- a/maps4fs/generator/dtm/dtm.py +++ b/maps4fs/generator/dtm/dtm.py @@ -30,7 +30,7 @@ class DTMProviderSettings(BaseModel): power_factor: int = 0 -# pylint: disable=too-many-public-methods +# pylint: disable=too-many-public-methods, too-many-instance-attributes class DTMProvider(ABC): """Base class for DTM providers.""" @@ -51,8 +51,9 @@ class DTMProvider(ABC): _instructions: str | None = None _base_instructions = ( - "ℹ️ Using **Easy mode** is recommended, as it automatically adjusts the values in the image, " - "so the terrain elevation in Giants Editor will match real world elevation in meters. \n" + "ℹ️ Using **Easy mode** is recommended, as it automatically adjusts the values in the " + "image, so the terrain elevation in Giants Editor will match real world " + "elevation in meters. \n" "ℹ️ If the terrain height difference in the real world is bigger than 255 meters, " "the [Height scale](https://github.com/iwatkot/maps4fs/blob/main/docs/dem.md#height-scale)" " parameter in the **map.i3d** file will be automatically adjusted. \n" diff --git a/maps4fs/generator/dtm/nrw.py b/maps4fs/generator/dtm/nrw.py index 25d5a5c9..d8a4929e 100644 --- a/maps4fs/generator/dtm/nrw.py +++ b/maps4fs/generator/dtm/nrw.py @@ -14,6 +14,7 @@ class NRWProviderSettings(DTMProviderSettings): """Settings for the USGS provider.""" +# pylint: disable=too-many-locals class NRWProvider(DTMProvider): """Generic provider of WCS sources.""" @@ -42,7 +43,10 @@ def download_tiles(self) -> list[str]: all_tif_files = self.download_all_tiles(tiles) return all_tif_files - def tile_bbox(self, bbox: tuple[float, float, float, float], tile_size: int) -> list[tuple[float, float, float, float]]: + def tile_bbox( + self, + bbox: tuple[float, float, float, float], + tile_size: int) -> list[tuple[float, float, float, float]]: """Tile the bounding box into smaller bounding boxes of a specified size. Arguments: @@ -95,7 +99,10 @@ def download_all_tiles(self, tiles: list[tuple[float, float, float, float]]) -> all_tif_files.append(file_path) return all_tif_files - def transform_bbox(self, bbox: tuple[float, float, float, float], to_crs: str) -> tuple[float, float, float, float]: + def transform_bbox( + self, + bbox: tuple[float, float, float, float], + to_crs: str) -> tuple[float, float, float, float]: """Transform the bounding box to a different coordinate reference system (CRS). Arguments: diff --git a/maps4fs/generator/dtm/srtm.py b/maps4fs/generator/dtm/srtm.py index f322b070..6725abb6 100644 --- a/maps4fs/generator/dtm/srtm.py +++ b/maps4fs/generator/dtm/srtm.py @@ -96,7 +96,7 @@ def get_or_download_tile(self, output_path: str, **kwargs) -> str | None: return None return output_path - def get_tile_parameters(self, *args, **kwargs) -> dict[str, str]: + def get_tile_parameters(self, *args) -> dict[str, str]: """Returns latitude band and tile name for SRTM tile from coordinates. Arguments: From aad0faa0affc52643c4c0d87d5750237652d7d83 Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 23:13:25 -0500 Subject: [PATCH 09/11] trying to fix pylint import error on GHA --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index c815c319..8be4a14c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -12,4 +12,4 @@ pympler pydantic maps4fs pygmdl -OWSLib \ No newline at end of file +owslib \ No newline at end of file From 7f2aa9c0f62eda340589791c78dce30a66e11645 Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Fri, 10 Jan 2025 23:14:50 -0500 Subject: [PATCH 10/11] trying to fix requirements install again --- dev/requirements.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dev/requirements.txt b/dev/requirements.txt index 12cbe635..4223d900 100644 --- a/dev/requirements.txt +++ b/dev/requirements.txt @@ -17,4 +17,5 @@ pympler pydoc-markdown streamlit-stl==0.0.5 pydantic -pygmdl \ No newline at end of file +pygmdl +owslib \ No newline at end of file From f06688a0f12d54a4254372a24b834f68a1e5c513 Mon Sep 17 00:00:00 2001 From: Kim Brandwijk <kim.brandwijk@gmail.com> Date: Sat, 11 Jan 2025 11:13:24 -0500 Subject: [PATCH 11/11] updated docs --- docs/dtm_providers.md | 118 ++++++++++++++++++------------------------ 1 file changed, 50 insertions(+), 68 deletions(-) diff --git a/docs/dtm_providers.md b/docs/dtm_providers.md index af2bf9e7..83b5df23 100644 --- a/docs/dtm_providers.md +++ b/docs/dtm_providers.md @@ -15,19 +15,29 @@ For obvious reasons, in our case we need the DTM and the DSM won't be useful, be ### What a DTM provider does? -A DTM provider is a service that provides elevation data for a given location. The generator will use this data to create a dem image for the map. While it's plenty of DTM providers available, only the ones that provide a free and open access to their data can be used by the generator. +A DTM provider is a service that provides elevation data for a given location. The generator will use this data to create a dem image for the map. While there's plenty of DTM providers available, only the ones that provide a free and open access to their data can be used by the generator. + +The base provider class, [DTMProvider](../mapsfs/generator/dtm/dtm.py) that all DTM providers inherit from, is responsible for all processing of DEM data. Individual DTM providers are responsible only for downloading the DTM tile(s) for the map area. + +The process for generating the elevation data is: + +- Download all DTM tiles for the desired map area (implemented by each DTM provider) +- If the DTM provider downloaded multiple tiles, merge these tiles into one +- If the tile uses a different projection, reproject it to EPSG:4326, which is used for all other data (like OSM) +- Extract the map area from the tile (some providers, like SRTM, return big tiles that are larger than just the desired area) +- Process the elevation data in the tile (optionally, using Easy Mode, the terrain is moved to ground level, so it appears at z=0 in Giants Editor, and processed so the elevation corresponds to real world meters using the default heightScale in your map, if the elevation of your terrain is more than 255 meters, the heightScale property of your map will automatically be adjusted to fit this elevation) ### How to implement a DTM provider? -So the DTM provider is a simple class, that receives coordinate of the center point, the size of the region of interest and should return a 16-bit single channeled numpy array with the elevation data. The elevation data should be in meters above the sea level. +So the DTM provider is a simple class, that receives coordinate of the center point, the size of the region of interest and should download all the needed DTM tiles and return a list of downloaded tiles for further processing by the base class. ### Example of a DTM provider -➡️ Base class and existing providers can be found in the [dtm.py](../maps4fs/generator/dtm.py) file. +➡️ Base class and existing providers can be found in the [dtm](../maps4fs/generator/dtm) folder. -Let's take a look at an example of a DTM provider implementation. +Let's take a look at an example of a DTM provider implementation. -**Step 1:** define description of the provider. +**Step 1:** define description of the provider. ```python class SRTM30Provider(DTMProvider): @@ -47,7 +57,7 @@ class SRTM30Provider(DTMProvider): _instructions = "When working with SRTM provider..." ``` -So, we inherit from the `DTMProvider` class, add some properties to identify the Provider (such as code and region). The most important part is the `_url` property, which is a template for the URL to download the elevation data. But if your provider uses some other approach, you can reimplement related methods. +So, we inherit from the `DTMProvider` class, add some properties to identify the Provider (such as code and region). The most important part is the `_url` property, which is a template for the URL to download the elevation data. But if your provider uses some other approach, you can reimplement related methods. Also, please provide MD-formatted author information, where in [] will be the name of the author and in () will be the link to the author's GitHub profile (or any other profile if you wish). @@ -55,7 +65,7 @@ Please, set the `_is_community` property to `True`, it means that it was develop If you want some message to be displayed when the user selects your provider, you can set the `_instructions` property. -**Step 3 (optional):** use the `DTMProviderSetting` class to define your own settings (if needed). +**Step 2 (optional):** use the `DTMProviderSetting` class to define your own settings (if needed). ```python class SRTM30ProviderSettings(DTMProviderSettings): @@ -65,7 +75,7 @@ class SRTM30ProviderSettings(DTMProviderSettings): input_something: int = 255 ``` -Also, you will need to add a new `_settings` property to the provider class. +Also, you will need to add a new `_settings` property to the provider class. ```python class SRTM30Provider(DTMProvider): @@ -80,65 +90,36 @@ enable_something = self.user_settings.enable_something input_something = self.user_settings.input_something ``` -**Step 3:** implement the `get_tile_parameters` method. - -```python - def get_tile_parameters(self, *args, **kwargs) -> dict[str, str]: - """Returns latitude band and tile name for SRTM tile from coordinates. - - Arguments: - lat (float): Latitude. - lon (float): Longitude. - - Returns: - dict: Tile parameters. - """ - lat, lon = args - - tile_latitude = math.floor(lat) - tile_longitude = math.floor(lon) - - latitude_band = f"N{abs(tile_latitude):02d}" if lat >= 0 else f"S{abs(tile_latitude):02d}" - if lon < 0: - tile_name = f"{latitude_band}W{abs(tile_longitude):03d}" - else: - tile_name = f"{latitude_band}E{abs(tile_longitude):03d}" - - self.logger.debug( - "Detected tile name: %s for coordinates: lat %s, lon %s.", tile_name, lat, lon - ) - return {"latitude_band": latitude_band, "tile_name": tile_name} -``` - -This method is required to understand how to format the download url. Of course different sources store data in different ways, so by default in the parent class this method is not implemented and you need to implement it in your provider. And if you're not using direct download, you obviously don't need this method. - -**Step 4:** implement the `get_numpy` method. +**Step 3:** implement the `download_tiles` method. ```python - def get_numpy(self) -> np.ndarray: - """Get numpy array of the tile. - - Returns: - np.ndarray: Numpy array of the tile. - """ - tile_parameters = self.get_tile_parameters(*self.coordinates) - tile_name = tile_parameters["tile_name"] - decompressed_tile_path = os.path.join(self.hgt_directory, f"{tile_name}.hgt") - - if not os.path.isfile(decompressed_tile_path): - compressed_tile_path = os.path.join(self.gz_directory, f"{tile_name}.hgt.gz") - if not self.get_or_download_tile(compressed_tile_path, **tile_parameters): - raise FileNotFoundError(f"Tile {tile_name} not found.") - - with gzip.open(compressed_tile_path, "rb") as f_in: - with open(decompressed_tile_path, "wb") as f_out: - shutil.copyfileobj(f_in, f_out) - - return self.extract_roi(decompressed_tile_path) + def download_tiles(self): + """Download SRTM tiles.""" + north, south, east, west = self.get_bbox() + + tiles = [] + # Look at each corner of the bbox in case the bbox spans across multiple tiles + for pair in [(north, east), (south, west), (south, east), (north, west)]: + tile_parameters = self.get_tile_parameters(*pair) + tile_name = tile_parameters["tile_name"] + decompressed_tile_path = os.path.join(self.hgt_directory, f"{tile_name}.hgt") + + if not os.path.isfile(decompressed_tile_path): + compressed_tile_path = os.path.join(self.gz_directory, f"{tile_name}.hgt.gz") + if not self.get_or_download_tile(compressed_tile_path, **tile_parameters): + raise FileNotFoundError(f"Tile {tile_name} not found.") + + with gzip.open(compressed_tile_path, "rb") as f_in: + with open(decompressed_tile_path, "wb") as f_out: + shutil.copyfileobj(f_in, f_out) + tiles.append(decompressed_tile_path) + + return tiles ``` -As you can see, we're using the `get_tile_parameters` method, that we've implemented earlier. Then we're downloading the tile, decompressing it and extracting the region of interest. The `get_or_download_tile` and -`extract_roi` methods are implemented in the parent class, so you don't need to reimplement them if you're using the same approach. +This method uses the helper method `get_bbox` to get the coordinates of the bounding box of the map area. If your DTM provider requires you to provide the coordinates in a different projection, you need to make sure you convert. For an example of this, see the `transform_bbox` method in [nrw.py](../maps4fs/generator/dtm/nrw.py). +Then, it determines which tiles are needed, downloads them all to a temporary folder and extracts them. The base class provides a `_tile_directory` property for convenience that points to a temp folder for your provider. +Finally, it returns a list of file paths to the downloaded tiles. As you can see, it's pretty simple to implement a DTM provider. You can use any source of elevation data, as long as it's free and open. NOTE: DTM Providers which require API keys, paid subscriptions, or any other form of payment will not be considered for implementation in the generator. @@ -156,9 +137,9 @@ mesh_z_scaling_factor = self.map.shared_settings.mesh_z_scaling_factor Here's the list of the shared settings, which directly related to the DTM Provider: -- `mesh_z_scaling_factor`: the scaling factor for the background terrain and water mesh. The simple explanation would be the following: to 3D model work properly it's Z coordinates should match the meters from real world. -Imagine the following: the highest point of your terrain is 200 meters, but in the 16-bit DEM file it's represented as 20000. So, the Z scaling factor should be 100.0. -Example of usage: +- `mesh_z_scaling_factor`: the scaling factor for the background terrain and water mesh. The simple explanation would be the following: to 3D model work properly it's Z coordinates should match the meters from real world. + Imagine the following: the highest point of your terrain is 200 meters, but in the 16-bit DEM file it's represented as 20000. So, the Z scaling factor should be 100.0. + Example of usage: ```python data: np.ndarray @@ -205,7 +186,8 @@ if some_condition: ``` ### Info sequence -If you want your provider to add some information to the `generataion_info.json` file, you can use the `data_info` property of the `DTMProvider` class. + +If you want your provider to add some information to the `generation_info.json` file, you can use the `data_info` property of the `DTMProvider` class. Note, that the `data_info` must me a correct JSON-serializable dictionary. @@ -233,4 +215,4 @@ The method in the example adds some basic information about the DEM image to the ### I implemented a DTM provider, what's next? -If you've implemented a DTM provider, you just need to create a pull request to the repository with the generator. After the review, your provider will be added to the generator and will be available for everyone to use. \ No newline at end of file +If you've implemented a DTM provider, you just need to create a pull request to the repository with the generator. After the review, your provider will be added to the generator and will be available for everyone to use.