[GRASS-user] GRASS SCRIPT Execute Error
Ayush Chauhan
ayush at huviair.com
Sun Jul 18 23:16:06 PDT 2021
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 at 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 at 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 at 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 at 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 at PERMANENT>...
Building areas...
100%
Removing dangles...
100%
Removing bridges...
100%
Attaching islands...
Building topology for vector map <region at 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 at 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20210719/b09cf959/attachment-0001.html>
More information about the grass-user
mailing list