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)
+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
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'
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
# 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