Simplify terrain module to run with raster tiles

Use GDAL VRT file to abstract over multiple raster tiles.

Use point_query to sample elevation profile.

data/S_AVE_DSM.vrt

@@ -0,0 +1,21 @@
+<VRTDataset rasterXSize="7200" rasterYSize="3600">
+  <SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>
+  <GeoTransform>  2.6000000000000000e+01,  2.7777777777777778e-04,  0.0000000000000000e+00, -3.0000000000000000e+00,  0.0000000000000000e+00, -2.7777777777777778e-04</GeoTransform>
+  <VRTRasterBand dataType="Int16" band="1">
+    <ColorInterp>Gray</ColorInterp>
+    <SimpleSource>
+      <SourceFilename relativeToVRT="1">S004E026_AVE_DSM.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SourceProperties RasterXSize="3600" RasterYSize="3600" DataType="Int16" BlockXSize="3600" BlockYSize="1" />
+      <SrcRect xOff="0" yOff="0" xSize="3600" ySize="3600" />
+      <DstRect xOff="0" yOff="0" xSize="3600" ySize="3600" />
+    </SimpleSource>
+    <SimpleSource>
+      <SourceFilename relativeToVRT="1">S004E027_AVE_DSM.tif</SourceFilename>
+      <SourceBand>1</SourceBand>
+      <SourceProperties RasterXSize="3600" RasterYSize="3600" DataType="Int16" BlockXSize="3600" BlockYSize="1" />
+      <SrcRect xOff="0" yOff="0" xSize="3600" ySize="3600" />
+      <DstRect xOff="3600" yOff="0" xSize="3600" ySize="3600" />
+    </SimpleSource>
+  </VRTRasterBand>
+</VRTDataset>

scripts/area.py

@@ -202,28 +202,22 @@ def csv_writer(data, directory, filename):
     dem_path = BASE_PATH
     directory = DATA_PROCESSED
     directory_shapes = os.path.join(DATA_PROCESSED, 'shapes')
-
-    old_crs = 'EPSG:4326'
-    new_crs = 'EPSG:3857'
     cell_range = 20000
 
-    #create new geojson for Crystal Palace radio transmitter
-    transmitter = {
-        'type': 'Feature',
-        'geometry': {
-            'type': 'Point',
-            'coordinates': (-0.07491679518573545, 51.42413477117786)
-            },
-        'properties': {
-            'id': 'Crystal Palace radio transmitter'
-            }
-        }
-
     #Terrain Irregularity Parameter delta h (in meters)
-    tip = terrain_area(dem_path, transmitter, cell_range, old_crs)
+    tip = terrain_area(
+        os.path.join(dem_path, 'ASTGTM2_N51W001_dem.tif'),
+        -0.074916,
+        51.424134,
+        cell_range)
+    print('TIP for AST DEM', tip)
 
     print('Running itmlogic')
     output = itmlogic_area(tip)
 
-    print('Writing results to .csv')
+    print('Writing results to ', os.path.join(directory, 'uarea_output.csv'))
     csv_writer(output, directory, 'uarea_output.csv')
+
+    tip = terrain_area(os.path.join(dem_path, 'S_AVE_DSM.vrt'), 26.9976, -3.5409, cell_range)
+    print('TIP for ALOS DEM', tip)
+    output = itmlogic_area(tip)

scripts/p2p.py

@@ -336,10 +336,9 @@ def straight_line_from_points(a, b):
     #create new geojson for terrain path
     line = straight_line_from_points(transmitter, receiver)
 
-    current_crs = 'EPSG:4326'
-
     #run terrain module
-    measured_terrain_profile, distance_km, points = terrain_p2p(dem_folder, line, current_crs)
+    measured_terrain_profile, distance_km, points = terrain_p2p(
+        os.path.join(dem_folder, 'ASTGTM2_N51W001_dem.tif'), line)
     print('Distance is {}km'.format(distance_km))
 
     #check (out of interest) how many measurements are in each profile
