<div dir="ltr">Hi,<div><br></div><div>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.</div><div>```</div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div># %module</div><div># % description: Calculate volume of area and prints the volume to stdout</div><div># %end</div><div># %option</div><div># % key: area_file</div><div># % type: string</div><div># % required: yes</div><div># % multiple: no</div><div># % description: Geospatial file containing the area to measure</div><div># %end</div><div># %option</div><div># % key: points_file</div><div># % type: string</div><div># % required: yes</div><div># % multiple: no</div><div># % description: Geospatial file containing the points defining the area</div><div># %end</div><div># %option</div><div># % key: dsm_file</div><div># % type: string</div><div># % required: yes</div><div># % multiple: no</div><div># % description: GeoTIFF DEM containing the surface</div><div># %end</div><div><br></div><div>import sys</div><div>from grass.pygrass.modules import Module</div><div>import grass.script as grass</div><div><br></div><div><br></div><div>def main():</div><div><br></div><div> Module(</div><div> "v.import",</div><div> input=opts['area_file'],</div><div> output="polygon_area",</div><div> overwrite=True</div><div> )</div><div> Module(</div><div> "v.import",</div><div> input=opts['points_file'],</div><div> output="polygon_points",</div><div> overwrite=True</div><div> )</div><div> Module(</div><div> "v.buffer",</div><div> input="polygon_area",</div><div> s=True,</div><div> type="area",</div><div> output="region",</div><div> distance=1,</div><div> minordistance=1,</div><div> overwrite=True</div><div> )</div><div> Module("r.external", input=opts['dsm_file'], output="dsm", overwrite=True)</div><div> # Set Grass region and resolution to DSM</div><div> Module("g.region", raster="dsm")</div><div> # Set Grass region to vector bbox</div><div> Module("g.region", vector="region")</div><div> # Create a mask to speed up computation</div><div> Module("r.mask", vector="region")</div><div> # Transfer dsm raster data to vector</div><div> Module("v.what.rast", map="polygon_points", raster="dsm", column="height")</div><div> # Decimate DSM and generate interpolation of new terrain</div><div> Module("v.surf.rst", input="polygon_points", zcolumn="height",</div><div> elevation="dsm_below_pile", overwrite=True)</div><div> # Compute difference between dsm and new dsm</div><div> Module("r.mapcalc",</div><div> expression='pile_height_above_dsm=dsm-dsm_below_pile', overwrite=True)</div><div> # Set region to polygon area to calculate volume</div><div> Module("g.region", vector="polygon_area")</div><div> # Volume output from difference</div><div> Module("r.volume", input="pile_height_above_dsm", f=True)</div><div><br></div><div> return 0</div><div><br></div><div><br></div><div>if __name__ == "__main__":</div><div> opts, _ = grass.parser()</div><div> sys.exit(main())</div></blockquote>```<div><br><div><div>The area file I am providing contains the below contents(area_file.geojson):</div><div><br></div><div>```</div></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div><div>{"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]]]}}</div></div></blockquote><div><div>```</div><div><br></div><div><div>The points file I am providing contains the below contents(points_file.geojson):</div><div><br></div><div>```</div></div></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div><div><div>{"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]}}]}</div></div></div></blockquote><div><div><div>```<br></div></div><div><br></div><div>When I run the below command,</div><div><br></div><div>```</div></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div><div><b>/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78 -c \</b></div><div><b>/home/ec2-user/volume_calc/media/tour_id/files/dsm_volume.tif location \</b></div><div><b>--exec python3 /home/ec2-user/volume_calc/grass_calc_vol.py \</b></div><div><b>area_file=/home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/area_file.geojson \</b></div><div><b>points_file=/home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/points_file.geojson \</b></div><div><b>dsm_file=/home/ec2-user/volume_calc/media/tour_id/files/dsm_volume.tif</b></div></div></blockquote><div><div>```</div><div><br></div><div>I get below output.</div><div>```</div></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div><div>(.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</div></div><div><div>Starting GRASS GIS...</div></div><div><div>Creating new GRASS GIS location <location>...</div></div><div><div>Cleaning up temporary files...</div></div><div><div>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> ...</div></div><div><div>Check if OGR layer <OGRGeoJSON> contains polygons...</div></div><div><div> 100%</div></div><div><div>Creating attribute table for layer <OGRGeoJSON>...</div></div><div><div>Importing 1 features (OGR layer <OGRGeoJSON>)...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Registering primitives...</div></div><div><div>-----------------------------------------------------</div></div><div><div>Cleaning polygons</div></div><div><div>-----------------------------------------------------</div></div><div><div>Breaking polygons...</div></div><div><div>Breaking polygons (pass 1: select break points)...</div></div><div><div> 100%</div></div><div><div>Breaking polygons (pass 2: break at selected points)...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Removing duplicates...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Breaking boundaries...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Removing duplicates...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Cleaning boundaries at nodes...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Merging boundaries...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Removing dangles...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Building areas...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Removing bridges...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Registering primitives...</div></div><div><div>Building areas...</div></div><div><div> 100%</div></div><div><div>Attaching islands...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Finding centroids for OGR layer <OGRGeoJSON>...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Writing centroids...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>1 input polygons</div></div><div><div>Total area: 4106.45 (1 areas)</div></div><div><div>-----------------------------------------------------</div></div><div><div>Copying features...</div></div><div><div> 100%</div></div><div><div>Building topology for vector map <polygon_area@PERMANENT>...</div></div><div><div>Registering primitives...</div></div><div><div>Building areas...</div></div><div><div> 100%</div></div><div><div>Attaching islands...</div></div><div><div> 100%</div></div><div><div>Attaching centroids...</div></div><div><div> 100%</div></div><div><div>Input</div></div><div><div></home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/area_file.geojson></div></div><div><div>successfully imported without reprojection</div></div><div><div>Check if OGR layer <OGRGeoJSON> contains polygons...</div></div><div><div> 100%</div></div><div><div>Creating attribute table for layer <OGRGeoJSON>...</div></div><div><div>Importing 5 features (OGR layer <OGRGeoJSON>)...</div></div><div><div> 100%</div></div><div><div>-----------------------------------------------------</div></div><div><div>Building topology for vector map <polygon_points@PERMANENT>...</div></div><div><div>Registering primitives...</div></div><div><div>Input</div></div><div><div></home/ec2-user/volume_calc/media/tour_id/tmpo99fe5bi_grass_engine/points_file.geojson></div></div><div><div>successfully imported without reprojection</div></div><div><div>Note: In latitude-longitude coordinate system specify distances in degree</div></div><div><div>unit</div></div><div><div>WARNING: Option 'minordistance' is not available with GEOS buffering</div></div><div><div>Buffering areas...</div></div><div><div> 100%</div></div><div><div>Cleaning buffers...</div></div><div><div>Building parts of topology...</div></div><div><div>Building topology for vector map <region@PERMANENT>...</div></div><div><div>Registering primitives...</div></div><div><div>Snapping boundaries...</div></div><div><div>Reading features...</div></div><div><div>Snap vertices Pass 1: select points</div></div><div><div> 100%</div></div><div><div>Snap vertices Pass 2: assign anchor vertices</div></div><div><div> 100%</div></div><div><div>Snap vertices Pass 3: snap to assigned points</div></div><div><div> 100%</div></div><div><div>Breaking polygons...</div></div><div><div>Breaking polygons (pass 1: select break points)...</div></div><div><div> 100%</div></div><div><div>Breaking polygons (pass 2: break at selected points)...</div></div><div><div> 100%</div></div><div><div>Removing duplicates...</div></div><div><div> 100%</div></div><div><div>Breaking boundaries...</div></div><div><div> 100%</div></div><div><div>Removing duplicates...</div></div><div><div> 100%</div></div><div><div>Cleaning boundaries at nodes</div></div><div><div> 100%</div></div><div><div>Building topology for vector map <region@PERMANENT>...</div></div><div><div>Building areas...</div></div><div><div> 100%</div></div><div><div>Removing dangles...</div></div><div><div> 100%</div></div><div><div>Removing bridges...</div></div><div><div> 100%</div></div><div><div>Attaching islands...</div></div><div><div>Building topology for vector map <region@PERMANENT>...</div></div><div><div>Attaching islands...</div></div><div><div> 100%</div></div><div><div>Calculating centroids for all areas...</div></div><div><div> 100%</div></div><div><div>Generating list of boundaries to be deleted...</div></div><div><div> 100%</div></div><div><div>Deleting boundaries...</div></div><div><div> 100%</div></div><div><div>Calculating centroids for areas...</div></div><div><div> 100%</div></div><div><div>Building topology for vector map <region@PERMANENT>...</div></div><div><div>Registering primitives...</div></div><div><div>Building areas...</div></div><div><div> 100%</div></div><div><div>Attaching islands...</div></div><div><div> 100%</div></div><div><div>Attaching centroids...</div></div><div><div> 100%</div></div><div><div>WARNING: Unable to open datum table file</div></div><div><div> </home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/datum.table></div></div><div><div>WARNING: Datum <WGS_1984> not recognised by GRASS and no parameters found</div></div><div><div>WARNING: Unable to open ellipsoid table file</div></div><div><div> </home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/ellipse.table></div></div><div><div>Projection of input dataset and current location appear to match</div></div><div><div>Importing band 1 of 1...</div></div><div><div>Link to raster map <dsm> created</div></div><div><div>Reading areas...</div></div><div><div> 100%</div></div><div><div>Writing raster map...</div></div><div><div> 100%</div></div><div><div>Reading areas...</div></div><div><div> 100%</div></div><div><div>Writing raster map...</div></div><div><div> 100%</div></div><div><div>Reading areas...</div></div><div><div> 100%</div></div><div><div>Writing raster map...</div></div><div><div> 100%</div></div><div><div>Reading areas...</div></div><div><div> 100%</div></div><div><div>Writing raster map...</div></div><div><div>^CTraceback (most recent call last):</div></div><div><div> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/scripts/r.mask", line 199, in <module></div></div><div><div> main()</div></div><div><div> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/scripts/r.mask", line 179, in main</div></div><div><div> type='area', cats=cats, where=where, env=env)</div></div><div><div> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/script/core.py", line 440, in run_command</div></div><div><div> returncode = ps.wait()</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1019, in wait</div></div><div><div> return self._wait(timeout=timeout)</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1653, in _wait</div></div><div><div> (pid, sts) = self._try_wait(0)</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1611, in _try_wait</div></div><div><div> (pid, sts) = os.waitpid(self.pid, wait_flags)</div></div><div><div>KeyboardInterrupt</div></div><div><div>Traceback (most recent call last):</div></div><div><div> File "/home/ec2-user/volume_calc/grass_calc_vol.py", line 80, in <module></div></div><div><div> sys.exit(main())</div></div><div><div> File "/home/ec2-user/volume_calc/grass_calc_vol.py", line 61, in main</div></div><div><div> Module("r.mask", vector="region")</div></div><div><div> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/pygrass/modules/interface/module.py", line 602, in __init__</div></div><div><div> self.__call__(*args, **kargs)</div></div><div><div> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/pygrass/modules/interface/module.py", line 660, in __call__</div></div><div><div> return self.run()</div></div><div><div> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/pygrass/modules/interface/module.py", line 783, in run</div></div><div><div> self.wait()</div></div><div><div> File "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/etc/python/grass/pygrass/modules/interface/module.py", line 796, in wait</div></div><div><div> stdout, stderr = self._popen.communicate(input=self.stdin)</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 956, in communicate</div></div><div><div> self.wait()</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1019, in wait</div></div><div><div> return self._wait(timeout=timeout)</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1653, in _wait</div></div><div><div> (pid, sts) = self._try_wait(0)</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1611, in _try_wait</div></div><div><div> (pid, sts) = os.waitpid(self.pid, wait_flags)</div></div><div><div>KeyboardInterrupt</div></div><div><div>Traceback (most recent call last):</div></div><div><div> File "/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78", line 2348, in <module></div></div><div><div> main()</div></div><div><div> File "/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78", line 2286, in main</div></div><div><div> returncode = run_batch_job(batch_job)</div></div><div><div> File "/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/grass78", line 1635, in run_batch_job</div></div><div><div> returncode = proc.wait()</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1019, in wait</div></div><div><div> return self._wait(timeout=timeout)</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1653, in _wait</div></div><div><div> (pid, sts) = self._try_wait(0)</div></div><div><div> File "/usr/lib64/python3.7/subprocess.py", line 1611, in _try_wait</div></div><div><div> (pid, sts) = os.waitpid(self.pid, wait_flags)</div></div><div><div>KeyboardInterrupt</div></div></blockquote><div><div>```</div><div><br></div><div>I had to hit <b>CTRL + C </b>to stop the infinite loop. </div><div><br></div><div>Few lines from my grass78 bin file inside <b>/home/ec2-user/grass/bin.x86_64-pc-linux-gnu/</b></div><div><br></div><div>```</div></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div>GISBASE = os.path.normpath("/home/ec2-user/grass/dist.x86_64-pc-linux-gnu")</div><div>os.environ['GISBASE'] = GISBASE</div><div>CMD_NAME = "grass78"</div><div>GRASS_VERSION = "7.8.6dev"</div><div>LD_LIBRARY_PATH_VAR = 'LD_LIBRARY_PATH'</div><div>CONFIG_PROJSHARE = os.environ.get('GRASS_PROJSHARE', "/usr/share/proj")</div></blockquote><div><div>```</div><div>Environment variables that I am setting:</div><div><br></div><div>```</div> env["LC_ALL"] = "C.UTF-8"<br> env["GISBASE"] = "/home/ec2-user/grass/dist.x86_64-pc-linux-gnu"<br> env["GISDBASE"] = self.tmpdir # in the above code eg.
<b>tmpo99fe5bi_grass_engine</b><br> env['PATH'] = "$PATH:/usr/bin:$GISBASE/bin:$GISBASE/scripts"<br> env['LD_LIBRARY_PATH'] = "$LD_LIBRARY_PATH:$GISBASE/lib"<br> env['GISRC'] = "$HOME/.grass7"<br> env['PYTHONPATH'] = "$HOME/grass/dist.x86_64-pc-linux-gnu/etc/python"<div>```</div><div><br></div><div>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</div></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div><div><b><br></b></div><div><b>/home/ec2-user/grass/dist.x86_64-pc-linux-gnu/bin/</b></div></div></blockquote><div><br></div>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)?<div>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.<br><div>Please let me know if you need more information.</div><div><br></div><div>--</div><div><div dir="ltr" data-smartmail="gmail_signature"><div dir="ltr"><div>Regards,</div><div>Ayush Chauhan</div><div><br></div></div></div></div></div></div></div>