From 280104369f14d3cb946bd17ea61fbff0cf2c1601 Mon Sep 17 00:00:00 2001
From: Philipp Spitzer <philipp@spitzer.priv.at>
Date: Mon, 25 Jan 2021 23:06:12 +0100
Subject: [PATCH] Create module wrdem and script get_elevation.

---
 requirements.txt         |  9 +++++++++
 scripts/get_elevation.py | 22 +++++++++++++++++++++
 setup.py                 |  1 +
 wrpylib/wrdem.py         | 42 ++++++++++++++++++++++++++++++++++++++++
 4 files changed, 74 insertions(+)
 create mode 100644 requirements.txt
 create mode 100644 scripts/get_elevation.py
 create mode 100644 wrpylib/wrdem.py

diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..ccff48c
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,9 @@
+mysqlclient~=1.3.10
+SQLAlchemy~=1.2.18
+wrpylib~=0.7.0
+mwparserfromhell~=0.5.2
+Fiona~=1.8.4
+Shapely~=1.6.4
+Pillow~=5.4.1
+rasterio~=1.0.21
+setuptools~=41.0.1
\ No newline at end of file
diff --git a/scripts/get_elevation.py b/scripts/get_elevation.py
new file mode 100644
index 0000000..94d3a6f
--- /dev/null
+++ b/scripts/get_elevation.py
@@ -0,0 +1,22 @@
+import argparse
+
+from wrpylib.wrdem import DemBasemap
+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()]
+    for dem in dems:
+        ele = dem.get_ele(lonlat)
+        if ele is not None:
+            print(f'Elevation for {lonlat_to_str(lonlat)}: {ele:.0f} m (source: {dem.get_name()})')
+            return
+    print(f'No elevation data available for {lonlat_to_str(lonlat)}')
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser('Get most accurate elevation for a given point')
+    parser.add_argument('coordinate', type=lonlat_from_str, help='coordinate like "47.275712 N 11.3456496 E"')
+    args = parser.parse_args()
+    main(args.coordinate)
diff --git a/setup.py b/setup.py
index 98ad7a4..3d2c12a 100644
--- a/setup.py
+++ b/setup.py
@@ -16,6 +16,7 @@ setup(name='wrpylib',
     packages=['wrpylib'],
     scripts=[
         'scripts/createthumbnails.py',
+        'scripts/get_elevation.py',
         'scripts/update_intermaps.py',
         'scripts/updatewrinncache.py',
         'scripts/updatewrmapcache.py',
diff --git a/wrpylib/wrdem.py b/wrpylib/wrdem.py
new file mode 100644
index 0000000..8e33ab0
--- /dev/null
+++ b/wrpylib/wrdem.py
@@ -0,0 +1,42 @@
+import math
+
+import fiona.crs
+import fiona.transform
+import rasterio
+import rasterio.sample
+
+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) -> float:
+    """
+    Examples of filenames:
+    * '/path/to/file.tif'
+    * 'file:///path/to/file.tif'
+    * 'zip:///path/to/file.zip!/folder/file.tif'
+    * 'zip+file:///path/to/file.zip!/folder/file.tif'
+    """
+    with rasterio.open(filename) as dataset:
+        xs, ys = fiona.transform.transform(LON_LAT_CRS, dataset.crs.to_proj4(), [lon_lat.lon], [lon_lat.lat])
+        xy_list = [(x, y) for x, y in zip(xs, ys)]
+        ele = list(rasterio.sample.sample_gen(dataset, xy_list, band))
+        return None if ele is None else float(ele[0])
+
+
+class DemBase:
+    def get_ele(self, lon_lat: LonLat) -> float:
+        return math.nan
+
+    def get_name(self) -> str:
+        return self.__class__.__name__
+
+
+class DemBasemap(DemBase):
+    def __init__(self):
+        self.filename = '/home/philipp/daten/GeoData/dem/oesterreich_10m/dhm_lamb_10m.tif'
+        # 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)
-- 
2.39.5