Use shapes for DTM coverage

dev/generate_extents.py

@@ -0,0 +1,80 @@
+""" Generate a shapefile with the extents of the regions for all providers.
+Source of this data: https://www.naturalearthdata.com/downloads/"""
+
+import pandas
+import shapely
+import geopandas
+
+# prepare state and country files for uniform reading
+states = geopandas.read_file("ne_10m_admin_1_states_provinces.zip")
+countries = geopandas.read_file("ne_50m_admin_0_countries.zip")
+countries.rename(columns={"NAME": "name"}, inplace=True)
+
+combined = pandas.concat([states, countries])
+
+# handling by name (for most of them)
+names = [
+    "Baden-Württemberg",
+    "Bayern",
+    "Canada",
+    "Czechia",
+    "Denmark",
+    "Finland",
+    "France",
+    "Hessen",
+    "Italy",
+    "Mecklenburg-Vorpommern",
+    "Niedersachsen",
+    "Norway",
+    "Nordrhein-Westfalen",
+    "Sachsen-Anhalt",
+    "Spain",
+    "Switzerland",
+    "United States of America",
+    "Antarctica",
+]
+names_set = combined[combined["name"].isin(names)]
+
+# handling by state/province
+england_set = combined[combined["geonunit"] == "England"]
+union = shapely.union_all(england_set["geometry"])
+england_union = geopandas.GeoDataFrame(columns=["name", "geometry"], data=[["England", union]])
+
+scotland_set = combined[combined["geonunit"] == "Scotland"]
+union = shapely.union_all(scotland_set["geometry"])
+scotland_union = geopandas.GeoDataFrame(columns=["name", "geometry"], data=[["Scotland", union]])
+
+flanders_set = combined[combined["region"] == "Flemish"]
+union = shapely.union_all(flanders_set["geometry"])
+flanders_union = geopandas.GeoDataFrame(columns=["name", "geometry"], data=[["Flanders", union]])
+
+# manual bounds that don't correspond to a specific country
+srtm = geopandas.GeoDataFrame({"name": ["SRTM"], "geometry": [shapely.box(-180, -90, 180, 90)]})
+
+# Arctic is a special case
+arctic_names = [
+    "Canada",
+    "Sweden",
+    "Finland",
+    "Iceland",
+    "Norway",
+    "Russia",
+    "United States of America",
+    "Greenland",
+]
+arctic_set = combined[combined["name"].isin(arctic_names)]
+union = shapely.union_all(arctic_set["geometry"])
+arctic_clipped = shapely.clip_by_rect(
+    union, -180, 60.7492704708152, 179.99698443265999, 83.98823036056658
+)
+arctic = geopandas.GeoDataFrame(columns=["name", "geometry"], data=[["Arctic", arctic_clipped]])
+
+# combine all
+result = pandas.concat([names_set, england_union, scotland_union, flanders_union, arctic, srtm])[
+    ["name", "geometry"]
+]
+
+result.to_file(
+    "../maps4fs/generator/dtm/extents.shz",
+    driver="ESRI Shapefile",
+)

dev/requirements.txt

@@ -23,4 +23,5 @@ tqdm
 types-tqdm
 scipy
 schedule
-streamlit-folium
+streamlit-folium
+geopandas

maps4fs/generator/dtm/arctic.py

