From 5c5aa188da65847103c7679e110566a03dc446a8 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Sat, 29 Feb 2020 22:01:15 -0500 Subject: [PATCH 1/8] Grass engine refactoring, volume calculations using python --- app/plugins/grass_engine.py | 40 +++++--------- plugins/measure/api.py | 2 +- plugins/measure/calc_volume.grass | 32 ----------- plugins/measure/calc_volume.py | 89 +++++++++++++++++++++++++++++++ 4 files changed, 103 insertions(+), 60 deletions(-) delete mode 100755 plugins/measure/calc_volume.grass create mode 100755 plugins/measure/calc_volume.py diff --git a/app/plugins/grass_engine.py b/app/plugins/grass_engine.py index be89043d..b7cd4854 100644 --- a/app/plugins/grass_engine.py +++ b/app/plugins/grass_engine.py @@ -4,8 +4,6 @@ import tempfile import subprocess import os -from string import Template - from webodm import settings logger = logging.getLogger('app.logger') @@ -30,12 +28,12 @@ class GrassEngine: class GrassContext: - def __init__(self, grass_binary, tmpdir = None, template_args = {}, location = None, auto_cleanup=True): + def __init__(self, grass_binary, tmpdir = None, script_opts = {}, location = None, auto_cleanup=True): self.grass_binary = grass_binary if tmpdir is None: tmpdir = os.path.basename(tempfile.mkdtemp('_grass_engine', dir=settings.MEDIA_TMP)) self.tmpdir = tmpdir - self.template_args = template_args + self.script_opts = script_opts self.location = location self.auto_cleanup = auto_cleanup @@ -48,15 +46,15 @@ class GrassContext: dst_path = os.path.abspath(os.path.join(self.get_cwd(), filename)) with open(dst_path, 'w') as f: f.write(source) - self.template_args[param] = dst_path + self.script_opts[param] = dst_path if use_as_location: - self.set_location(self.template_args[param]) + self.set_location(self.script_opts[param]) return dst_path def add_param(self, param, value): - self.template_args[param] = value + self.script_opts[param] = value def set_location(self, location): """ @@ -75,25 +73,12 @@ class GrassContext: script = os.path.abspath(script) - # Create grass script via template substitution - try: - with open(script) as f: - script_content = f.read() - except FileNotFoundError: - raise GrassEngineException("Script does not exist: {}".format(script)) - - tmpl = Template(script_content) - - # Write script to disk - if not os.path.exists(self.get_cwd()): - os.mkdir(self.get_cwd()) - - with open(os.path.join(self.get_cwd(), 'script.sh'), 'w') as f: - f.write(tmpl.substitute(self.template_args)) + # Create param list + params = ["{}={}".format(opt,value) for opt,value in self.script_opts.items()] # Execute it - logger.info("Executing grass script from {}: {} -c {} location --exec sh script.sh".format(self.get_cwd(), self.grass_binary, self.location)) - p = subprocess.Popen([self.grass_binary, '-c', self.location, 'location', '--exec', 'sh', 'script.sh'], + logger.info("Executing grass script from {}: {} -c {} location --exec python {} {}".format(self.get_cwd(), self.grass_binary, self.location, script, " ".join(params))) + p = subprocess.Popen([self.grass_binary, '-c', self.location, 'location', '--exec', 'python', script] + params, cwd=self.get_cwd(), stdout=subprocess.PIPE, stderr=subprocess.PIPE) out, err = p.communicate() @@ -108,14 +93,15 @@ class GrassContext: def serialize(self): return { 'tmpdir': self.tmpdir, - 'template_args': self.template_args, + 'script_opts': self.script_opts, 'location': self.location, 'auto_cleanup': self.auto_cleanup } def cleanup(self): - if os.path.exists(self.get_cwd()): - shutil.rmtree(self.get_cwd()) + pass + # if os.path.exists(self.get_cwd()): + # shutil.rmtree(self.get_cwd()) def __del__(self): if self.auto_cleanup: diff --git a/plugins/measure/api.py b/plugins/measure/api.py index b45f5e4e..a58ed41d 100644 --- a/plugins/measure/api.py +++ b/plugins/measure/api.py @@ -41,7 +41,7 @@ class TaskVolume(TaskView): celery_task_id = execute_grass_script.delay(os.path.join( os.path.dirname(os.path.abspath(__file__)), - "calc_volume.grass" + "calc_volume.py" ), context.serialize()).task_id return Response({'celery_task_id': celery_task_id}, status=status.HTTP_200_OK) diff --git a/plugins/measure/calc_volume.grass b/plugins/measure/calc_volume.grass deleted file mode 100755 index cfe9886c..00000000 --- a/plugins/measure/calc_volume.grass +++ /dev/null @@ -1,32 +0,0 @@ -# area_file: Geospatial file containing the area to measure -# points_file: Geospatial file containing the points defining the area -# dsm_file: GeoTIFF DEM containing the surface -# ------ -# output: prints the volume to stdout - -#Import raster and vector -v.import input=${area_file} output=polygon_area --overwrite -v.import input=${points_file} output=polygon_points --overwrite -v.buffer -s --overwrite input=polygon_area type=area output=region distance=1 minordistance=1 -r.external input=${dsm_file} output=dsm --overwrite - -# Set Grass region to vector bbox -g.region vector=region - -# Create a mask to speed up computation -r.mask vect=region - -# Transfer dsm raster data to vector -v.what.rast map=polygon_points raster=dsm column=height - -# Decimate DSM and generate interpolation of new terrain -v.surf.rst --overwrite input=polygon_points zcolumn=height elevation=dsm_below_pile - -# Compute difference between dsm and new dsm -r.mapcalc expression='pile_height_above_dsm=dsm-dsm_below_pile' --overwrite - -# Set region to polygon area to calculate volume -g.region vect=polygon_area - -# Volume output from difference -r.volume -f input=pile_height_above_dsm diff --git a/plugins/measure/calc_volume.py b/plugins/measure/calc_volume.py new file mode 100755 index 00000000..16cb2c6e --- /dev/null +++ b/plugins/measure/calc_volume.py @@ -0,0 +1,89 @@ +#%module +#% description: Calculate volume of area and prints the volume to stdout +#%end +#%option +#% key: area_file +#% type: string +#% required: yes +#% multiple: no +#% description: Geospatial file containing the area to measure +#%end +#%option +#% key: points_file +#% type: string +#% required: yes +#% multiple: no +#% description: Geospatial file containing the points defining the area +#%end +#%option +#% key: dsm_file +#% type: string +#% required: yes +#% multiple: no +#% description: GeoTIFF DEM containing the surface +#%end + +import sys +from grass.pygrass.modules import Module +import grass.script as grass + +def main(): + # Import raster and vector + Module("v.import", input=opts['area_file'], output="polygon_area", overwrite=True) + Module("v.import", input=opts['points_file'], output="polygon_points", overwrite=True) + Module("v.buffer", input="polygon_area", s=True, type="area", output="region", distance=1, minordistance=1, overwrite=True) + Module("r.external", input=opts['dsm_file'], output="dsm", overwrite=True) + + # Set Grass region to vector bbox + Module("g.region", vector="region") + + # Create a mask to speed up computation + Module("r.mask", vector="region") + + # Transfer dsm raster data to vector + Module("v.what.rast", map="polygon_points", raster="dsm", column="height") + + # Decimate DSM and generate interpolation of new terrain + Module("v.surf.rst", input="polygon_points", zcolumn="height", elevation="dsm_below_pile", overwrite=True) + + # Compute difference between dsm and new dsm + Module("r.mapcalc", expression='pile_height_above_dsm=dsm-dsm_below_pile', overwrite=True) + + # Set region to polygon area to calculate volume + Module("g.region", vector="polygon_area") + + # Volume output from difference + Module("r.volume", input="pile_height_above_dsm", f=True) + + return 0 + +if __name__ == "__main__": + opts, _ = grass.parser() + sys.exit(main()) + +#Import raster and vector +# v.import input=${area_file} output=polygon_area --overwrite +# v.import input=${points_file} output=polygon_points --overwrite +# v.buffer -s --overwrite input=polygon_area type=area output=region distance=1 minordistance=1 +# r.external input=${dsm_file} output=dsm --overwrite +# +# # Set Grass region to vector bbox +# g.region vector=region +# +# # Create a mask to speed up computation +# r.mask vect=region +# +# # Transfer dsm raster data to vector +# v.what.rast map=polygon_points raster=dsm column=height +# +# # Decimate DSM and generate interpolation of new terrain +# v.surf.rst --overwrite input=polygon_points zcolumn=height elevation=dsm_below_pile +# +# # Compute difference between dsm and new dsm +# r.mapcalc expression='pile_height_above_dsm=dsm-dsm_below_pile' --overwrite +# +# # Set region to polygon area to calculate volume +# g.region vect=polygon_area +# +# # Volume output from difference +# r.volume -f input=pile_height_above_dsm From 40f370457d73792ab3a96065ddab6739d1d60055 Mon Sep 17 00:00:00 2001 From: Sylvain POULAIN Date: Sun, 1 Mar 2020 13:16:53 +0400 Subject: [PATCH 2/8] Update region resolution Update region resolution to have same resolution as DSM resolution. Allows more accurate measurement --- plugins/measure/calc_volume.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/measure/calc_volume.py b/plugins/measure/calc_volume.py index 16cb2c6e..94a2fad5 100755 --- a/plugins/measure/calc_volume.py +++ b/plugins/measure/calc_volume.py @@ -35,7 +35,7 @@ def main(): Module("r.external", input=opts['dsm_file'], output="dsm", overwrite=True) # Set Grass region to vector bbox - Module("g.region", vector="region") + Module("g.region", vector="region", res="dsm") # Create a mask to speed up computation Module("r.mask", vector="region") From 5fd4ce0f538af23e64cc47ded2135ab418e83e47 Mon Sep 17 00:00:00 2001 From: Sylvain POULAIN Date: Sun, 1 Mar 2020 13:19:23 +0400 Subject: [PATCH 3/8] Update comment --- plugins/measure/calc_volume.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/measure/calc_volume.py b/plugins/measure/calc_volume.py index 94a2fad5..c4e3f9db 100755 --- a/plugins/measure/calc_volume.py +++ b/plugins/measure/calc_volume.py @@ -34,7 +34,7 @@ def main(): Module("v.buffer", input="polygon_area", s=True, type="area", output="region", distance=1, minordistance=1, overwrite=True) Module("r.external", input=opts['dsm_file'], output="dsm", overwrite=True) - # Set Grass region to vector bbox + # Set Grass region to vector bbox and resolution to DSM Module("g.region", vector="region", res="dsm") # Create a mask to speed up computation From 2fce053f2407e140b359c0588c542b0315433843 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Sun, 1 Mar 2020 22:33:13 -0500 Subject: [PATCH 4/8] Started rewriting contours grass script --- plugins/contours/calc_contours.grass | 50 --------------- plugins/contours/calc_contours.py | 95 ++++++++++++++++++++++++++++ plugins/measure/calc_volume.py | 27 -------- 3 files changed, 95 insertions(+), 77 deletions(-) delete mode 100755 plugins/contours/calc_contours.grass create mode 100755 plugins/contours/calc_contours.py diff --git a/plugins/contours/calc_contours.grass b/plugins/contours/calc_contours.grass deleted file mode 100755 index 870f94bd..00000000 --- a/plugins/contours/calc_contours.grass +++ /dev/null @@ -1,50 +0,0 @@ -# dem_file: GeoTIFF DEM containing the surface to calculate contours -# interval: Contours interval -# format: OGR output format -# simplify: Simplify value -# epsg: target EPSG code -# -# ------ -# output: If successful, prints the full path to the contours file. Otherwise it prints "error" - -ext="" -if [ "${format}" = "GeoJSON" ]; then - ext="json" -elif [ "${format}" = "GPKG" ]; then - ext="gpkg" -elif [ "${format}" = "DXF" ]; then - ext="dxf" -elif [ "${format}" = "ESRI Shapefile" ]; then - ext="shp" -fi - -#gdal_contour -a elevation -i ${interval} -f GPKG "${dem_file}" contours.gpkg > /dev/null -#ogr2ogr -dialect SQLite -where "ST_Length(geom) > 4" -simplify ${simplify} -t_srs EPSG:${epsg} -overwrite -f "${format}" output.$$ext contours.gpkg > /dev/null - -MIN_CONTOUR_LENGTH=5 - -r.external input="${dem_file}" output=dem --overwrite -g.region raster=dem -r.contour input=dem output=contours step=${interval} --overwrite -v.generalize --overwrite input=contours output=contours_smooth method=douglas threshold=${simplify} -v.generalize --overwrite input=contours_smooth output=contours_simplified method=chaiken threshold=1 -v.generalize --overwrite input=contours_simplified output=contours_final method=douglas threshold=${simplify} -v.edit map=contours_final tool=delete threshold=-1,0,-$$MIN_CONTOUR_LENGTH query=length -v.out.ogr input=contours_final output=temp.gpkg format="GPKG" - -ogr2ogr -t_srs EPSG:${epsg} -overwrite -f "${format}" output.$$ext temp.gpkg > /dev/null - -if [ -e "output.$$ext" ]; then - # ESRI ShapeFile extra steps to compress into a zip archive - # we leverage Python's shutil in this case - if [ "${format}" = "ESRI Shapefile" ]; then - ext="zip" - mkdir contours/ - mv output* contours/ - echo "import shutil;shutil.make_archive('output', 'zip', 'contours/')" | python - fi - - echo "$$(pwd)/output.$$ext" -else - echo "error" -fi diff --git a/plugins/contours/calc_contours.py b/plugins/contours/calc_contours.py new file mode 100755 index 00000000..39004838 --- /dev/null +++ b/plugins/contours/calc_contours.py @@ -0,0 +1,95 @@ +#%module +#% description: Calculate contours +#%end +#%option +#% key: dem_file +#% type: string +#% required: yes +#% multiple: no +#% description: GeoTIFF DEM containing the surface to calculate contours +#%end +#%option +#% key: interval +#% type: double +#% required: yes +#% multiple: no +#% description: Contours interval +#%end +#%option +#% key: format +#% type: string +#% required: yes +#% multiple: no +#% description: OGR output format +#%end +#%option +#% key: simplify +#% type: double +#% required: yes +#% multiple: no +#% description: OGR output format +#%end +#%option +#% key: epsg +#% type: string +#% required: yes +#% multiple: no +#% description: target EPSG code +#%end + +# output: If successful, prints the full path to the contours file. Otherwise it prints "error" + +import sys +from grass.pygrass.modules import Module +import grass.script as grass + +def main(): + ext = "" + if opts['format'] == "GeoJSON": + ext = "json" + elif opts['format'] == "GPKG": + ext = "gpkg" + elif opts['format'] == "DXF": + ext = "dxf" + elif opts['format'] == "ESRI Shapefile": + ext = "shp" + + MIN_CONTOUR_LENGTH = 5 + Module("r.external", input=opts['dem_file'], output="dem", overwrite=True) + Module("g.region", raster="dem") + Module("r.contour", input="dem", output="contours", step=opts["interval"], overwrite=True) + Module("r.generalize", input="contours", output="contours_smooth", method="douglas", threshold=opts["simplify"], overwrite=True) + Module("r.generalize", input="contours_smooth", output="contours_simplified", method="chaiken", threshold=1, overwrite=True) + Module("r.generalize", input="contours_simplified", output="contours_final", method="douglas", threshold=opts["simplify"], overwrite=True) + Module("v.edit", input="contours_final", tool="delete", threshold="-1,0,-%s" % MIN_CONTOUR_LENGTH, query="length") + Module("v.out.ogr", input="contours_final", output="temp.gpkg", format="GPKG") + + return 0 + +if __name__ == "__main__": + opts, _ = grass.parser() + sys.exit(main()) + + +# TODO +# Running external commands from Python +# For information on running external commands from Python, see: http://docs.python.org/lib/module-subprocess.html + +# Avoid using the older os.* functions. Section 17.1.3 lists equivalents using the Popen() interface, which is more robust (particularly on Windows). +# +# ogr2ogr -t_srs EPSG:${epsg} -overwrite -f "${format}" output.$$ext temp.gpkg > /dev/null +# +# if [ -e "output.$$ext" ]; then +# # ESRI ShapeFile extra steps to compress into a zip archive +# # we leverage Python's shutil in this case +# if [ "${format}" = "ESRI Shapefile" ]; then +# ext="zip" +# mkdir contours/ +# mv output* contours/ +# echo "import shutil;shutil.make_archive('output', 'zip', 'contours/')" | python +# fi +# +# echo "$$(pwd)/output.$$ext" +# else +# echo "error" +# fi diff --git a/plugins/measure/calc_volume.py b/plugins/measure/calc_volume.py index c4e3f9db..8deb2673 100755 --- a/plugins/measure/calc_volume.py +++ b/plugins/measure/calc_volume.py @@ -60,30 +60,3 @@ def main(): if __name__ == "__main__": opts, _ = grass.parser() sys.exit(main()) - -#Import raster and vector -# v.import input=${area_file} output=polygon_area --overwrite -# v.import input=${points_file} output=polygon_points --overwrite -# v.buffer -s --overwrite input=polygon_area type=area output=region distance=1 minordistance=1 -# r.external input=${dsm_file} output=dsm --overwrite -# -# # Set Grass region to vector bbox -# g.region vector=region -# -# # Create a mask to speed up computation -# r.mask vect=region -# -# # Transfer dsm raster data to vector -# v.what.rast map=polygon_points raster=dsm column=height -# -# # Decimate DSM and generate interpolation of new terrain -# v.surf.rst --overwrite input=polygon_points zcolumn=height elevation=dsm_below_pile -# -# # Compute difference between dsm and new dsm -# r.mapcalc expression='pile_height_above_dsm=dsm-dsm_below_pile' --overwrite -# -# # Set region to polygon area to calculate volume -# g.region vect=polygon_area -# -# # Volume output from difference -# r.volume -f input=pile_height_above_dsm From b542baca16b7b060ea0a4e5a53f5f7af3d63c7d3 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 11 Mar 2020 22:06:04 -0400 Subject: [PATCH 5/8] Contours script rewrite --- plugins/contours/api.py | 2 +- plugins/contours/calc_contours.py | 52 +++++++++++++++---------------- 2 files changed, 26 insertions(+), 28 deletions(-) diff --git a/plugins/contours/api.py b/plugins/contours/api.py index b1a3b926..def7b764 100644 --- a/plugins/contours/api.py +++ b/plugins/contours/api.py @@ -42,7 +42,7 @@ class TaskContoursGenerate(TaskView): celery_task_id = execute_grass_script.delay(os.path.join( os.path.dirname(os.path.abspath(__file__)), - "calc_contours.grass" + "calc_contours.py" ), context.serialize(), 'file').task_id return Response({'celery_task_id': celery_task_id}, status=status.HTTP_200_OK) diff --git a/plugins/contours/calc_contours.py b/plugins/contours/calc_contours.py index 39004838..987ab5bb 100755 --- a/plugins/contours/calc_contours.py +++ b/plugins/contours/calc_contours.py @@ -40,8 +40,12 @@ # output: If successful, prints the full path to the contours file. Otherwise it prints "error" import sys +import glob +import os +import shutil from grass.pygrass.modules import Module import grass.script as grass +import subprocess def main(): ext = "" @@ -58,38 +62,32 @@ def main(): Module("r.external", input=opts['dem_file'], output="dem", overwrite=True) Module("g.region", raster="dem") Module("r.contour", input="dem", output="contours", step=opts["interval"], overwrite=True) - Module("r.generalize", input="contours", output="contours_smooth", method="douglas", threshold=opts["simplify"], overwrite=True) - Module("r.generalize", input="contours_smooth", output="contours_simplified", method="chaiken", threshold=1, overwrite=True) - Module("r.generalize", input="contours_simplified", output="contours_final", method="douglas", threshold=opts["simplify"], overwrite=True) - Module("v.edit", input="contours_final", tool="delete", threshold="-1,0,-%s" % MIN_CONTOUR_LENGTH, query="length") + Module("v.generalize", input="contours", output="contours_smooth", method="douglas", threshold=opts["simplify"], overwrite=True) + Module("v.generalize", input="contours_smooth", output="contours_simplified", method="chaiken", threshold=1, overwrite=True) + Module("v.generalize", input="contours_simplified", output="contours_final", method="douglas", threshold=opts["simplify"], overwrite=True) + Module("v.edit", map="contours_final", tool="delete", threshold=[-1,0,-MIN_CONTOUR_LENGTH], query="length") Module("v.out.ogr", input="contours_final", output="temp.gpkg", format="GPKG") + subprocess.check_call(["ogr2ogr", "-t_srs", "EPSG:%s" % opts['epsg'], + '-overwrite', '-f', opts["format"], "output.%s" % ext, "temp.gpkg"], stdout=subprocess.DEVNULL) + + if os.path.isfile("output.%s" % ext): + if opts["format"] == "ESRI Shapefile": + ext="zip" + os.makedirs("contours") + contour_files = glob.glob("output.*") + for cf in contour_files: + shutil.move(cf, os.path.join("contours", os.path.basename(cf))) + + shutil.make_archive('output', 'zip', 'contours/') + + print(os.path.join(os.getcwd(), "output.%s" % ext)) + else: + print("error") + return 0 if __name__ == "__main__": opts, _ = grass.parser() sys.exit(main()) - -# TODO -# Running external commands from Python -# For information on running external commands from Python, see: http://docs.python.org/lib/module-subprocess.html - -# Avoid using the older os.* functions. Section 17.1.3 lists equivalents using the Popen() interface, which is more robust (particularly on Windows). -# -# ogr2ogr -t_srs EPSG:${epsg} -overwrite -f "${format}" output.$$ext temp.gpkg > /dev/null -# -# if [ -e "output.$$ext" ]; then -# # ESRI ShapeFile extra steps to compress into a zip archive -# # we leverage Python's shutil in this case -# if [ "${format}" = "ESRI Shapefile" ]; then -# ext="zip" -# mkdir contours/ -# mv output* contours/ -# echo "import shutil;shutil.make_archive('output', 'zip', 'contours/')" | python -# fi -# -# echo "$$(pwd)/output.$$ext" -# else -# echo "error" -# fi From d6f0f47a0a09cea52d0ffd967e56fe8ee4c3df4a Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 11 Mar 2020 22:31:37 -0400 Subject: [PATCH 6/8] Adapted elevation map plugin to new GRASS engine format --- plugins/elevationmap/api.py | 11 +- plugins/elevationmap/calc_elevation_map.grass | 45 ------- plugins/elevationmap/elevationmap.py | 119 ++++++++++++++---- 3 files changed, 97 insertions(+), 78 deletions(-) delete mode 100755 plugins/elevationmap/calc_elevation_map.grass diff --git a/plugins/elevationmap/api.py b/plugins/elevationmap/api.py index 13e9914e..5dd49e02 100644 --- a/plugins/elevationmap/api.py +++ b/plugins/elevationmap/api.py @@ -37,21 +37,18 @@ class TaskElevationMapGenerate(TaskView): noise_filter_size = float(request.data.get('noise_filter_size', 2)) current_dir = os.path.dirname(os.path.abspath(__file__)) - context.add_param('dsm_file', dsm) + context.add_param('dsm', dsm) context.add_param('interval', interval) context.add_param('format', format) context.add_param('noise_filter_size', noise_filter_size) context.add_param('epsg', epsg) - context.add_param('python_script_path', os.path.join(current_dir, "elevationmap.py")) - context.add_param('python_path', plugin.get_python_packages_path()) if dtm != None: - context.add_param('dtm', '--dtm {}'.format(dtm)) - else: - context.add_param('dtm', '') + context.add_param('dtm', dtm) + context.set_location(dsm) - celery_task_id = execute_grass_script.delay(os.path.join(current_dir, "calc_elevation_map.grass"), context.serialize()).task_id + celery_task_id = execute_grass_script.delay(os.path.join(current_dir, "elevationmap.py"), context.serialize()).task_id return Response({'celery_task_id': celery_task_id}, status=status.HTTP_200_OK) except GrassEngineException as e: diff --git a/plugins/elevationmap/calc_elevation_map.grass b/plugins/elevationmap/calc_elevation_map.grass deleted file mode 100755 index a23accd5..00000000 --- a/plugins/elevationmap/calc_elevation_map.grass +++ /dev/null @@ -1,45 +0,0 @@ -# dsm_file: GeoTIFF DSM containing the surface to calculate the levels of elevation -# interval: Interval that will split the different levels -# format: OGR output format -# noise_filter_size: Area in meters where we will clean up noise in the contours -# epsg: target EPSG code -# python_script_path: Path of the python script -# python_path: Path to python modules -# dtm: Optional text to include the GeoTIFF DTM -# -# ------ -# output: If successful, prints the full path to the elevation map file. If an error.txt file is found, print its content. Otherwise, print 'error'" - -ext="" -if [ "${format}" = "GeoJSON" ]; then - ext="json" -elif [ "${format}" = "GPKG" ]; then - ext="gpkg" -elif [ "${format}" = "DXF" ]; then - ext="dxf" -elif [ "${format}" = "ESRI Shapefile" ]; then - ext="shp" -fi - -PYTHONPATH="${python_path}" "${python_script_path}" "${dsm_file}" ${interval} --epsg ${epsg} --noise_filter_size ${noise_filter_size} ${dtm} -o output.json - -if [ $$ext != "json" ]; then - ogr2ogr -f "${format}" output.$$ext output.json > /dev/null -fi - -if [ -e "output.$$ext" ]; then - # ESRI ShapeFile extra steps to compress into a zip archive - # we leverage Python's shutil in this case - if [ "${format}" = "ESRI Shapefile" ]; then - ext="zip" - mkdir elevationmap/ - mv output* elevationmap/ - echo "import shutil;shutil.make_archive('output', 'zip', 'elevationmap/')" | python - fi - - echo "$$(pwd)/output.$$ext" -elif [ -e "error.txt" ]; then - cat "error.txt" -else - echo "error" -fi diff --git a/plugins/elevationmap/elevationmap.py b/plugins/elevationmap/elevationmap.py index ee0107e8..dd00bf4d 100755 --- a/plugins/elevationmap/elevationmap.py +++ b/plugins/elevationmap/elevationmap.py @@ -1,23 +1,84 @@ -#!/usr/bin/env python3 +#%module +#% description: This script takes a GeoTIFF file, calculates its heighmap, and outputs it as a GeoJSON +#%end +#%option +#% key: dsm +#% type: string +#% required: yes +#% multiple: no +#% description: The path for the dsm file +#%end +#%option +#% key: intervals +#% type: double +#% required: yes +#% multiple: no +#% description: The intervals used to generate the diferent elevation levels +#%end +#%option +#% key: format +#% type: string +#% required: yes +#% multiple: no +#% description: OGR output format +#%end +#%option +#% key: dtm +#% type: string +#% required: no +#% multiple: no +#% description: The path for the dtm file +#%end +#%option +#% key: epsg +#% type: string +#% required: yes +#% multiple: no +#% description: The epsg code that will be used for output +#%end +#%option +#% key: noise_filter_size +#% type: double +#% required: yes +#% multiple: no +#% description: Area in meters where we will clean up noise in the contours +#%end + import cv2, math, argparse import numpy as np import rasterio as rio from rasterio import warp, transform -from geojson import Feature, FeatureCollection, MultiPolygon, dumps +from geojson import Feature, FeatureCollection, MultiPolygon, dumps +import subprocess +import os +import glob +import shutil +import sys +import grass.script as grass + +def main(): + ext = "" + if opts['format'] == "GeoJSON": + ext = "json" + elif opts['format'] == "GPKG": + ext = "gpkg" + elif opts['format'] == "DXF": + ext = "dxf" + elif opts['format'] == "ESRI Shapefile": + ext = "shp" -def main(args): # Open dsm - dsm = rio.open(args.dsm) + dsm = rio.open(opts['dsm']) # Read the tiff as an numpy masked array dsm_array = dsm.read(1, masked = True) # Create a kernel based on the parameter 'noise_filter_size' and the tiff resolution - kernel = get_kernel(args.noise_filter_size, dsm) + kernel = get_kernel(float(opts['noise_filter_size']), dsm) # Check if we want to use the dtm also - if args.dtm != None: + if opts['dtm'] != '': # Open the dtm - dtm = rio.open(args.dtm) + dtm = rio.open(opts['dtm']) # Assert that the dtm and dsm have the same bounds and resolution assert_same_bounds_and_resolution(dsm, dtm) # Calculate the different between the dsm and dtm @@ -26,7 +87,7 @@ def main(args): array = dsm_array # Calculate the ranges based on the parameter 'intervals' and the elevation array - ranges = calculate_ranges(args.intervals, array) + ranges = calculate_ranges(opts['intervals'], array) features = [] @@ -40,15 +101,32 @@ def main(args): # Check if we found something if len(contours) > 0: # Transform contours from pixels to coordinates - mapped_contours = [map_pixels_to_coordinates(dsm, args.epsg, to_pixel_format(contour)) for contour in contours] + mapped_contours = [map_pixels_to_coordinates(dsm, opts['epsg'], to_pixel_format(contour)) for contour in contours] # Build the MultiPolygon for based on the contours and their hierarchy built_multi_polygon = LevelBuilder(bottom, top, mapped_contours, hierarchy[0]).build_multi_polygon() features.append(built_multi_polygon) # Write the GeoJSON to a file dump = dumps(FeatureCollection(features)) - with open(args.output, 'w+') as output: - output.write(dump) + with open("output.json", 'w+') as output: + output.write(dump) + + if ext != "json": + subprocess.check_call(["ogr2ogr", "-f", opts['format'], "output.%s" % ext, "output.json"], stdout=subprocess.DEVNULL) + + if os.path.isfile("output.%s" % ext): + if opts['format'] == "ESRI Shapefile": + ext="zip" + os.makedirs("contours") + contour_files = glob.glob("output.*") + for cf in contour_files: + shutil.move(cf, os.path.join("contours", os.path.basename(cf))) + + shutil.make_archive('output', 'zip', 'contours/') + + print(os.path.join(os.getcwd(), "output.%s" % ext)) + else: + print("error") def get_kernel(noise_filter_size, dsm): """Generate a kernel for noise filtering. Will return none if the noise_filter_size isn't positive""" @@ -151,18 +229,7 @@ class LevelBuilder: multi_polygon = MultiPolygon(polygons) return Feature(geometry = multi_polygon, properties = { 'bottom': int(self.bottom), 'top': int(self.top) }) -if __name__ == '__main__': - parser = argparse.ArgumentParser(description = 'This script takes a GeoTIFF file, calculates its heighmap, and outputs it as a GeoJSON') - parser.add_argument('dsm', type = str, help = 'The path for the dsm file') - parser.add_argument('intervals', type = str, help = 'The intervals used to generate the diferent elevation levels') - parser.add_argument('-d', '--dtm', type = str, help = 'The path for the dtm file') - parser.add_argument('-e', '--epsg', type = int, help = 'The epsg code that will be used for output', default = 4326) - parser.add_argument('-k', '--noise_filter_size', type = float, help = 'Area in meters where we will clean up noise in the contours', default = 2) - parser.add_argument('-o', '--output', type = str, help = 'The path for the output file', default = "output.json") - - args = parser.parse_args() - try: - main(args) - except Exception as e: - with open('error.txt', 'w+') as output: - output.write(str(e)) \ No newline at end of file + +if __name__ == "__main__": + opts, _ = grass.parser() + sys.exit(main()) From f3b761e46d9475762872b9e2ae9754b588507dbf Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 12 Mar 2020 11:16:32 -0400 Subject: [PATCH 7/8] Fixed tests --- app/plugins/grass_engine.py | 9 +++++--- app/tests/grass_scripts/simple_test.grass | 6 ------ app/tests/grass_scripts/simple_test.py | 25 +++++++++++++++++++++++ app/tests/test_plugins.py | 13 ++++++------ 4 files changed, 38 insertions(+), 15 deletions(-) delete mode 100644 app/tests/grass_scripts/simple_test.grass create mode 100644 app/tests/grass_scripts/simple_test.py diff --git a/app/plugins/grass_engine.py b/app/plugins/grass_engine.py index b7cd4854..7695766b 100644 --- a/app/plugins/grass_engine.py +++ b/app/plugins/grass_engine.py @@ -73,6 +73,10 @@ class GrassContext: script = os.path.abspath(script) + # Make sure working directory exists + if not os.path.exists(self.get_cwd()): + os.mkdir(self.get_cwd()) + # Create param list params = ["{}={}".format(opt,value) for opt,value in self.script_opts.items()] @@ -99,9 +103,8 @@ class GrassContext: } def cleanup(self): - pass - # if os.path.exists(self.get_cwd()): - # shutil.rmtree(self.get_cwd()) + if os.path.exists(self.get_cwd()): + shutil.rmtree(self.get_cwd()) def __del__(self): if self.auto_cleanup: diff --git a/app/tests/grass_scripts/simple_test.grass b/app/tests/grass_scripts/simple_test.grass deleted file mode 100644 index 732810c2..00000000 --- a/app/tests/grass_scripts/simple_test.grass +++ /dev/null @@ -1,6 +0,0 @@ -# test: Geospatial test file -# ------ -# output: greets the user and prints the information of a spatial file - -v.in.ogr input=${test} layer=test output=test --overwrite -v.info map=test diff --git a/app/tests/grass_scripts/simple_test.py b/app/tests/grass_scripts/simple_test.py new file mode 100644 index 00000000..ed7d3e29 --- /dev/null +++ b/app/tests/grass_scripts/simple_test.py @@ -0,0 +1,25 @@ +#%module +#% description: greets the user and prints the information of a spatial file +#%end +#%option +#% key: test +#% type: string +#% required: yes +#% multiple: no +#% description: Geospatial test file +#%end + +import sys +from grass.pygrass.modules import Module +import grass.script as grass + +def main(): + # Import raster and vector + Module("v.in.ogr", input=opts['test'], layer="test", output="test", overwrite=True) + info = grass.vector_info("test") + print("Number of points: %s" % info['points']) + +if __name__ == "__main__": + opts, _ = grass.parser() + sys.exit(main()) + diff --git a/app/tests/test_plugins.py b/app/tests/test_plugins.py index 4cb1b427..651887f5 100644 --- a/app/tests/test_plugins.py +++ b/app/tests/test_plugins.py @@ -149,10 +149,11 @@ class TestPlugins(BootTestCase): ctx.set_location("EPSG:4326") result = execute_grass_script.delay( - os.path.join(grass_scripts_dir, "simple_test.grass"), + os.path.join(grass_scripts_dir, "simple_test.py"), ctx.serialize() ).get() - self.assertTrue("Number of points: 1" in result.get('output')) + + self.assertEqual("Number of points: 1", result.get('output')) self.assertTrue(result.get('context') == ctx.serialize()) @@ -160,24 +161,24 @@ class TestPlugins(BootTestCase): self.assertFalse(os.path.exists(ctx.get_cwd())) error = execute_grass_script.delay( - os.path.join(grass_scripts_dir, "nonexistant_script.grass"), + os.path.join(grass_scripts_dir, "nonexistant_script.py"), ctx.serialize() ).get() self.assertIsInstance(error, dict) self.assertIsInstance(error['error'], str) with self.assertRaises(GrassEngineException): - ctx.execute(os.path.join(grass_scripts_dir, "nonexistant_script.grass")) + ctx.execute(os.path.join(grass_scripts_dir, "nonexistant_script.py")) ctx = grass.create_context({"auto_cleanup": False}) ctx.add_file('test.geojson', points) ctx.set_location("EPSG:4326") result = execute_grass_script.delay( - os.path.join(grass_scripts_dir, "simple_test.grass"), + os.path.join(grass_scripts_dir, "simple_test.py"), ctx.serialize() ).get() - self.assertTrue("Number of points: 1" in result.get('output')) + self.assertEqual("Number of points: 1", result.get('output')) # Path still there self.assertTrue(os.path.exists(ctx.get_cwd())) From cb6fbca5ef08a13e2ffc2a9955ba701325c23d3c Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Fri, 13 Mar 2020 18:03:21 -0400 Subject: [PATCH 8/8] Fix #833 --- plugins/measure/calc_volume.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/measure/calc_volume.py b/plugins/measure/calc_volume.py index 8deb2673..4a3fbfef 100755 --- a/plugins/measure/calc_volume.py +++ b/plugins/measure/calc_volume.py @@ -35,7 +35,7 @@ def main(): Module("r.external", input=opts['dsm_file'], output="dsm", overwrite=True) # Set Grass region to vector bbox and resolution to DSM - Module("g.region", vector="region", res="dsm") + Module("g.region", vector="region") # Create a mask to speed up computation Module("r.mask", vector="region")