Omit empty parameters in key/value templates.

[philipp/winterrodeln/wrpylib.git] / wrpylib / wrdem.py

diff --git a/wrpylib/wrdem.py b/wrpylib/wrdem.py
index 4cc9a937aa311795f9f947c4bd0db9f268f8fbe9..ce0163fbf4a9420d6382110775f62b247ce24cb2 100644 (file)
--- a/wrpylib/wrdem.py
+++ b/wrpylib/wrdem.py
@@ -3,7 +3,7 @@ import math
 from abc import abstractmethod
 from typing import Optional, Tuple
 
-import fiona.crs
+import fiona.crs  # pip install Fiona
 import fiona.transform
 import owslib.crs
 import owslib.wms
@@ -16,12 +16,10 @@ from wrpylib.wrvalidators import LonLat
 LON_LAT_CRS = fiona.crs.from_epsg(4326)
 
 
-def get_bbox_from_raster(dataset: DatasetReader) -> Tuple[Tuple[float, float, float, float], str]:
-    """Returns the bounding box in dataset CRS."""
-    return dataset.bounds, dataset.crs.to_proj4()
+XY = Tuple[float, float]
 
 
-def get_ele_from_raster(lon_lat: LonLat, dataset: DatasetReader, band: int = 1) -> Optional[float]:
+def get_ele_from_raster(dataset: DatasetReader, x_y: XY, band: int = 1, crs: str = None) -> Optional[float]:
     """
     Examples of filenames:
     * '/path/to/file.tif'
@@ -30,11 +28,13 @@ def get_ele_from_raster(lon_lat: LonLat, dataset: DatasetReader, band: int = 1)
     * 'zip+file:///path/to/file.zip!/folder/file.tif'
 
     Usage:
+    >>> lon_lat = LonLat(11.2, 47.8)
     >>> filename = '/path/to/file.tif'
     >>> with rasterio.open(filename) as dataset:
-    >>>     ele = get_ele_from_raster(lon_lat, dataset)
+    >>>     crs = dataset.crs.to_proj4()
+    >>>     x_y = transform_lon_lat(lon_lat, crs)
+    >>>     ele = get_ele_from_raster(dataset, x_y)
     """
-    x_y = transform_lon_lat(lon_lat, dataset.crs.to_proj4())
     ele = list(rasterio.sample.sample_gen(dataset, [x_y], band))[0]
     return None if ele == dataset.nodata else float(ele)
 
@@ -77,7 +77,7 @@ def get_ele_from_wms(lon_lat: LonLat, wms: owslib.wms.WebMapService, layer_name:
     return ele
 
 
-def transform_lon_lat(lon_lat: LonLat, crs: str) -> Tuple[float, float]:
+def transform_lon_lat(lon_lat: LonLat, crs: str) -> XY:
     """Transform the coordinates to dataset CRS."""
     xs, ys = fiona.transform.transform(LON_LAT_CRS, crs, [lon_lat.lon], [lon_lat.lat])
     return [(x, y) for x, y in zip(xs, ys)][0]
@@ -112,20 +112,22 @@ class DemBase:
 
 
 class DemRasterBase(DemBase):
-    def __init__(self, filename: str, band: int = 1):
+    def __init__(self, filename: str, band: int = 1, crs: str = None):
         self.filename = filename
         self.band = band
         self.dataset = rasterio.open(filename)
+        self.crs = crs if crs else self.dataset.crs.to_proj4()
 
     def close(self):
         self.dataset.close()
 
     def get_ele(self, lon_lat: LonLat) -> Optional[float]:
-        return get_ele_from_raster(lon_lat, self.dataset, self.band)
+        x_y = transform_lon_lat(lon_lat, self.crs)
+        return get_ele_from_raster(self.dataset, x_y, self.band, self.crs)
 
     def is_in_bbox(self, lon_lat: LonLat) -> bool:
-        bbox, crs = get_bbox_from_raster(self.dataset)
-        x_y = transform_lon_lat(lon_lat, crs)
+        bbox = self.dataset.bounds
+        x_y = transform_lon_lat(lon_lat, self.crs)
         return _is_in_bbox(x_y, bbox)
 
 
@@ -135,6 +137,18 @@ class DemBasemap(DemRasterBase):
         super().__init__(filename)
 
 
+class DemSwitzerland(DemRasterBase):
+    def __init__(self):
+        filename = 'zip:///home/philipp/daten/GeoData/dem/schweiz_200m/dem_schweiz_200m.zip!data/DHM200.asc'
+        super().__init__(filename, crs='EPSG:21781')
+
+
+class DemBavaria(DemRasterBase):
+    def __init__(self):
+        filename = 'zip:///home/philipp/daten/GeoData/dem/bayern_50m/dgm50_epsg25832.tif.zip!dgm50_epsg25832.tif'
+        super().__init__(filename)
+
+
 class DemWmsBase(DemBase):
     def __init__(self, url: str, wms_version: str, layer_name: str):
         self.url = url
@@ -169,7 +183,12 @@ class DemSouthTyrol(DemWmsBase):
 
 class MultiDem:
     def __init__(self):
-        self.dems = [DemBasemap(), DemSouthTyrol()]
+        self.dems = [
+            DemBasemap(),
+            DemSouthTyrol(),
+            # DemBavaria(),  # too coarse - we should have a 25 m grid or denser - this is 50 m
+            # DemSwitzerland(), # too coarse - we should have a 25 m grid or denser - this is 200 m
+        ]
 
     def close(self):
         for dem in self.dems: