pull/825/head
Piero Toffanin 2020-03-15 16:16:51 -04:00
commit f3c6ac6a48
13 zmienionych plików z 297 dodań i 196 usunięć

Wyświetl plik

@ -4,8 +4,6 @@ import tempfile
import subprocess import subprocess
import os import os
from string import Template
from webodm import settings from webodm import settings
logger = logging.getLogger('app.logger') logger = logging.getLogger('app.logger')
@ -30,12 +28,12 @@ class GrassEngine:
class GrassContext: 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 self.grass_binary = grass_binary
if tmpdir is None: if tmpdir is None:
tmpdir = os.path.basename(tempfile.mkdtemp('_grass_engine', dir=settings.MEDIA_TMP)) tmpdir = os.path.basename(tempfile.mkdtemp('_grass_engine', dir=settings.MEDIA_TMP))
self.tmpdir = tmpdir self.tmpdir = tmpdir
self.template_args = template_args self.script_opts = script_opts
self.location = location self.location = location
self.auto_cleanup = auto_cleanup self.auto_cleanup = auto_cleanup
@ -48,15 +46,15 @@ class GrassContext:
dst_path = os.path.abspath(os.path.join(self.get_cwd(), filename)) dst_path = os.path.abspath(os.path.join(self.get_cwd(), filename))
with open(dst_path, 'w') as f: with open(dst_path, 'w') as f:
f.write(source) f.write(source)
self.template_args[param] = dst_path self.script_opts[param] = dst_path
if use_as_location: if use_as_location:
self.set_location(self.template_args[param]) self.set_location(self.script_opts[param])
return dst_path return dst_path
def add_param(self, param, value): def add_param(self, param, value):
self.template_args[param] = value self.script_opts[param] = value
def set_location(self, location): def set_location(self, location):
""" """
@ -75,25 +73,16 @@ class GrassContext:
script = os.path.abspath(script) script = os.path.abspath(script)
# Create grass script via template substitution # Make sure working directory exists
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()): if not os.path.exists(self.get_cwd()):
os.mkdir(self.get_cwd()) os.mkdir(self.get_cwd())
with open(os.path.join(self.get_cwd(), 'script.sh'), 'w') as f: # Create param list
f.write(tmpl.substitute(self.template_args)) params = ["{}={}".format(opt,value) for opt,value in self.script_opts.items()]
# Execute it # Execute it
logger.info("Executing grass script from {}: {} -c {} location --exec sh script.sh".format(self.get_cwd(), self.grass_binary, self.location)) 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', 'sh', 'script.sh'], p = subprocess.Popen([self.grass_binary, '-c', self.location, 'location', '--exec', 'python', script] + params,
cwd=self.get_cwd(), stdout=subprocess.PIPE, stderr=subprocess.PIPE) cwd=self.get_cwd(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = p.communicate() out, err = p.communicate()
@ -108,7 +97,7 @@ class GrassContext:
def serialize(self): def serialize(self):
return { return {
'tmpdir': self.tmpdir, 'tmpdir': self.tmpdir,
'template_args': self.template_args, 'script_opts': self.script_opts,
'location': self.location, 'location': self.location,
'auto_cleanup': self.auto_cleanup 'auto_cleanup': self.auto_cleanup
} }

Wyświetl plik

@ -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

Wyświetl plik

@ -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())

Wyświetl plik

@ -149,10 +149,11 @@ class TestPlugins(BootTestCase):
ctx.set_location("EPSG:4326") ctx.set_location("EPSG:4326")
result = execute_grass_script.delay( 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() ctx.serialize()
).get() ).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()) self.assertTrue(result.get('context') == ctx.serialize())
@ -160,24 +161,24 @@ class TestPlugins(BootTestCase):
self.assertFalse(os.path.exists(ctx.get_cwd())) self.assertFalse(os.path.exists(ctx.get_cwd()))
error = execute_grass_script.delay( 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() ctx.serialize()
).get() ).get()
self.assertIsInstance(error, dict) self.assertIsInstance(error, dict)
self.assertIsInstance(error['error'], str) self.assertIsInstance(error['error'], str)
with self.assertRaises(GrassEngineException): 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 = grass.create_context({"auto_cleanup": False})
ctx.add_file('test.geojson', points) ctx.add_file('test.geojson', points)
ctx.set_location("EPSG:4326") ctx.set_location("EPSG:4326")
result = execute_grass_script.delay( 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() ctx.serialize()
).get() ).get()
self.assertTrue("Number of points: 1" in result.get('output')) self.assertEqual("Number of points: 1", result.get('output'))
# Path still there # Path still there
self.assertTrue(os.path.exists(ctx.get_cwd())) self.assertTrue(os.path.exists(ctx.get_cwd()))

Wyświetl plik

