GRASS script to create DEM, contour and textured relief maps.

Signed-off-by: Mark Hale <mark.hale@physics.org>
Former-commit-id: 9ab6d4a43c

contrib/grasss/odm_grass.py

@@ -0,0 +1,131 @@
+#!/usr/bin/env python
+
+# To run, set the following env variables:
+# PYTHONHOME location of Python
+# PYTHONPATH location of GRASS Python libs
+# PATH include GRASS bin and lib
+# GISBASE location of GRASS
+
+import os
+import sys
+import grass.script as gscript
+import grass.script.core
+import grass.script.setup
+
+vlidarName = 'odm_vlidar'
+rsurfName = 'odm_rsurf'
+contourName = 'odm_contour'
+orthophotoName = 'odm_orthophoto'
+reliefName = 'odm_relief'
+shadedReliefName = reliefName + '_shaded'
+
+overwrite = True
+
+def main():
+    if len(sys.argv) < 2:
+       sys.exit('Please provide the ODM project path.')
+
+    projectHome = sys.argv[1]
+
+    gisdb = projectHome+'/grassdata'
+    location = 'odm'
+    gisrc = gscript.setup.init(os.environ['GISBASE'], gisdb, location)
+
+    # get srs and initial extents
+    with open(projectHome+'/odm_georeferencing/coords.txt') as f:
+        srs = f.readline().split()
+        mean = f.readline().split()
+        meanX = float(mean[0])
+        meanY = float(mean[1])
+        minX = float('inf')
+        maxX = float('-inf')
+        minY = float('inf')
+        maxY = float('-inf')
+        for line in f:
+           xy = line.split()
+           x = float(xy[0])
+           y = float(xy[1])
+           minX = min(x, minX)
+           maxX = max(x, maxX)
+           minY = min(y, minY)
+           maxY = max(y, maxY)
+
+    datum = srs[0]
+    proj = srs[1]
+    zone = srs[2]
+    gscript.core.create_location(gisdb, location, datum=datum, proj4='+proj='+proj+' +zone='+zone, overwrite=overwrite)
+
+    n = meanY + maxY
+    s = meanY + minY
+    e = meanX + maxX
+    w = meanX + minX
+    gscript.run_command('g.region', flags='s', n=n, s=s, e=e, w=w, res=0.01, res3=0.01, overwrite=overwrite)
+
+    dem(projectHome)
+    contour(projectHome)
+    relief(projectHome)
+
+    os.remove(gisrc)
+
+
+
+def dem(projectHome):
+    """
+    Creates a DEM in GeoTIFF format.
+    NB: this is a data raster, not an RGBA raster and so is normally only readable by GIS and not image software.
+    """
+    print 'Creating DEM'
+
+    step = 0.5
+    gscript.run_command('v.in.lidar', flags='beo', input=projectHome+'/odm_georeferencing/odm_georeferenced_model.ply.las', output=vlidarName, overwrite=overwrite)
+
+    gscript.run_command('v.surf.bspline', input=vlidarName, raster_output=rsurfName, ew_step=step, ns_step=step, method='bicubic', memory=4096, overwrite=overwrite)
+
+    gscript.run_command('r.out.gdal', flags='cf', input=rsurfName, output=projectHome+'/odm_georeferencing/odm_dem.tif', format='GTiff', type='Float32', createopt='TILED=yes,COMPRESS=DEFLATE,PREDICTOR=2,BLOCKXSIZE=512,BLOCKYSIZE=512', nodata=0, overwrite=overwrite)
+
+
+
+def contour(projectHome):
+    """
+    Creates a contour map.
+    """
+    print 'Creating contour map'
+
+    step = 0.25
+    gscript.run_command('r.contour', input=rsurfName, output=contourName, step=step, overwrite=overwrite)
+
+    gscript.run_command('v.out.ogr', input=contourName, output=projectHome+'/odm_georeferencing/odm_contour.shp', overwrite=overwrite)
+
+
+
+def relief(projectHome):
+    """
+    Creates a textured relief map in GeoTIFF format.
+    NB: this is an RGBA raster and so is readable by image software.
+    """
+    print 'Create relief map'
+
+    gscript.run_command('r.in.gdal', flags='e', input=projectHome+'/odm_orthophoto/odm_orthophoto.tif', output=orthophotoName, memory=2047, overwrite=overwrite)
+
+    gscript.run_command('r.composite', red=orthophotoName+'.red', green=orthophotoName+'.green', blue=orthophotoName+'.blue', output=orthophotoName+'.rgb', overwrite=overwrite)
+
+    gscript.run_command('r.relief', input=rsurfName, output=reliefName, overwrite=overwrite)
+
+    gscript.run_command('r.shade', shade=reliefName, color=orthophotoName+'.rgb', output=shadedReliefName, overwrite=overwrite)
+
+    calc = ';'.join([
+        '$shadedRelief.red = if(isnull($orthophoto.red), 0, r#$shadedRelief)',
+        '$shadedRelief.green = if(isnull($orthophoto.green), 0, g#$shadedRelief)',
+        '$shadedRelief.blue = if(isnull($orthophoto.blue), 0, b#$shadedRelief)',
+        '$shadedRelief.alpa = if(isnull($orthophoto.alpha), 0, 255)'
+    ])
+    gscript.mapcalc(calc, shadedRelief=shadedReliefName, orthophoto=orthophotoName, overwrite=overwrite)
+
+    gscript.run_command('i.group', group=shadedReliefName+'.group', input=shadedReliefName+'.red,'+shadedReliefName+'.green,'+shadedReliefName+'.blue,'+shadedReliefName+'.alpha')
+
+    gscript.run_command('r.out.gdal', flags='c', input=shadedReliefName+'.group', output=projectHome+'/odm_orthophoto/odm_relief.tif', format='GTiff', type='Byte', createopt='TILED=yes,COMPRESS=DEFLATE,PREDICTOR=2,BLOCKXSIZE=512,BLOCKYSIZE=512', nodata=0, overwrite=overwrite)
+
+
+
+if __name__ == '__main__':
+    main()