scripts/dem_tmp.py

   1 #!/usr/bin/python3.4
   2 """Get most accurate elevation for a given point from predefined GDAL data sources or WMS sources.
   3
   4 Terms:
   5  * wrgeo: string with coordinates in Winterrodeln notation, like '47.275712 N 11.3456496 E'
   6  * lonlat: tuple with WGS84 referenced longitude, latitude pair in degrees.
   7  * xy: projected pair of coordinates
   8
   9 Useful resource
  10 https://pcjericks.github.io/py-gdalogr-cookbook/
  11
  12 Corner Coordinates:
  13     Upper Left  (   73120.000,  247300.000)
  14     Lower Left  (   73120.000,  230670.000)
  15     Upper Right (   85220.000,  247300.000)
  16     Lower Right (   85220.000,  230670.000)
  17     Center      (   79170.000,  238985.000)
  18
  19 gdaltransform -t_srs WGS84 -s_srs EPSG:31254
  20 80000 240000
  21 11.3866966549679 47.2939095182383 31.7259142374596
  22
  23
  24 Plan:
  25 (1) Get coordinates from command line in "Winterrodeln" format
  26 (2) Transform them to x/y float values
  27 (3) Transform WGS84 x/y values to EPSG:31254 x/y value
  28 (4) Get the nearest elevation from the DEM map
  29
  30 See: file:///usr/share/doc/libgdal-doc/ogr/osr_tutorial.html
  31
  32
  33 >>> import dem
  34 >>> print(dem.wrgeo_to_ele_tirol('47.275712 N 11.3456496 E'))
  35 844
  36 >>> print(dem.wrgeo_to_ele_vbg('47.435870 N 10.085163 E'))
  37 1426.52
  38
  39 """
  40 import re
  41 import struct
  42 import urllib.request, urllib.parse, urllib.error
  43 from urllib.parse import urlencode
  44 import xml.etree.ElementTree as et
  45 import os.path
  46 import zipfile
  47 import shutil
  48 import argparse
  49 from collections import namedtuple
  50 import osgeo.osr
  51 import osgeo.gdal
  52 import osgeo.ogr
  53 from osgeo.gdalconst import *
  54
  55 osgeo.gdal.UseExceptions()
  56
  57
  58 BBox = namedtuple('BBox', ['minx', 'maxx', 'miny', 'maxy'])  # bounding box. same order as returned by .getEnvelope() of osgeo.ogr geometries
  59
  60
  61 # common functions
  62 # ----------------
  63
  64 def wrgeo_to_lonlat(wrgeo):
  65     """Transfers coordinates in the format '47.175712 N 11.1456496 E' to the tuple (11.1456496, 47.175712)."""
  66     match = re.match('(\d+\.\d+) N (\d+\.\d+) E', wrgeo)
  67     y, x = match.groups()
  68     return float(x), float(y)
  69
  70
  71 def lonlat_to_epsg(lon, lat, epsg):
  72     """lon/lat in WGS84 to specific epsg system like 31254."""
  73     # source
  74     srs = osgeo.osr.SpatialReference()
  75     srs.SetWellKnownGeogCS('WGS84')
  76     # dest
  77     srd = osgeo.osr.SpatialReference()
  78     srd.ImportFromEPSG(epsg)
  79     # transform
  80     ct = osgeo.osr.CoordinateTransformation(srs, srd)
  81     x, y, z = ct.TransformPoint(lon, lat)
  82     return x, y
  83
  84
  85 def lonlat_to_epsg31254(lon, lat):
  86     """lon/lat in WGS84 to MGI_Austria_GK_West (EPSG:31254) conversion.
  87
  88     >>> lonlat_to_epsg31254(11.39, 47.29)
  89     (80255.74881170085, 239568.73826148175)"""
  90     return lonlat_to_epsg(lon, lat, 31254)
  91
  92
  93 def lonlat_to_utm32n(lon, lat):
  94     """lon/lat in WGS84 to 'UTM Zone 32 North' (EPSG:32632) conversion.
  95
  96     >>> lonlat_to_utm32n(11.39, 47.29)
  97     (680711.5152423368, 5240161.832015872)"""
  98     return lonlat_to_epsg(lon, lat, 32632)
  99
 100
 101 def interpolate_bilinear(fx, fy, z00, z10, z01, z11):
 102     """Interpolates the elevation between the values z00, z01, z10 and z11.
 103
 104     :param fx: fraction between "left" (0) and "right" (1) in x-direction
 105     :param fy: fraction between "bottom" (0) and "top" (1) in y-direction
 106     :param z00: elevation in the "left bottom" edge (where fx=0, fy=0)
 107     :param z10: elevation in the "right bottom" edge (where fx=1, fy=0)
 108     :param z01: elevation in the "left upper" edge (where fx=0, fy=1)
 109     :param z11: elevation in the "right upper" edge (where fx=1, fy=1)
 110     """
 111     return z00*(1-fx)*(1-fy) + z10*fx*(1-fy) + z01*(1-fx)*fy + z11*fx*fy
 112
 113
 114 def degree_to_dms(degree):
 115     d = int(degree)
 116     r = (degree - d) * 60
 117     m = int(r)
 118     s = (r - m) * 60
 119     return '{:2d}d{:02d}\'{:4.2f}"'.format(d, m, s)
 120
 121
 122 # GDAL supported raster files
 123 # ---------------------------
 124
 125 # http://www.gdal.org/gdal_datamodel.html says:
 126 # "Note that the pixel/line coordinates in the above are from (0.0,0.0) at the top left corner of the top left pixel
 127 # to (width_in_pixels, height_in_pixels) at the bottom right corner of the bottom right pixel.
 128 # The pixel/line location of the center of the top left pixel would therefore be (0.5,0.5)."
 129
 130 def gdal_info(filename):
 131     """Returns similar information as the command line program gdalinfo. This can be used to assist adding
 132     a new raster based source file to this library."""
 133     dataset = osgeo.gdal.Open(filename, GA_ReadOnly)
 134     print('raster size X: {} pixel'.format(dataset.RasterXSize))
 135     print('raster size Y: {} pixel'.format(dataset.RasterYSize))
 136
 137     s_srs = osgeo.osr.SpatialReference()
 138     s_srs.ImportFromWkt(dataset.GetProjectionRef())
 139     t_srs = osgeo.osr.SpatialReference()
 140     t_srs.ImportFromEPSG(4326)
 141     # note that gdal_info reports lon/lat coordinates with the reference system t_srs = s_srs.CloneGeogCS()
 142     # In general, this gives slightly different results than t_srs.ImportFromEPSG(4326) used here.
 143     coord_trans = osgeo.osr.CoordinateTransformation(s_srs, t_srs)
 144
 145     gt = dataset.GetGeoTransform()
 146
 147     # create multipoint containing the four edges in lon/lat
 148     multipoint = osgeo.ogr.Geometry(osgeo.ogr.wkbMultiPoint)
 149     for e in range(4):
 150         i = [0, dataset.RasterXSize][e % 2]   # horizontal pixel
 151         j = [0, dataset.RasterYSize][e // 2]  # vertical pixel
 152         x, y = osgeo.gdal.ApplyGeoTransform(gt, i, j)  # map units (probably meter or something like that)
 153         lon, lat, ele = coord_trans.TransformPoint(x, y)
 154
 155         point = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
 156         point.AddPoint(lon, lat)
 157         multipoint.AddGeometry(point)
 158     env = multipoint.GetEnvelope()
 159     print('bounding box: east {}°, west {}°, south {}°, north {}°'.format(*env))
 160
 161
 162 def gdal_xy_to_ele(x, y, filename, band=1):
 163     """Returns the (linearly interpolated) elevation for the given point in the datasource specified by its filename.
 164
 165     :param x: horizontal value in the coordinate system specified by the datasource
 166     :param y: vertical value in the coordinate system specified by the datasource
 167     :param filename: filename of datasource
 168     :param band: band of the datasource containing the elevation data
 169     :return: Elevation as float or None
 170
 171     Example:
 172     >>> gdal_xy_to_ele(76665, 241185, 'ibk_10m.asc')
 173     2253.0
 174     """
 175     dataset = osgeo.gdal.Open(filename, GA_ReadOnly)
 176     assert dataset is not None
 177     assert band <= dataset.RasterCount
 178     band = dataset.GetRasterBand(band)
 179     gt = dataset.GetGeoTransform() # gt == (xleft, xscale, xrot, ytop, yscale, yrot)
 180     igt = osgeo.gdal.InvGeoTransform(gt)[1]
 181     pixel, line = osgeo.gdal.ApplyGeoTransform(igt, x, y) # pixel and line are float values
 182     if pixel < 0 or pixel >= dataset.RasterXSize or line < 0 or line >= dataset.RasterYSize:
 183         return None # out of bounding box
 184
 185     # nearest pixel (spixel, sline) left & above of the position specified by (pixel, line)
 186     # the 0.5 in the formula is because the center of pixels is in the middle of the pixels.
 187     spixel, sline = int(pixel-0.5), int(line-0.5)
 188     ele_str = band.ReadRaster(spixel, sline, 2, 2)
 189     raster_format = {GDT_Int32: '=4i', GDT_Float32: '=4f'}[band.DataType]
 190     ele_edges = struct.unpack(raster_format, ele_str)  # [ele[spixel, sline], ele[spixel+1, sline], ele[spixel, sline+1], ele[spixel+1, sline+2]]
 191     if band.GetNoDataValue() in ele_edges:
 192         return None
 193     return interpolate_bilinear(pixel-0.5-spixel, line-0.5-sline, *ele_edges)
 194
 195
 196 # WMS
 197 # ---
 198
 199 # The WMS 1.3 specification says about bounding boxes:
 200 # The relation of the Bounding Box to the map pixel matrix is that the bounding box goes around the “outside”
 201 # of the pixels of the map rather than through the centres of the map’s border pixels.
 202 # In this context, individual pixels represent an area on the ground.
 203
 204 def wms_get_capabilities(url):
 205     # create URL
 206     # url = 'http://vogis.cnv.at/mapserver/mapserv?map=i_hoehen_und_gelaende_r_wms.map&'
 207     us = urllib.parse.urlsplit(url)
 208     query = urllib.parse.parse_qs(us.query)
 209     query.update({'VERSION': '1.0.3', 'REQUEST': 'GetCapabilities', 'SERVICE': 'WMS'})
 210     gc_url = urllib.parse.SplitResult(us.scheme, us.netloc, us.path, urllib.parse.urlencode(query, doseq=True), us.fragment).geturl()
 211
 212     # request
 213     f = urllib.request.urlopen(gc_url)
 214     xml = f.read()
 215     f.close()
 216
 217     # parse response XML
 218     wms_cap = et.fromstring(xml)
 219     return wms_cap
 220
 221
 222 def wms_xy_to_qbbox(x, y, bbox, width, height):
 223     """Returns the bounding box of the 4 pixel surrounding the given coordinates."""
 224     x_idx = width * (x - bbox.minx) / (bbox.maxx - bbox.minx) - 0.5  # index of pixel of total map image in x direction (float)
 225     y_idx = height * (y - bbox.miny) / (bbox.maxy - bbox.miny) - 0.5  # index of pixel of total map image in y direction (float)
 226     x0_idx = int(x_idx)
 227     x1_idx = x0_idx + 2
 228     y0_idx = int(y_idx)
 229     y1_idx = y0_idx + 2
 230     qbbox = BBox(minx=x0_idx * (bbox.maxx - bbox.minx) / width + bbox.minx,
 231                  maxx=x1_idx * (bbox.maxx - bbox.minx) / width + bbox.minx,
 232                  miny=y0_idx * (bbox.maxy - bbox.miny) / height + bbox.miny,
 233                  maxy=y1_idx * (bbox.maxy - bbox.miny) / height + bbox.miny)
 234     return qbbox
 235
 236
 237 def wms_ij_to_ele(url, layer, crs, bbox, width, height, i, j):
 238     """Returns the elevation via a getfeatureinfo request of the pixels specified by i and j within the range
 239     specified by bbox, width and height."""
 240     us = urllib.parse.urlsplit(url)
 241     query = urllib.parse.parse_qs(us.query)
 242     query.update({'VERSION': '1.3.0', 'REQUEST': 'GetFeatureInfo', 'LAYERS': layer, 'STYLES': '',
 243                   'CRS': crs, 'BBOX': ','.join(map(str, [bbox.minx, bbox.miny, bbox.maxx, bbox.maxy])),
 244                   'WIDTH': str(width), 'HEIGHT': str(height), 'FORMAT': 'image/png',
 245                   'QUERY_LAYERS': layer, 'INFO_FORMAT': 'application/vnd.ogc.gml', 'I': str(i), 'J': str(j)})
 246     fi_url = urllib.parse.SplitResult(us.scheme, us.netloc, us.path, urllib.parse.urlencode(query, doseq=True), us.fragment).geturl()
 247
 248     # request
 249     f = urllib.request.urlopen(fi_url)
 250     xml = f.read()
 251     f.close()
 252
 253     # parse XML
 254     root = et.fromstring(xml)
 255     return float(root.find('.//value_0').text)
 256
 257
 258 def wms_xy_to_ele(x, y, url, layer, crs, bbox, width, height):
 259     """Returns the (linearly interpolated) elevation for the given point from the WMS.
 260
 261     :param x: horizontal value in the native coordinate system of used layer of the WMS service
 262     :param y: vertical value in the native coordinate system of used layer of the WMS service
 263     :param url: WMS URL
 264     :param layer: WMS layer name
 265     :param crs: native coordinate reference system, e.g. 'EPSG:31254'
 266     :param bbox: layer bounding box in native coordinate system of used layer of the WMS service (instance of BBox)
 267     :param width: number of "pixel" in bbox in x direction
 268     :param height: number of "pixel" in bbox in y direction
 269     :return: Elevation as float or None
 270
 271     Example:
 272     >>> wms_xy_to_ele(80000.0, 240000.0, url)
 273     1129.
 274     """
 275     # crs=EPSG:31254&dpiMode=7&format=image/png&layers=hoehenmodell_2011_gelaende_50cm&styles=&url=http://vogis.cnv.at/mapserver/mapserv?map%3Di_hoehen_und_gelaende_r_wms.map%26version%3D1.3.0%26
 276
 277     # find qbbox (query bounding box) containing 4 pixels surrounding x/y and centers of 4 surrounding pixels.
 278     # <BoundingBox CRS="EPSG:31254" minx="180000" miny="-75000" maxx="280000" maxy="0"/>
 279     qbbox = wms_xy_to_qbbox(x, y, bbox, width, height)
 280     ele_list = []
 281     for i in range(2):
 282         for j in range(1, -1, -1):
 283             ele = wms_ij_to_ele(
 284                 'http://vogis.cnv.at/mapserver/mapserv?map=i_hoehen_und_gelaende_r_wms.map',
 285                 'hoehenmodell_2011_gelaende_50cm',
 286                 crs, qbbox, 2, 2, i, j)
 287             if ele is None:
 288                 return None
 289             ele_list.append(ele)
 290     cbbox = BBox(
 291             minx = qbbox.minx + (qbbox.maxx-qbbox.minx)/4,
 292             maxx = qbbox.minx + (qbbox.maxx-qbbox.minx)*3/4,
 293             miny = qbbox.miny + (qbbox.maxy-qbbox.miny)/4,
 294             maxy = qbbox.miny + (qbbox.maxy-qbbox.miny)*3/4)  # this bbox has the edges at the pixel centers
 295     return interpolate_bilinear((x-cbbox.minx)/(cbbox.maxx-cbbox.minx), (y-cbbox.miny)/(cbbox.maxy-cbbox.miny), *ele_list)
 296
 297
 298 # Map sources
 299 # -----------
 300
 301 class MapSource:
 302     region = None  # Human readable string shortly describing the covered region
 303     source = None  # Short description of data provider (string)
 304     attribution = None  # attribution string
 305     resolution_m = None  # (mean) distance between two samples in meter (float)
 306     bbox = None  # bounding box of data expressed in longitude/latitude (BBox)
 307     epsg = None  # native EPSG number, e.g. 31287 for 'EPSG:31287'
 308
 309
 310 class MapSourceFile(MapSource):
 311     filename = None  # filename including path of map data
 312
 313
 314 class MapSourceWmf(MapSource):
 315     url = None  # url of WMS service
 316     layer = None  # WMS layer Name
 317
 318
 319 # Austria
 320 # -------
 321
 322 class MapSourceAustriaBasemap(MapSourceFile):
 323     region = 'Austria'
 324     source = 'basemap.at'
 325     resolution_m = 10.
 326     bbox = BBox(9.314479, 17.358544, 46.297834, 49.022520)
 327     filename = '/daten/GeoData/dem/oesterreich_10m/dhm_lamb_10m.tif'
 328     epsg = 31287
 329
 330
 331 def lonlat_to_ele_austria(lon, lat):
 332     """DEM Austria 10x10m EPSG:31287
 333     C-BY-3.0: Land Kaernten - data.ktn.gv.at"""
 334     filename = '/daten/GeoData/dem/oesterreich_10m/dhm_lamb_10m.tif'
 335     x, y = lonlat_to_epsg(lon, lat, 31287)
 336     return gdal_xy_to_ele(x, y, filename)
 337
 338
 339 # Vorarlberg
 340 # ----------
 341
 342 class MapSourceWmfVorarlberg(MapSourceWmf):
 343     region = 'Vorarlberg'
 344     source = 'Landesamt fuer Vermessung und Geoinformation Vorarlberg'
 345     resolution_m = 0.5
 346     bbox = BBox(9.33463, 10.333, 46.7549, 47.6586)
 347     url = 'http://vogis.cnv.at/mapserver/mapserv?map=i_hoehen_und_gelaende_r_wms.map&'
 348     layer = 'hoehenmodell_2011_oberflaeche_50cm'
 349     epsg = 31254
 350
 351
 352 def wms_getfeatureinfo_vbg(x, y):
 353     """Uses a WMS GetFeatureInfo call to retrieve the DEM data.
 354
 355     >>> print wms_getfeatureinfo_vbg(9.783502, 47.335013)
 356     """
 357     # Note: a better (or at least alternative) way to implement this would be to use GDAL to do it for us.
 358     #
 359     # http://www.gdal.org/frmt_wms.html
 360     # "Starting with GDAL >= 2.0, WMS layers can be queried (through a GetFeatureInfo request) with the
 361     # gdallocationinfo utility, or with a GetMetadataItem("Pixel_iCol_iLine", "LocationInfo") call on a band object.
 362     #
 363     # # This would mean:
 364     # filename='dem_vogis.xml'
 365     # dataset = osgeo.gdal.Open(filename, GA_ReadOnly)
 366     # assert dataset is not None
 367     # assert band <= dataset.RasterCount
 368     # band = dataset.GetRasterBand(1)
 369     # band.GetMetadataItem("GeoPixel_{}_{}".format(1000, 1000), "LocationInfo")
 370     # # No query is executed and None is returned.
 371     #
 372     # Using gdallocationinfo just returns the RGB values 255.
 373
 374     # The third alternative also just returns the RGB value:
 375     # print dem.gdal_xy_to_ele(-100, 200000, 'dem_vogis.xml', 1)
 376     # 255
 377
 378     delta = 0.000001
 379     xl = str(x-delta)
 380     xr = str(x+delta)
 381     yb = str(y-delta)
 382     yt = str(y+delta)
 383     parameter = [
 384         ('map', 'i_hoehen_und_gelaende_r_wms.map'),
 385         ('SERVICE', 'WMS'),
 386         ('VERSION', '1.1.1'),
 387         ('REQUEST', 'GetFeatureInfo'),
 388         ('BBOX', ','.join([xl, yb, xr, yt])),
 389         ('SRS', 'EPSG:4326'),
 390         ('WIDTH', '3'),
 391         ('HEIGHT', '3'),
 392         ('LAYERS', 'hoehenmodell_terrain_1m'),
 393         ('QUERY_LAYERS', 'hoehenmodell_terrain_1m'),
 394         ('INFO_FORMAT', 'gml'),
 395         ('X', '1'),
 396         ('Y', '1')]
 397     url = 'http://vogis.cnv.at/mapserver/mapserv?' + urllib.parse.urlencode(parameter)
 398     f = urllib.request.urlopen(url)
 399     xml = f.read()
 400     f.close()
 401     root = et.fromstring(xml)
 402     elevation = float(root.find('hoehenmodell_terrain_1m_layer/hoehenmodell_terrain_1m_feature/value_0').text)
 403     return elevation
 404
 405
 406 def wrgeo_to_ele_vbg(wrgeo):
 407     """Example:
 408     >>> wrgeo_to_ele_vbg('47.435870 N 10.085163 E')
 409     1426.52
 410     """
 411     x, y = wrgeo_to_lonlat(wrgeo)
 412     return wms_getfeatureinfo_vbg(x, y)
 413
 414
 415 # Suedtirol
 416 # ---------
 417
 418 # http://geoservices.buergernetz.bz.it/geoserver/p_bz-elevation/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&BBOX=2922671.81412236206233501,1209060.08528055297210813,2926426.12536783283576369,1212054.986808244837448&CRS=EPSG:2056&WIDTH=1142&HEIGHT=911&LAYERS=DTM-2p5m&STYLES=&FORMAT=image/png&QUERY_LAYERS=DTM-2p5m&INFO_FORMAT=application/json&I=79&J=30
 419
 420 def utm32n_to_mapcode(utm_x, utm_y):
 421     """The mapcode is the filename of the map part without extension.
 422
 423     >>> utm32n_to_mapcode(705033, 5175081)
 424     'dtm_utm10161'
 425     """
 426     # East
 427     # 604998.750 # left  edge of 09013
 428     # 608201.250 # right edge of 09013
 429     # 608198.750 # left  edge of 09012
 430     # 611401.250 # right edge of 09012
 431     # 611398.750 # left  edge of 09023
 432     # 614601.250 # right edge of 09023
 433     east_i = int((utm_x-605000) // 3200 + 2)
 434     east_major = east_i // 2
 435     east_minor = east_i %  2
 436
 437     # North
 438     # 5128401.250 # upper edge of 02113
 439     # 5125598.750 # lower edge of 02113
 440     # 5125601.250 # upper edge of 01114
 441     # 5122798.750 # lower edge of 01114
 442     # 5122801.250 # upper edge of 01113
 443     # 5119998.750 # lower edge of 01113
 444     north_i = int((utm_y-5120000) // 2800 + 2)
 445     north_major = north_i // 2
 446     north_minor = north_i %  2
 447
 448     minor = {(1, 1): 1, (1, 0): 2, (0, 0): 3, (0, 1): 4}[(east_minor, north_minor)]
 449
 450     return 'dtm_utm{:02d}{:02d}{:d}'.format(north_major, east_major, minor)
 451
 452
 453 def wrgeo_to_ele_stirol(wrgeo):
 454     lon, lat = wrgeo_to_lonlat(wrgeo)
 455     utm_x, utm_y = lonlat_to_utm32n(lon, lat)
 456     mapcode = utm32n_to_mapcode(utm_x, utm_y)
 457
 458     # Unpack zip archive in /tmp
 459     if not os.path.exists(os.path.join('/tmp', mapcode + '.ascii')):
 460         shutil.copy(os.path.join('/daten/GeoData/dem/suedtirol', mapcode + '.zip'), '/tmp')
 461         zip = zipfile.ZipFile(os.path.join('/tmp', mapcode + '.zip'), 'r')
 462         zip.extractall('/tmp')
 463
 464     # Get elevation
 465     return gdal_xy_to_ele(utm_x, utm_y, os.path.join('/tmp', mapcode + '.ascii'))
 466
 467
 468 # main program
 469 # ------------
 470
 471 def main(wrgeo, no_wms):
 472     lon, lat = wrgeo_to_lonlat(wrgeo)
 473
 474     ele = lonlat_to_ele_austria(lon, lat)
 475
 476     ele_str = str(ele) if ele is None else '{:.0f} m'.format(ele)
 477     print('lat: {:.6f}° lon: {:0.6f}° ele: {}'.format(lat, lon, ele_str))
 478
 479
 480 if __name__ == '__main__':
 481     parser = argparse.ArgumentParser('Get most accurate elevation for a given point from predefined GDAL data sources or WMS sources')
 482     parser.add_argument('--no-wms', help='don\'t use WMF services')
 483     parser.add_argument('wrgeo', help='coordinate like "47.275712 N 11.3456496 E"')
 484     args = parser.parse_args()
 485     main(args.wrgeo, args.no_wms)