@ -42,7 +42,7 @@ class TaskContoursGenerate(TaskView):
celery_task_id = execute_grass_script.delay(os.path.join( celery_task_id = execute_grass_script.delay(os.path.join(
os.path.dirname(os.path.abspath(__file__)), os.path.dirname(os.path.abspath(__file__)),
"calc_contours.grass" "calc_contours.py"
), context.serialize(), 'file').task_id ), context.serialize(), 'file').task_id
return Response({'celery_task_id': celery_task_id}, status=status.HTTP_200_OK) return Response({'celery_task_id': celery_task_id}, status=status.HTTP_200_OK)

Wyświetl plik

@ -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

Wyświetl plik

@ -0,0 +1,93 @@
#%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
import glob
import os
import shutil
from grass.pygrass.modules import Module
import grass.script as grass
import subprocess
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("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())

Wyświetl plik

@ -37,21 +37,18 @@ class TaskElevationMapGenerate(TaskView):
noise_filter_size = float(request.data.get('noise_filter_size', 2)) noise_filter_size = float(request.data.get('noise_filter_size', 2))
current_dir = os.path.dirname(os.path.abspath(__file__)) 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('interval', interval)
context.add_param('format', format) context.add_param('format', format)
context.add_param('noise_filter_size', noise_filter_size) context.add_param('noise_filter_size', noise_filter_size)
context.add_param('epsg', epsg) 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: if dtm != None:
context.add_param('dtm', '--dtm {}'.format(dtm)) context.add_param('dtm', dtm)
else:
context.add_param('dtm', '')
context.set_location(dsm) 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) return Response({'celery_task_id': celery_task_id}, status=status.HTTP_200_OK)
except GrassEngineException as e: except GrassEngineException as e:

Wyświetl plik

@ -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

Wyświetl plik

@ -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 cv2, math, argparse
import numpy as np import numpy as np
import rasterio as rio import rasterio as rio
from rasterio import warp, transform 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 # Open dsm
dsm = rio.open(args.dsm) dsm = rio.open(opts['dsm'])
# Read the tiff as an numpy masked array # Read the tiff as an numpy masked array
dsm_array = dsm.read(1, masked = True) dsm_array = dsm.read(1, masked = True)
# Create a kernel based on the parameter 'noise_filter_size' and the tiff resolution # 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 # Check if we want to use the dtm also
if args.dtm != None: if opts['dtm'] != '':
# Open the 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 that the dtm and dsm have the same bounds and resolution
assert_same_bounds_and_resolution(dsm, dtm) assert_same_bounds_and_resolution(dsm, dtm)
# Calculate the different between the dsm and dtm # Calculate the different between the dsm and dtm
@ -26,7 +87,7 @@ def main(args):
array = dsm_array array = dsm_array
# Calculate the ranges based on the parameter 'intervals' and the elevation 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 = [] features = []
@ -40,15 +101,32 @@ def main(args):
# Check if we found something # Check if we found something
if len(contours) > 0: if len(contours) > 0:
# Transform contours from pixels to coordinates # 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 # Build the MultiPolygon for based on the contours and their hierarchy
built_multi_polygon = LevelBuilder(bottom, top, mapped_contours, hierarchy[0]).build_multi_polygon() built_multi_polygon = LevelBuilder(bottom, top, mapped_contours, hierarchy[0]).build_multi_polygon()
features.append(built_multi_polygon) features.append(built_multi_polygon)
# Write the GeoJSON to a file # Write the GeoJSON to a file
dump = dumps(FeatureCollection(features)) dump = dumps(FeatureCollection(features))
with open(args.output, 'w+') as output: with open("output.json", 'w+') as output:
output.write(dump) 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): def get_kernel(noise_filter_size, dsm):
"""Generate a kernel for noise filtering. Will return none if the noise_filter_size isn't positive""" """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) multi_polygon = MultiPolygon(polygons)
return Feature(geometry = multi_polygon, properties = { 'bottom': int(self.bottom), 'top': int(self.top) }) 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') if __name__ == "__main__":
parser.add_argument('dsm', type = str, help = 'The path for the dsm file') opts, _ = grass.parser()
parser.add_argument('intervals', type = str, help = 'The intervals used to generate the diferent elevation levels') sys.exit(main())
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))

Wyświetl plik

@ -41,7 +41,7 @@ class TaskVolume(TaskView):
celery_task_id = execute_grass_script.delay(os.path.join( celery_task_id = execute_grass_script.delay(os.path.join(
os.path.dirname(os.path.abspath(__file__)), os.path.dirname(os.path.abspath(__file__)),
"calc_volume.grass" "calc_volume.py"
), context.serialize()).task_id ), context.serialize()).task_id
return Response({'celery_task_id': celery_task_id}, status=status.HTTP_200_OK) return Response({'celery_task_id': celery_task_id}, status=status.HTTP_200_OK)

Wyświetl plik

@ -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

Wyświetl plik

@ -0,0 +1,62 @@
#%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 and resolution to DSM
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())