Add DEM Switzerland.

diff --git a/tests/test_wrdem.py b/tests/test_wrdem.py
index 05700d57109e907c4aedc450545055a9d57aea9c..ebdc8bcf9337f10edbc0d44770a14f748d19f5ce 100644 (file)
--- a/tests/test_wrdem.py
+++ b/tests/test_wrdem.py
@@ -3,7 +3,7 @@ import unittest
 import owslib.wms
 import rasterio
 
-from wrpylib.wrdem import get_ele_from_raster, get_ele_from_wms, MultiDem
+from wrpylib.wrdem import get_ele_from_raster, get_ele_from_wms, MultiDem, transform_lon_lat, DemSwitzerland, DemBasemap
 from wrpylib.wrvalidators import LonLat
 
 
@@ -17,17 +17,19 @@ class TestWrDemRaster(unittest.TestCase):
         self.dataset.close()
 
     def test_get_ele_from_raster_inside_valid(self):
-        value = get_ele_from_raster(LonLat(13.682109, 47.934012), self.dataset)
-        ele = round(int(value))
-        self.assertEqual(ele, 557)
+        x_y = transform_lon_lat(LonLat(13.682109, 47.934012), self.dataset.crs.to_proj4())
+        ele = get_ele_from_raster(self.dataset, x_y)
+        self.assertAlmostEqual(ele, 557, delta=0.5)
 
     def test_get_ele_from_raster_inside_invalid(self):
-        value = get_ele_from_raster(LonLat(14.67656, 48.16182), self.dataset)
-        self.assertIsNone(value)
+        x_y = transform_lon_lat(LonLat(14.67656, 48.16182), self.dataset.crs.to_proj4())
+        ele = get_ele_from_raster(self.dataset, x_y)
+        self.assertIsNone(ele)
 
     def test_get_ele_from_raster_outside(self):
-        value = get_ele_from_raster(LonLat(8.67656, 47.16182), self.dataset)
-        self.assertIsNone(value)
+        x_y = transform_lon_lat(LonLat(8.67656, 47.16182), self.dataset.crs.to_proj4())
+        ele = get_ele_from_raster(self.dataset, x_y)
+        self.assertIsNone(ele)
 
 
 class TestWrDemWms(unittest.TestCase):
@@ -50,6 +52,46 @@ class TestWrDemWms(unittest.TestCase):
         self.assertIsNone(ele)
 
 
+class TestDemBasemap(unittest.TestCase):
+    def setUp(self):
+        self.dem = DemBasemap()
+
+    def tearDown(self):
+        self.dem.close()
+
+    def test_get_ele_inside_valid(self):
+        ele = self.dem.get_ele(LonLat(13.682109, 47.934012))
+        self.assertAlmostEqual(ele, 557, delta=0.5)
+
+    def test_get_ele_inside_invalid(self):
+        ele = self.dem.get_ele(LonLat(14.67656, 48.16182))
+        self.assertIsNone(ele)
+
+    def test_get_ele_outside(self):
+        ele = self.dem.get_ele(LonLat(8.67656, 47.16182))
+        self.assertIsNone(ele)
+
+
+class TestDemSwitzerland(unittest.TestCase):
+    def setUp(self):
+        self.dem = DemSwitzerland()
+
+    def tearDown(self):
+        self.dem.close()
+
+    def test_get_ele_inside_valid(self):
+        ele = self.dem.get_ele(LonLat(8.452512, 47.064715))
+        self.assertAlmostEqual(ele, 818, delta=0.5)
+
+    def test_get_ele_inside_invalid(self):
+        ele = self.dem.get_ele(LonLat(8.319, 45.970))
+        self.assertIsNone(ele)
+
+    def test_get_ele_outside(self):
+        ele = self.dem.get_ele(LonLat(5.182008, 46.799233))
+        self.assertIsNone(ele)
+
+
 class TestWrDemMultiDem(unittest.TestCase):
     def setUp(self):
         self.dem = MultiDem()

diff --git a/wrpylib/wrdem.py b/wrpylib/wrdem.py
index 4cc9a937aa311795f9f947c4bd0db9f268f8fbe9..9184247d412beed292c5f8cdd220ec366f7650ae 100644 (file)
--- a/wrpylib/wrdem.py
+++ b/wrpylib/wrdem.py
@@ -1,6 +1,7 @@
 import json
 import math
 from abc import abstractmethod
+from collections import namedtuple
 from typing import Optional, Tuple
 
 import fiona.crs
@@ -16,12 +17,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 +29,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 +78,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 +113,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 +138,14 @@ class DemBasemap(DemRasterBase):
         super().__init__(filename)
 
 
+class DemSwitzerland(DemRasterBase):
+    def __init__(self):
+        filename = '/home/philipp/daten/GeoData/dem/schweiz_200m/data/DHM200.xyz'
+        filename = '/home/philipp/daten/GeoData/dem/schweiz_200m/data/DHM200.asc'
+        # https://www.swisstopo.admin.ch/de/swisstopo/kostenlose-geobasisdaten.html
+        super().__init__(filename, crs='EPSG:21781')
+
+
 class DemWmsBase(DemBase):
     def __init__(self, url: str, wms_version: str, layer_name: str):
         self.url = url
@@ -169,7 +180,7 @@ class DemSouthTyrol(DemWmsBase):
 
 class MultiDem:
     def __init__(self):
-        self.dems = [DemBasemap(), DemSouthTyrol()]
+        self.dems = [DemBasemap(), DemSouthTyrol(), DemSwitzerland()]
 
     def close(self):
         for dem in self.dems: