Implement DEM for South Tyrol.

diff --git a/scripts/get_elevation.py b/scripts/get_elevation.py
index 94d3a6fefe99228d9f1b9f4c7450f0215ed4c098..c7f8cc0acde33a33c9073ca18dc5545878aa5483 100644 (file)
--- a/scripts/get_elevation.py
+++ b/scripts/get_elevation.py
@@ -1,12 +1,12 @@
 import argparse
 
-from wrpylib.wrdem import DemBasemap
+from wrpylib.wrdem import DemBasemap, DemSouthTyrol
 from wrpylib.wrvalidators import lonlat_from_str, LonLat, lonlat_to_str
 
 
 def main(lonlat: LonLat):
     """Coordinate like "47.275712 N 11.3456496 E"."""
-    dems = [DemBasemap()]
+    dems = [DemBasemap(), DemSouthTyrol()]
     for dem in dems:
         ele = dem.get_ele(lonlat)
         if ele is not None:

diff --git a/tests/test_wrdem.py b/tests/test_wrdem.py
index 8fbd1201697228c0f867ff704ec1bbbab68d2f55..8f166db49b4c2d0cec990ad045d600d2bc7885d6 100644 (file)
--- a/tests/test_wrdem.py
+++ b/tests/test_wrdem.py
@@ -1,23 +1,45 @@
 import unittest
 
-from wrpylib.wrdem import get_ele_from_raster
+import owslib.wms
+
+from wrpylib.wrdem import get_ele_from_raster, get_ele_from_wms
 from wrpylib.wrvalidators import LonLat
 
 
-class TestWrDem(unittest.TestCase):
+class TestWrDemRaster(unittest.TestCase):
     def setUp(self):
         self.filename = '/home/philipp/daten/GeoData/dem/oesterreich_10m/dhm_lamb_10m.tif'
         # https://www.data.gv.at/katalog/dataset/dgm
 
     def test_get_ele_from_raster_inside_valid(self):
-        value = get_ele_from_raster(self.filename, LonLat(13.682109, 47.934012))
+        value = get_ele_from_raster(LonLat(13.682109, 47.934012), self.filename)
         ele = round(int(value))
         self.assertEqual(ele, 557)
 
     def test_get_ele_from_raster_inside_invalid(self):
-        value = get_ele_from_raster(self.filename, LonLat(14.67656, 48.16182))
+        value = get_ele_from_raster(LonLat(14.67656, 48.16182), self.filename)
         self.assertIsNone(value)
 
     def test_get_ele_from_raster_outside(self):
-        value = get_ele_from_raster(self.filename, LonLat(8.67656, 47.16182))
+        value = get_ele_from_raster(LonLat(8.67656, 47.16182), self.filename)
         self.assertIsNone(value)
+
+
+class TestWrDemWms(unittest.TestCase):
+    def setUp(self):
+        url = 'http://geoservices.buergernetz.bz.it/geoserver/p_bz-elevation/ows?' \
+              'SERVICE=WMS&request=getCapabilities'
+        self.wms = owslib.wms.WebMapService(url, '1.3.0')
+        self.layer = 'DTM-2p5m'
+
+    def test_get_ele_from_wms_inside_valid(self):
+        ele = get_ele_from_wms(LonLat(11.31786, 46.68018), self.wms, self.layer)
+        self.assertAlmostEqual(ele, 2140.460, delta=0.1)
+
+    def test_get_ele_from_wms_inside_invalid(self):
+        ele = get_ele_from_wms(LonLat(10.893336, 46.849965), self.wms, self.layer)
+        self.assertIsNone(ele)
+
+    def test_get_ele_from_wms_outside(self):
+        ele = get_ele_from_wms(LonLat(9.908062, 45.948257), self.wms, self.layer)
+        self.assertIsNone(ele)

diff --git a/wrpylib/wrdem.py b/wrpylib/wrdem.py
index ead79e1e4c2a65b86914ec8a81bba5957a383078..b43cec89df37691014354dd4f8fc160871cddb39 100644 (file)
--- a/wrpylib/wrdem.py
+++ b/wrpylib/wrdem.py
@@ -1,8 +1,11 @@
+import json
 import math
-from typing import Optional
+from typing import Optional, Tuple
 
 import fiona.crs
 import fiona.transform
+import owslib.crs
+import owslib.wms
 import rasterio
 import rasterio.sample
 