@@ -397,10 +396,9 @@ def straight_line_from_points(a, b):
     #create new geojson for terrain path
     line = straight_line_from_points(transmitter, receiver)
 
-    current_crs = 'EPSG:4326'
-
     #run terrain module
-    measured_terrain_profile, distance_km, points = terrain_p2p(dem_folder, line, current_crs)
+    measured_terrain_profile, distance_km, points = terrain_p2p(
+        os.path.join(dem_folder, 'S_AVE_DSM.vrt'), line)
 
     print("Profile [", measured_terrain_profile[0], ",  ... ,", measured_terrain_profile[-1], "]")
     print("Distance is {}km".format(distance_km))

scripts/terrain_module.py

@@ -6,169 +6,146 @@
 June 2019
 
 """
-import configparser
 import glob
-
-import os
 import math
+import os
+import sys
+from collections import OrderedDict
 from functools import partial, lru_cache
 
 import fiona
 import pyproj
 import rasterio
 import numpy as np
+from fiona.crs import from_epsg
+from rasterstats import gen_zonal_stats, point_query
 from shapely.ops import transform
 from shapely.geometry import Point, LineString, mapping
-from rasterstats import zonal_stats
-import configparser
-from collections import OrderedDict
 
-CONFIG = configparser.ConfigParser()
-CONFIG.read(os.path.join(os.path.dirname(__file__), 'script_config.ini'))
-BASE_PATH = CONFIG['file_locations']['base_path']
 
-
-def terrain_area(dem_folder, transmitter, cell_range, current_crs):
+def terrain_area(dem, lon, lat, cell_range):
     """
     This module takes a single set of point coordinates for a site
     along with an estimate of the cell range. The irregular terrain
     parameter is returned.
 
     Parameters
     ----------
-    dem_folder : string
-        Folder path to the available Digital Elevation Model tiles.
-    transmitter : dict
-        Geojson. Must be in WGS84 / EPSG: 4326
+    dem : str
+        Path to the available Digital Elevation Model as single raster file or vrt.
+    lon : float
+        Longitude of cell point
+    lat : float
+        Latitude of cell point
     cell_range : int
         Radius of cell area in meters.
-    current_crs : string
-        The coordinate reference system the transmitter coordinates
-        are in (must be WGS84 / EPSG: 4326).
 
     Returns
     -------
     Inter-decile range : int
         The terrain irregularity parameter.
 
     """
-    point_geometry = Point(transmitter['geometry']['coordinates'])
+    # Buffer around cell point
+    cell_area = geodesic_point_buffer(lon, lat, cell_range)
 
-    to_projected = pyproj.Transformer.from_proj(
-        pyproj.Proj('epsg:4326'), # source coordinate system
-        pyproj.Proj('epsg:3857')) # destination coordinate system
+    # Calculate raster stats
+    stats = next(gen_zonal_stats(
+        [cell_area],
+        dem,
+        add_stats={
+            'interdecile_range': interdecile_range
+        },
+        nodata=-9999
+    ))
 
-    point_geometry = transform(to_projected.transform, point_geometry)
+    id_range = stats['interdecile_range']
 
-    cell_area_projected = point_geometry.buffer(cell_range)
+    return id_range
 
-    to_unprojected = pyproj.Transformer.from_proj(
-        pyproj.Proj('epsg:3857'), # source coordinate system
-        pyproj.Proj('epsg:4326')) # destination coordinate system
 
-    cell_area_unprojected = transform(to_unprojected.transform, cell_area_projected)
+def geodesic_point_buffer(lon, lat, distance_m):
+    """Calculate a buffer a specified number of metres around a lat/lon point
+    """
+    # Azimuthal equidistant projection around lat/lon point
+    aeqd = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0'
+    crs = aeqd.format(lat=lat, lon=lon)
 
-    stats = zonal_stats([cell_area_unprojected],
-                        os.path.join(BASE_PATH,'ASTGTM2_N51W001_dem.tif'),
-                        add_stats={'interdecile_range':interdecile_range})
+    # To be transformed to WGS84 / EPSG:4326
+    transformer = pyproj.Transformer.from_crs(crs, "epsg:4326", always_xy=True)
 
-    return stats[0]['interdecile_range']
+    # Buffer origin by distance in metres
+    buf = Point(0, 0).buffer(distance_m)
 
+    # Return transformed
+    return transform(transformer.transform, buf)
 
-def interdecile_range(x):
-    """
-    Get range between bottom 10% and top 10% of values.
 
+def interdecile_range(x):
+    """Get range between bottom 10% and top 10% of values.
     """
     q90, q10 = np.percentile(x, [90, 10])
+    return round(q90 - q10)
 
-    return int(round(q90 - q10, 0))
+
+def all_data(x):
+    """Get all data points within mask, with values greater than zero
+    """
+    data = x.compressed()
+    return data[data > 0]
 
 
-def terrain_p2p(dem_folder, line, current_crs):
+def terrain_p2p(dem, line):
     """
     This module takes a set of point coordinates and returns
     the elevation profile.
 
-    Line : dict
+    dem : str
+        Path to the available Digital Elevation Model as single raster file or vrt.
+    line : dict
         Geojson. Must be in WGS84 / EPSG: 4326
-
     """
-    extents = load_extents(dem_folder)
     line_geom = LineString(line['geometry']['coordinates'])
 
     # Geographic distance
     geod = pyproj.Geod(ellps="WGS84")
     distance_m = geod.geometry_length(line_geom)
     distance_km = distance_m / 1e3
 
+    # Interpolate along line to get sampling points
     num_samples = determine_num_samples(distance_m)
-
     steps = np.interp(range(num_samples), [0,num_samples], [0, line_geom.length])
+    point_geoms = [line_geom.interpolate(currentdistance) for currentdistance in steps]
 
-    elevation_profile = []
-    points = []
-
-    for currentdistance in steps:
-        point = line_geom.interpolate(currentdistance)
-        xp, yp = point.x, point.y
-        tile_path = get_tile_path_for_point(extents, xp, yp)
-        z = get_value_from_dem_tile(tile_path, xp, yp)
-
-        elevation_profile.append(z)
+    # Sample elevation profile
+    elevation_profile = point_query(point_geoms, dem)
 
-        points.append({
+    # Put together point features with raster values
+    points = [
+        {
             'type': 'Feature',
             'geometry': mapping(point),
             'properties': {
                 'elevation': float(z),
             }
-        })
-
+        }
+        for point, z in zip(point_geoms, elevation_profile)
+    ]
 
     return elevation_profile, distance_km, points
 
 
-def load_extents(dem_folder):
-    """
-    Check the extent of each DEM tile, save to dict for future reference.
-
-    """
-    extents = {}
-    for tile_path in glob.glob(os.path.join(dem_folder, "*.tif")):
-        dataset = rasterio.open(tile_path)
-        extents[tuple(dataset.bounds)] = tile_path
-
-    return extents
-
-
-def get_tile_path_for_point(extents, x, y):
-    for (left, bottom, right, top), path in extents.items():
-        if x >= left and x <= right and y <= top and y >= bottom:
-            return path
-    raise ValueError("No tile includes x {} y {}".format(x, y))
-
-
-def get_value_from_dem_tile(tile_path, x, y):
-    """
-    Read all tile extents, load value from relevant tile.
-
-    """
-    dataset = rasterio.open(tile_path)
-    row, col = dataset.index(x, y)
-    band = dataset.read(1)
-    dataset.close()
-    return band[row, col]
-
-
 def determine_num_samples(distance_m):
     """Guarantee a number of samples between 2 and 600.
 
     Longley-Rice Irregular Terrain Model is limited to only 600
     elevation points, so this function ensures this number is not
     passed.
-
     """
+    # This is -1/x translated and rescaled:
+    # - to hit 2 at x=0
+    # - to approach 600 as x->inf
+
     # Online plot to get a feel for the function:
     # https://www.wolframalpha.com/input/?i=plot+%28-1%2F%280.0001x%2B%281%2F598%29%29%29+%2B+600+from+0+to+100