@@ -18,7 +18,7 @@ class ArcticProvider(DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
 
-    _extents = [(83.98823036056658, 50.7492704708152, 179.99698443265999, -180)]
+    _extents_identifier = "Arctic"
 
     _instructions = (
         "This provider source includes 2 meter DEM data for the entire Arctic region above 50 "

maps4fs/generator/dtm/baden.py

@@ -16,7 +16,7 @@ class BadenWurttembergProvider(WCSProvider, DTMProvider):
     _is_community = True
     _instructions = None
     _is_base = False
-    _extents = [(49.79645444804715, 47.52877040346605, 10.54203149250156, 7.444081717803481)]
+    _extents_identifier = "Baden-Württemberg"
 
     _url = "https://owsproxy.lgl-bw.de/owsproxy/wcs/WCS_INSP_BW_Hoehe_Coverage_DGM1"
     _wcs_version = "2.0.1"

maps4fs/generator/dtm/bavaria.py

@@ -23,7 +23,7 @@ class BavariaProvider(DTMProvider):
     _author = "[H4rdB4se](https://github.com/H4rdB4se)"
     _is_community = True
     _instructions = None
-    _extents = [(50.56, 47.25, 13.91, 8.95)]
+    _extents_identifier = "Bayern"
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)

maps4fs/generator/dtm/canada.py

@@ -15,7 +15,8 @@ class CanadaProvider(WCSProvider, DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
     _is_base = False
-    _extents = [(76.49491845750764, 33.66564101989275, -26.69697497450798, -157.7322455868316)]
+    _extents_identifier = "Canada"
+
     _instructions = (
         "HRDEM coverage for Canada is limited. Make sure to check the "
         "[coverage map](https://geo.ca/imagery/high-resolution-digital"

maps4fs/generator/dtm/czech.py

@@ -16,9 +16,8 @@ class CzechProvider(WCSProvider, DTMProvider):
     _is_community = True
     _instructions = None
     _is_base = False
-    _extents = [
-        (51.0576876059846754, 48.4917065572081754, 18.9775933665038821, 12.0428143585602161)
-    ]
+
+    _extents_identifier = "Czechia"
 
     _url = "https://ags.cuzk.cz/arcgis2/services/INSPIRE_Nadmorska_vyska/ImageServer/WCSServer"  # pylint: disable=line-too-long
     _wcs_version = "1.0.0"

maps4fs/generator/dtm/denmark.py

@@ -22,7 +22,7 @@ class DenmarkProvider(WCSProvider, DTMProvider):
     _is_community = True
     _is_base = False
     _settings = DenmarkProviderSettings
-    _extents = [(57.7690657013977, 54.4354651516217, 15.5979112056959, 8.00830949937517)]
+    _extents_identifier = "Denmark"
 
     _instructions = (
         "ℹ️ This provider requires an access token. See [here](https://confluence"

maps4fs/generator/dtm/dtm.py

@@ -9,6 +9,8 @@
 from typing import TYPE_CHECKING, Type
 from zipfile import ZipFile
 
+import geopandas
+from geopandas import GeoDataFrame
 import numpy as np
 import osmnx as ox
 import rasterio
@@ -17,6 +19,7 @@
 from rasterio.enums import Resampling
 from rasterio.merge import merge
 from rasterio.warp import calculate_default_transform, reproject
+import shapely
 from tqdm import tqdm
 
 from maps4fs.logger import Logger
@@ -46,8 +49,8 @@ class DTMProvider(ABC):
     _is_base: bool = False
     _settings: Type[DTMProviderSettings] | None = DTMProviderSettings
 
-    """Bounding box of the provider in the format (north, south, east, west)."""
-    _extents: list[tuple[float, float, float, float]] | None = None
+    """Name of the entry in the shapes dataset that corresponds to the provider."""
+    _extents_identifier: str | None = None
 
     _instructions: str | None = None
 
@@ -225,35 +228,46 @@ def get_provider_by_code(cls, code: str) -> Type[DTMProvider] | None:
         return None
 
     @classmethod
-    def get_valid_provider_descriptions(cls, lat_lon: tuple[float, float]) -> dict[str, str]:
+    def get_valid_provider_descriptions(
+        cls, lat_lon: tuple[float, float], size: int
+    ) -> dict[str, str]:
         """Get descriptions of all providers, where keys are provider codes and
         values are provider descriptions.
 
         Returns:
             dict: Provider descriptions.
         """
+
+        df = geopandas.read_file(os.path.join(os.getcwd(), "maps4fs/generator/dtm/extents.shz"))
+
         providers = {}
         for provider in cls.__subclasses__():
             # pylint: disable=W0212
-            if not provider._is_base and provider.inside_bounding_box(lat_lon):
+            if not provider._is_base and provider.inside_bounding_box(df, lat_lon, size):
                 providers[provider._code] = provider.description()  # pylint: disable=W0212
         return providers  # type: ignore
 
     @classmethod
-    def inside_bounding_box(cls, lat_lon: tuple[float, float]) -> bool:
+    def inside_bounding_box(
+        cls, extents: GeoDataFrame, lat_lon: tuple[float, float], size: int
+    ) -> bool:
         """Check if the coordinates are inside the bounding box of the provider.
 
         Returns:
             bool: True if the coordinates are inside the bounding box, False otherwise.
         """
         lat, lon = lat_lon
-        extents = cls._extents
-        if extents is None:
+        bbox = ox.utils_geo.bbox_from_point((lat, lon), dist=size // 2, project_utm=False)
+        box = shapely.geometry.box(*bbox)
+
+        provider_extents = extents[extents["name"] == cls._extents_identifier]
+        if provider_extents is None:
             return True
-        for extent in extents:
-            if extent[0] >= lat >= extent[1] and extent[2] >= lon >= extent[3]:
-                return True
-        return False
+
+        geom = provider_extents["geometry"]
+        is_inside = geom.contains_properly(box).any()
+
+        return is_inside
 
     @abstractmethod
     def download_tiles(self) -> list[str]:
@@ -316,6 +330,7 @@ def get_bbox(self) -> tuple[float, float, float, float]:
             self.coordinates, dist=self.size // 2, project_utm=False
         )
         bbox = float(north), float(south), float(east), float(west)
+
         return bbox
 
     def download_tif_files(self, urls: list[str], output_path: str) -> list[str]:

maps4fs/generator/dtm/england.py

@@ -16,7 +16,7 @@ class England1MProvider(WCSProvider, DTMProvider):
     _is_community = True
     _instructions = None
     _is_base = False
-    _extents = [(55.87708724246775, 49.85060473351981, 2.0842821419111135, -7.104775741839742)]
+    _extents_identifier = "England"
 
     _url = "https://environment.data.gov.uk/geoservices/datasets/13787b9a-26a4-4775-8523-806d13af58fc/wcs"  # pylint: disable=line-too-long
     _wcs_version = "2.0.1"

maps4fs/generator/dtm/finland.py

@@ -24,7 +24,7 @@ class FinlandProvider(WCSProvider, DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
     _is_base = False
-    _extents = [(70.09, 59.45, 31.59, 19.08)]
+    _extents_identifier = "Finland"
 
     _url = "https://avoin-karttakuva.maanmittauslaitos.fi/ortokuvat-ja-korkeusmallit/wcs/v2"
     _wcs_version = "2.0.1"

maps4fs/generator/dtm/flanders.py

@@ -15,7 +15,7 @@ class FlandersProvider(WCSProvider, DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
     _is_base = False
-    _extents = [(51.5150730375579684, 50.6694992827160817, 5.9444417082210812, 2.5170092434134252)]
+    _extents_identifier = "Flanders"
 
     _url = "https://geo.api.vlaanderen.be/el-dtm/wcs"
     _wcs_version = "1.0.0"

maps4fs/generator/dtm/france.py

@@ -19,18 +19,8 @@ class FranceProvider(DTMProvider):
 
     _url = "https://data.cquest.org/ign/rgealti/repack/cog/RGEALTI_2-0_1M_COG_LAMB93-IGN69_FXX.vrt"
     _is_base = False
-    # no extents, because it also has a few colonies throughout the world
-    _extents = [
-        (51.2, 41.333, 9.55, -5.225),  # France
-        (8.6038842, 1.1710017, -61.414905, -56.4689543),  # Guyana
-        (16.5144664, 15.8320085, -61.809764, -61.0003663),  # Guadeloupe
-        (14.8787029, 14.3948596, -61.2290815, -60.8095833),  # Martinique
-        (-12.6365902, -13.0210119, 45.0183298, 45.2999917),  # Mayotte
-        (-20.8717136, -21.3897308, 55.2164268, 55.8366924),  # Reunion
-        (18.1375569, 17.670931, -63.06639, -62.5844019),  # Saint Barthelemy
-        (18.1902778, 17.8963535, -63.3605643, -62.7644063),  # Saint Martin
-        (47.365, 46.5507173, -56.6972961, -55.9033333),  # Saint Pierre and Miquelon
-    ]
+
+    _extents_identifier = "France"
 
     def download_tiles(self) -> list[str]:
         with rasterio.open(

maps4fs/generator/dtm/hessen.py

@@ -15,7 +15,7 @@ class HessenProvider(WCSProvider, DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
     _is_base = False
-    _extents = [(51.66698, 49.38533, 10.25780, 7.72773)]
+    _extents_identifier = "Hessen"
 
     _url = "https://inspire-hessen.de/raster/dgm1/ows"
     _wcs_version = "2.0.1"

maps4fs/generator/dtm/italy.py

@@ -18,7 +18,7 @@ class ItalyProvider(WCSProvider, DTMProvider):
     _is_community = True
     _instructions = None
     _is_base = False
-    _extents = [(47.15570815704503, 35.177652867276855, 19.720144130809693, 6.527697471770745)]
+    _extents_identifier = "Italy"
 
     _url = "http://tinitaly.pi.ingv.it/TINItaly_1_1/wcs"
     _wcs_version = "2.0.1"

maps4fs/generator/dtm/mv.py

@@ -27,7 +27,7 @@ class MecklenburgVorpommernProvider(WCSProvider, DTMProvider):
     _instructions = None
     _is_base = False
     _settings = MecklenburgVorpommernProviderSettings
-    _extents = [(54.8, 53, 14.5, 10.5)]
+    _extents_identifier = "Mecklenburg-Vorpommern"
 
     _url = "https://www.geodaten-mv.de/dienste/dgm_wcs"
     _wcs_version = "2.0.1"

maps4fs/generator/dtm/niedersachsen.py

@@ -20,7 +20,7 @@ class NiedersachsenProvider(WMSProvider, DTMProvider):
         "to smooth the data."
     )
     _is_base = False
-    _extents = [(54.148101, 51.153098, 11.754046, 6.505772)]
+    _extents_identifier = "Niedersachsen"
 
     _url = "https://opendata.lgln.niedersachsen.de/doorman/noauth/dgm_wms"
     _source_crs = "EPSG:25832"

maps4fs/generator/dtm/norway.py

@@ -16,9 +16,8 @@ class NorwayProvider(WCSProvider, DTMProvider):
     _is_community = True
     _instructions = None
     _is_base = False
-    _extents = [
-        (72.1016879476356962, 57.2738836442695103, 33.3365910058243742, -2.0075617181675725)
-    ]
+
+    _extents_identifier = "Norway"
 
     _instructions = (
         "This is a topobathy dataset which means it includes water depth information as well. "

maps4fs/generator/dtm/nrw.py

@@ -15,7 +15,7 @@ class NRWProvider(WCSProvider, DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
     _is_base = False
-    _extents = [(52.6008271, 50.1506045, 9.5315425, 5.8923538)]
+    _extents_identifier = "Nordrhein-Westfalen"
 
     _url = "https://www.wcs.nrw.de/geobasis/wcs_nw_dgm"
     _wcs_version = "2.0.1"

maps4fs/generator/dtm/rema.py

@@ -18,7 +18,7 @@ class REMAProvider(DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
 
-    _extents = [(-53.5443873459092, -53.5443873459092, 179.99698443265999, -180)]
+    _extents_identifier = "Antarctica"
 
     _instructions = (
         "This provider source includes 2 meter DEM data for the entire Antarctic region below 53 "

maps4fs/generator/dtm/sachsenanhalt.py

@@ -15,9 +15,8 @@ class SachsenAnhaltProvider(WCSProvider, DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
     _is_base = False
-    _extents = [
-        (53.0769416826493412, 50.8927195980075453, 13.3232545527125836, 10.5092298520646867)
-    ]
+
+    _extents_identifier = "Sachsen-Anhalt"
 
     _url = "https://www.geodatenportal.sachsen-anhalt.de/wss/service/ST_LVermGeo_DGM1_WCS_OpenData/guest"  # pylint: disable=line-too-long
     _wcs_version = "1.0.0"

maps4fs/generator/dtm/scotland.py

@@ -58,9 +58,8 @@ class ScotlandProvider(DTMProvider):
         "Coverage for Scotland is very limited. "
         "Make sure to check the [coverage map](https://remotesensingdata.gov.scot/data#/map)."
     )
-    _extents = [
-        (60.2151105070992756, 54.5525982243521881, -1.1045617513147328, -6.7070796770431951)
-    ]
+
+    _extents_identifier = "Scotland"
 
     _url = "https://srsp-catalog.jncc.gov.uk/search/product"
 

maps4fs/generator/dtm/spain.py

@@ -15,9 +15,8 @@ class SpainProvider(WCSProvider, DTMProvider):
     _author = "[kbrandwijk](https://github.com/kbrandwijk)"
     _is_community = True
     _is_base = False
-    _extents = [
-        (43.9299999999999997, 27.6299999999999990, 4.9400000000000004, -18.2100000000000009)
-    ]
+
+    _extents_identifier = "Spain"
 
     _url = "https://servicios.idee.es/wcs-inspire/mdt"
     _wcs_version = "2.0.1"