@@ -11,7 +14,7 @@ from wrpylib.wrvalidators import LonLat
 LON_LAT_CRS = fiona.crs.from_epsg(4326)
 
 
-def get_ele_from_raster(filename: str, lon_lat: LonLat, band: int = 1) -> Optional[float]:
+def get_ele_from_raster(lon_lat: LonLat, filename: str, band: int = 1) -> Optional[float]:
     """
     Examples of filenames:
     * '/path/to/file.tif'
@@ -26,6 +29,44 @@ def get_ele_from_raster(filename: str, lon_lat: LonLat, band: int = 1) -> Option
         return None if ele == dataset.nodata else float(ele)
 
 
+def _shift_lon_lat(orig: LonLat, west_m: float, north_m: float) -> LonLat:
+    """Shifts a given geographic coordinate by a given metric amount"""
+    meter_per_degree = 10_000_000 / 90  # distance equator to pole is 10^7 m
+    return LonLat(orig.lon + (west_m / meter_per_degree) * math.cos(orig.lat * math.pi / 180.),
+                  orig.lat + north_m / meter_per_degree)
+
+
+def get_ele_from_wms(lon_lat: LonLat, wms: owslib.wms.WebMapService, layer_name: str) -> Optional[float]:
+    nan_value = -9999
+    bbox_width_px = 3
+    bbox_height_px = 3
+    bbox_width_m = 5.
+    bbox_height_m = bbox_width_m
+    left_bottom = _shift_lon_lat(lon_lat, -bbox_height_m / 2., -bbox_height_m / 2.)
+    right_top = _shift_lon_lat(lon_lat, bbox_height_m / 2., bbox_height_m / 2.)
+    response = wms.getfeatureinfo(
+        layers=[layer_name],
+        srs='CRS:84',
+        bbox=(left_bottom.lon, left_bottom.lat, right_top.lon, right_top.lat),
+        size=(bbox_width_px, bbox_height_px),
+        info_format='application/json',
+        xy=(bbox_width_px // 2, bbox_height_px // 2),
+    )
+    info_json_str = response.read().decode()
+    info_json = json.loads(info_json_str)
+    assert info_json['type'] == 'FeatureCollection'
+    features = info_json['features']
+    assert len(features) <= 1
+    if len(features) == 0:
+        return None
+    feature = features[0]
+    assert feature['type'] == 'Feature'
+    ele = feature['properties']['GRAY_INDEX']
+    if ele == nan_value:
+        return None
+    return ele
+
+
 class DemBase:
     def get_ele(self, lon_lat: LonLat) -> float:
         return math.nan
@@ -40,4 +81,36 @@ class DemBasemap(DemBase):
         # self.filename = 'zip:///home/philipp/daten/GeoData/dem/oesterreich_10m/ogd-10m-at.zip!dhm_lamb_10m.tif'
 
     def get_ele(self, lon_lat: LonLat) -> float:
-        return get_ele_from_raster(self.filename, lon_lat)
+        return get_ele_from_raster(lon_lat, self.filename)
+
+
+class DemWmsBase(DemBase):
+    def __init__(self, url: str, wms_version: str, layer_name: str):
+        self.url = url
+        self.wms_version = wms_version
+        self.layer_name = layer_name
+        self.wms = None
+
+    def get_wms(self) -> owslib.wms.WebMapService:
+        if self.wms is None:
+            self.wms = owslib.wms.WebMapService(self.url, version=self.wms_version)
+        return self.wms
+
+    def get_bbox(self) -> Tuple[float, float, float, float]:
+        """Return the bounding box in Lon/Lat WGS84"""
+        layer = self.wms.contents[self.layer_name]
+        return layer.boundingBoxWGS84  # (10.36155279782108, 46.18071678348443, 12.53362295003934, 47.132345455408874)
+
+    def get_ele(self, lon_lat: LonLat) -> float:
+        wms = self.get_wms()
+        return get_ele_from_wms(lon_lat, wms, self.layer_name)
+
+
+class DemSouthTyrol(DemWmsBase):
+    def __init__(self):
+        url = 'http://geoservices.buergernetz.bz.it/geoserver/p_bz-elevation/ows?' \
+              'SERVICE=WMS&request=getCapabilities'
+        super().__init__(url, '1.3.0', 'DTM-2p5m')
+
+    def get_bbox(self) -> Tuple[float, float, float, float]:
+        return 10.36155279782108, 46.18071678348443, 12.53362295003934, 47.132345455408874