From 96d10d91f80032235750f0d29a1152a6d7a68ec0 Mon Sep 17 00:00:00 2001 From: Mark Hale Date: Wed, 15 Feb 2017 23:16:45 +0000 Subject: [PATCH 1/5] GRASS script to create DEM, contour and textured relief maps. Signed-off-by: Mark Hale Former-commit-id: 9ab6d4a43cd27c6cb63fbc06044a0c17b02db740 --- contrib/grasss/odm_grass.py | 131 ++++++++++++++++++++++++++++++++++++ 1 file changed, 131 insertions(+) create mode 100644 contrib/grasss/odm_grass.py diff --git a/contrib/grasss/odm_grass.py b/contrib/grasss/odm_grass.py new file mode 100644 index 00000000..0b8d989e --- /dev/null +++ b/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() From 2bb74969e7c1dac0f6ca861c52202129cda9495f Mon Sep 17 00:00:00 2001 From: Mark Hale Date: Wed, 15 Feb 2017 23:45:09 +0000 Subject: [PATCH 2/5] Corrected print msg. Former-commit-id: d1f7b2fc7609e11333b596a6fb3eefc73c7d5251 --- contrib/grasss/odm_grass.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contrib/grasss/odm_grass.py b/contrib/grasss/odm_grass.py index 0b8d989e..c7792734 100644 --- a/contrib/grasss/odm_grass.py +++ b/contrib/grasss/odm_grass.py @@ -103,7 +103,7 @@ 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' + print 'Creating relief map' gscript.run_command('r.in.gdal', flags='e', input=projectHome+'/odm_orthophoto/odm_orthophoto.tif', output=orthophotoName, memory=2047, overwrite=overwrite) From 97c509bf56ec948cb7a70dcc54b1420f844f5587 Mon Sep 17 00:00:00 2001 From: Mark Hale Date: Wed, 15 Feb 2017 23:46:53 +0000 Subject: [PATCH 3/5] Corrected file name. Signed-off-by: Mark Hale Former-commit-id: 22f16d872b0e4d414361561d7aea9202a4a7891b --- contrib/grasss/odm_grass.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contrib/grasss/odm_grass.py b/contrib/grasss/odm_grass.py index c7792734..38423cfc 100644 --- a/contrib/grasss/odm_grass.py +++ b/contrib/grasss/odm_grass.py @@ -117,7 +117,7 @@ def relief(projectHome): '$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)' + '$shadedRelief.alpha = if(isnull($orthophoto.alpha), 0, 255)' ]) gscript.mapcalc(calc, shadedRelief=shadedReliefName, orthophoto=orthophotoName, overwrite=overwrite) From f9c3c67cb1041a6dee41a68858c766cdd5e21ffd Mon Sep 17 00:00:00 2001 From: Mark Hale Date: Wed, 15 Feb 2017 23:55:55 +0000 Subject: [PATCH 4/5] Reduced geotiff file sizes. Signed-off-by: Mark Hale Former-commit-id: cfb9cae9ac39c10fe53db528f99d94a363694c7e --- contrib/grasss/odm_grass.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/contrib/grasss/odm_grass.py b/contrib/grasss/odm_grass.py index 38423cfc..46ed0e0b 100644 --- a/contrib/grasss/odm_grass.py +++ b/contrib/grasss/odm_grass.py @@ -81,7 +81,7 @@ def dem(projectHome): 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) + gscript.run_command('r.out.gdal', flags='cfm', 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) @@ -123,7 +123,7 @@ def relief(projectHome): 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) + gscript.run_command('r.out.gdal', flags='cm', 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) From a8f486b2a5314d8cf479f9df10242ef6e444e3ea Mon Sep 17 00:00:00 2001 From: Mark Hale Date: Tue, 21 Feb 2017 17:45:27 +0000 Subject: [PATCH 5/5] PEP8 compliant. Signed-off-by: Mark Hale Former-commit-id: 55781f1f33baa6231075fb9f2fbccbd026325e8b --- contrib/grass/odm_grass.py | 172 ++++++++++++++++++++++++++++++++++++ contrib/grasss/odm_grass.py | 131 --------------------------- 2 files changed, 172 insertions(+), 131 deletions(-) create mode 100644 contrib/grass/odm_grass.py delete mode 100644 contrib/grasss/odm_grass.py diff --git a/contrib/grass/odm_grass.py b/contrib/grass/odm_grass.py new file mode 100644 index 00000000..8cce77fc --- /dev/null +++ b/contrib/grass/odm_grass.py @@ -0,0 +1,172 @@ +#!/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='cfm', 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 'Creating 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.alpha = ' + + '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='cm', + 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() diff --git a/contrib/grasss/odm_grass.py b/contrib/grasss/odm_grass.py deleted file mode 100644 index 46ed0e0b..00000000 --- a/contrib/grasss/odm_grass.py +++ /dev/null @@ -1,131 +0,0 @@ -#!/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='cfm', 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 'Creating 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.alpha = 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='cm', 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()