Hi,
I am using the below code in my grass script to calculate volume from the DSM file input into the script and the polygon points on the DSM.
> # %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():
>
> 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 and resolution to DSM
> Module("g.region", raster="dsm")
> # 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())
The area file I am providing contains the below contents(area_file.geojson):
> {"type": "Feature", "properties": {}, "geometry": {"type": "Polygon", "coordinates": [[[73.31231666938123, 26.22297239890824], [73.31241572405983, 26.222379996744422], [73.3134282830031, 26.222577464467804], [73.31281194277611, 26.222906576595037], [73.31231666938123, 26.22297239890824]]]}}
The points file I am providing contains the below contents(points_file.geojson):
> {"type": "FeatureCollection", "features": [{"type": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [73.31231666938123, 26.22297239890824]}}, {"type": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [73.31241572405983, 26.222379996744422]}}, {"type": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [73.3134282830031, 26.222577464467804]}}, {"type": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [73.31281194277611, 26.222906576595037]}}, {"type": "Feature", "properties": {}, "geometry": {"type": "Point", "coordinates": [73.31231666938123, 26.22297239890824]}}]}
When I run the below command,
> **/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78 -c \**
> **/home/ec2-user/volume_calc/media/tour_id/files/dsm_volume.tif location \**
> **--exec python3 /home/ec2-user/volume_calc/grass_calc_vol.py \**
> **area_file=/home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/area_file.geojson \**
> **points_file=/home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/points_file.geojson \**
> **dsm_file=/home/ec2-user/volume_calc/media/tour_id/files/dsm_volume.tif**
I get below output.
> (.venv) [ec2-user@ip-10-0-3-14 volume_calc]$ /home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78 -c /home/ec2-user/volume_calc/media/tour_id/files/dsm_volume.tif location --exec python3 /home/ec2-user/volume_calc/grass_calc_vol.py area_file=/home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/area_file.geojson points_file=/home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/points_file.geojson dsm_file=/home/ec2-user/volume_calc/media/tour_id/files/dsm_volume.tif
>
> Starting GRASS GIS...
>
> Creating new GRASS GIS location <location>...
>
> Cleaning up temporary files...
>
> Executing <python3 /home/ec2-user/volume_calc/grass_calc_vol.py area_file=/home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/area_file.geojson points_file=/home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/points_file.geojson dsm_file=/home/ec2-user/volume_calc/media/tour_id/files/dsm_volume.tif> ...
>
> Check if OGR layer <OGRGeoJSON> contains polygons...
>
> 100%
>
> Creating attribute table for layer <OGRGeoJSON>...
>
> Importing 1 features (OGR layer <OGRGeoJSON>)...
>
> 100%
>
> -----------------------------------------------------
>
> Registering primitives...
>
> -----------------------------------------------------
>
> Cleaning polygons
>
> -----------------------------------------------------
>
> Breaking polygons...
>
> Breaking polygons (pass 1: select break points)...
>
> 100%
>
> Breaking polygons (pass 2: break at selected points)...
>
> 100%
>
> -----------------------------------------------------
>
> Removing duplicates...
>
> 100%
>
> -----------------------------------------------------
>
> Breaking boundaries...
>
> 100%
>
> -----------------------------------------------------
>
> Removing duplicates...
>
> 100%
>
> -----------------------------------------------------
>
> Cleaning boundaries at nodes...
>
> 100%
>
> -----------------------------------------------------
>
> Merging boundaries...
>
> 100%
>
> -----------------------------------------------------
>
> Removing dangles...
>
> 100%
>
> -----------------------------------------------------
>
> Building areas...
>
> 100%
>
> -----------------------------------------------------
>
> Removing bridges...
>
> 100%
>
> -----------------------------------------------------
>
> Registering primitives...
>
> Building areas...
>
> 100%
>
> Attaching islands...
>
> 100%
>
> -----------------------------------------------------
>
> Finding centroids for OGR layer <OGRGeoJSON>...
>
> 100%
>
> -----------------------------------------------------
>
> Writing centroids...
>
> 100%
>
> -----------------------------------------------------
>
> 1 input polygons
>
> Total area: 4106.45 (1 areas)
>
> -----------------------------------------------------
>
> Copying features...
>
> 100%
>
> Building topology for vector map <polygon_area@PERMANENT>...
>
> Registering primitives...
>
> Building areas...
>
> 100%
>
> Attaching islands...
>
> 100%
>
> Attaching centroids...
>
> 100%
>
> Input
>
> </home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/area_file.geojson>
>
> successfully imported without reprojection
>
> Check if OGR layer <OGRGeoJSON> contains polygons...
>
> 100%
>
> Creating attribute table for layer <OGRGeoJSON>...
>
> Importing 5 features (OGR layer <OGRGeoJSON>)...
>
> 100%
>
> -----------------------------------------------------
>
> Building topology for vector map <polygon_points@PERMANENT>...
>
> Registering primitives...
>
> Input
>
> </home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/points_file.geojson>
>
> successfully imported without reprojection
>
> Note: In latitude-longitude coordinate system specify distances in degree
>
> unit
>
> WARNING: Option 'minordistance' is not available with GEOS buffering
>
> Buffering areas...
>
> 100%
>
> Cleaning buffers...
>
> Building parts of topology...
>
> Building topology for vector map <region@PERMANENT>...
>
> Registering primitives...
>
> Snapping boundaries...
>
> Reading features...
>
> Snap vertices Pass 1: select points
>
> 100%
>
> Snap vertices Pass 2: assign anchor vertices
>
> 100%
>
> Snap vertices Pass 3: snap to assigned points
>
> 100%
>
> Breaking polygons...
>
> Breaking polygons (pass 1: select break points)...
>
> 100%
>
> Breaking polygons (pass 2: break at selected points)...
>
> 100%
>
> Removing duplicates...
>
> 100%
>
> Breaking boundaries...
>
> 100%
>
> Removing duplicates...
>
> 100%
>
> Cleaning boundaries at nodes
>
> 100%
>
> Building topology for vector map <region@PERMANENT>...
>
> Building areas...
>
> 100%
>
> Removing dangles...
>
> 100%
>
> Removing bridges...
>
> 100%
>
> Attaching islands...
>
> Building topology for vector map <region@PERMANENT>...
>
> Attaching islands...
>
> 100%
>
> Calculating centroids for all areas...
>
> 100%
>
> Generating list of boundaries to be deleted...
>
> 100%
>
> Deleting boundaries...
>
> 100%
>
> Calculating centroids for areas...
>
> 100%
>
> Building topology for vector map <region@PERMANENT>...
>
> Registering primitives...
>
> Building areas...
>
> 100%
>
> Attaching islands...
>
> 100%
>
> Attaching centroids...
>
> 100%
>
> WARNING: Unable to open datum table file
>
> </home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/datum.table>
>
> WARNING: Datum <WGS_1984> not recognised by GRASS and no parameters found
>
> WARNING: Unable to open ellipsoid table file
>
> </home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/ellipse.table>
>
> Projection of input dataset and current location appear to match
>
> Importing band 1 of 1...
>
> Link to raster map <dsm> created
>
> Reading areas...
>
> 100%
>
> Writing raster map...
>
> 100%
>
> Reading areas...
>
> 100%
>
> Writing raster map...
>
> 100%
>
> Reading areas...
>
> 100%
>
> Writing raster map...
>
> 100%
>
> Reading areas...
>
> 100%
>
> Writing raster map...
>
> ^CTraceback (most recent call last):
>
> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/scripts/r.mask", line 199, in <module>
>
> main()
>
> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/scripts/r.mask", line 179, in main
>
> type='area', cats=cats, where=where, env=env)
>
> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/script/core.py", line 440, in run_command
>
> returncode = ps.wait()
>
> File "/usr/lib64/python3.7/subprocess.py", line 1019, in wait
>
> return self._wait(timeout=timeout)
>
> File "/usr/lib64/python3.7/subprocess.py", line 1653, in _wait
>
> (pid, sts) = self._try_wait(0)
>
> File "/usr/lib64/python3.7/subprocess.py", line 1611, in _try_wait
>
> (pid, sts) = os.waitpid(self.pid, wait_flags)
>
> KeyboardInterrupt
>
> Traceback (most recent call last):
>
> File "/home/ec2-user/volume_calc/grass_calc_vol.py", line 80, in <module>
>
> sys.exit(main())
>
> File "/home/ec2-user/volume_calc/grass_calc_vol.py", line 61, in main
>
> Module("r.mask", vector="region")
>
> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/pygrass/modules/interface/module.py", line 602, in __init__
>
> self.__call__(*args, **kargs)
>
> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/pygrass/modules/interface/module.py", line 660, in __call__
>
> return self.run()
>
> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/pygrass/modules/interface/module.py", line 783, in run
>
> self.wait()
>
> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/pygrass/modules/interface/module.py", line 796, in wait
>
> stdout, stderr = self._popen.communicate(input=self.stdin)
>
> File "/usr/lib64/python3.7/subprocess.py", line 956, in communicate
>
> self.wait()
>
> File "/usr/lib64/python3.7/subprocess.py", line 1019, in wait
>
> return self._wait(timeout=timeout)
>
> File "/usr/lib64/python3.7/subprocess.py", line 1653, in _wait
>
> (pid, sts) = self._try_wait(0)
>
> File "/usr/lib64/python3.7/subprocess.py", line 1611, in _try_wait
>
> (pid, sts) = os.waitpid(self.pid, wait_flags)
>
> KeyboardInterrupt
>
> Traceback (most recent call last):
>
> File "/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78", line 2348, in <module>
>
> main()
>
> File "/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78", line 2286, in main
>
> returncode = run_batch_job(batch_job)
>
> File "/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78", line 1635, in run_batch_job
>
> returncode = proc.wait()
>
> File "/usr/lib64/python3.7/subprocess.py", line 1019, in wait
>
> return self._wait(timeout=timeout)
>
> File "/usr/lib64/python3.7/subprocess.py", line 1653, in _wait
>
> (pid, sts) = self._try_wait(0)
>
> File "/usr/lib64/python3.7/subprocess.py", line 1611, in _try_wait
>
> (pid, sts) = os.waitpid(self.pid, wait_flags)
>
> KeyboardInterrupt
I had to hit CTRL + C to stop the infinite loop.
Few lines from my grass78 bin file inside /home/ec2-user/grass/bin.x86_64-pc-linux-gnu/
> GISBASE = os.path.normpath("/home/ec2-user/grass/dist.x86_64-pc-linux-gnu")
> os.environ['GISBASE'] = GISBASE
> CMD_NAME = "grass78"
> GRASS_VERSION = "7.8.6dev"
> LD_LIBRARY_PATH_VAR = 'LD_LIBRARY_PATH'
> CONFIG_PROJSHARE = os.environ.get('GRASS_PROJSHARE', "/usr/share/proj")
Environment variables that I am setting:
env["LC_ALL"] = "C.UTF-8"
env["GISBASE"] = "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu"
env["GISDBASE"] = self.tmpdir # in the above code eg. **tmpo99fe5bi_grass_engine**
env['PATH'] = "$PATH:/usr/bin:$GISBASE/bin:$GISBASE/scripts"
env['LD_LIBRARY_PATH'] = "$LD_LIBRARY_PATH:$GISBASE/lib"
env['GISRC'] = "$HOME/.grass7"
env['PYTHONPATH'] = "$HOME/grass/dist.x86_64-pc-linux-gnu/etc/python"
Also, I had copied the r.external file from version grass-6.4.4 that was installed automatically when I ran sudo yum install grass. Because it was not present in the bin folder of this location
/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/bin/
Is there a way to get the file from outside, thus removing the grass-6.4.4 version and using the grass-7.8.5 version r.external file (can that be a cause of infinite loop)?Any suggestions as to what is going wrong? or anything wrong in the Grass SCRIPT mentioned above, that could potentially lead to an infinite loop.
Please let me know if you need more information.
–
Regards,
Ayush Chauhan