Create module wrdem and script get_elevation.
[philipp/winterrodeln/wrpylib.git] / wrpylib / wrdem.py
1 import math
2
3 import fiona.crs
4 import fiona.transform
5 import rasterio
6 import rasterio.sample
7
8 from wrpylib.wrvalidators import LonLat
9
10 LON_LAT_CRS = fiona.crs.from_epsg(4326)
11
12
13 def get_ele_from_raster(filename: str, lon_lat: LonLat, band: int = 1) -> float:
14     """
15     Examples of filenames:
16     * '/path/to/file.tif'
17     * 'file:///path/to/file.tif'
18     * 'zip:///path/to/file.zip!/folder/file.tif'
19     * 'zip+file:///path/to/file.zip!/folder/file.tif'
20     """
21     with rasterio.open(filename) as dataset:
22         xs, ys = fiona.transform.transform(LON_LAT_CRS, dataset.crs.to_proj4(), [lon_lat.lon], [lon_lat.lat])
23         xy_list = [(x, y) for x, y in zip(xs, ys)]
24         ele = list(rasterio.sample.sample_gen(dataset, xy_list, band))
25         return None if ele is None else float(ele[0])
26
27
28 class DemBase:
29     def get_ele(self, lon_lat: LonLat) -> float:
30         return math.nan
31
32     def get_name(self) -> str:
33         return self.__class__.__name__
34
35
36 class DemBasemap(DemBase):
37     def __init__(self):
38         self.filename = '/home/philipp/daten/GeoData/dem/oesterreich_10m/dhm_lamb_10m.tif'
39         # self.filename = 'zip:///home/philipp/daten/GeoData/dem/oesterreich_10m/ogd-10m-at.zip!dhm_lamb_10m.tif'
40
41     def get_ele(self, lon_lat: LonLat) -> float:
42         return get_ele_from_raster(self.filename, lon_lat)