Hi GRASS List,
[Sorry if you see this 2x. Sent to grass-dev list by mistake].
I'm having trouble getting v.patch to patch some vectors. I've tested this on
Linux grass 7.4.0 (default Ubuntu 18.04 version) and SVN trunk from this
evening.
I'm trying to "accumulate" vectors, incrementally, via v.patch with "-a" flag.
In general things appear to work fine, but every once in a while on vector is
simply not included. I am having trouble figuring out why. I can send an MWE
with ~30 lines of code and a few small input rasters to re-create the issue if
that would help. The MWE runs in ~30 seconds. For now I begin with just this
email.
I've run r.watershed. Wherever the direction raster is negative (outlets), I
run r.water.outlet to find the upstream catchment. I use the following
algorithm to accumulate the catchments into one vector:
For each outlet:
1) Set region to outlet +- 1 km
2) Run r.water.outlet
3) Check if any edge cells are non-null (catchments meets edge of region). If
so, expand region and repeat until check passes
4) Set output raster values to a unique (incrementing) number
5) Convert to area vector "b"
6) Patch, with "v.patch -a input=b output=basin --o"
where "basin" was created initially empty with "v.edit tool=create map=basin".
When I do all of this, things seem OK but I cannot display the vector until I
clean it:
v.build map=basin option=build
v.extract input=basin output=basin_clean type=area --o
g.rename vector=basin_clean,basin --o
Then I can display it. When I do, some catchments are missing. The missing
catchments are valid through step (5) above, but are not patched with step (6).
Most are, some are not, and there doesn't appear to be anything unique or
special about the problematic ones.
Attached is a figure showing basins ("d.vect -c basins"). The white vertical
chunk in the middle is two missing (unpatched) vectors. The vectors are
produced and displayed after step 5.
Oddly, I notice now while looking at this, that the burgundy colors bordering
the white (1 cell near the outlet, 3 more on the right lower, and 2 more on the
left lower again) are all assigned the same value, implying that is one
contiguous outlet. Presumably this bug is the cause or caused by the v.patch
issue.
The actual code follows (I can provide the 'dir' input upon request). Any
advice on fixing this bug, or general code comments will be much appreciated.
Thanks,
-k.
r.mapcalc "outlets = if(dir@routing < 0, 1, null())" --o
r.out.xyz input=outlets separator=comma output=- | cut -d"," -f1,2 > outlets.csv
count=0
v.edit tool=create map=basin --o --quiet
coords=228905,-669905 # DEBUG
for coords in $(cat outlets.csv); do
ew=$(echo ${coords} | cut -d, -f1)
ns=$(echo ${coords} | cut -d, -f2)
# make a 1x1 km region around the starting grid cell
g.region align=dir@routing res=30 e=$(( $ew + 1000 )) w=$(( $ew - 1000 ))
n=$(( $ns + 1000 )) s=$(( $ns - 1000 ))
# define the basin that feeds this outlet
r.water.outlet input=dir@routing output=basin coordinates=${coords} --o --q
# check if any edge cells for this sub-region have data
# If so, we hit the edge and 1x1 km was too small
# Keep enlarging area until outlet doesn't reach edge.
r.mapcalc "mask = basin * if((row() == 1) | (row() == nrows()) | (col() == 1)
| (col() == ncols()), 1, null())" --o --q
edge_cells=$(r.info mask | grep -o "max =...." | cut -d" " -f3)
# While basin is in contact w/ the edge, expand by X km in each direction
while [[ ${edge_cells} != "NUL" ]]; do
# echo "expanding ${coords}"
g.region align=dir@routing res=30 e=e+1000 w=w-1000 n=n+1000 s=s-1000 --o
r.water.outlet input=dir@routing output=basin coordinates=${coords} --o --q
# check again
r.mapcalc "mask = basin * if((row() == 1) | (row() == nrows()) | (col() ==
1) | (col() == ncols()), 1, null())" --o --q
edge_cells=$(r.info mask | grep -o "max =...." | cut -d" " -f3)
done
r.mapcalc "basin = if(basin, ${count}, null())" --o
count=$(( ${count} + 1 ))
r.to.vect -v input=basin output=b type=area --o
# clean up new vector
v.build map=b option=build &>> /dev/null
v.extract input=b output=b_clean type=area --o
g.rename vector=b_clean,b --o &>> /dev/null
# clean up accumulation vector
v.build map=basin option=build &>> /dev/null
v.extract input=basin output=basin_clean type=area --o
g.rename vector=basin_clean,basin --o
# patch new vector to existing vectors
v.patch -a -b -n input=b output=basin --o
done
# clean up again
v.build map=basin option=build
v.extract input=basin output=basin_clean type=area --o
g.rename vector=basin_clean,basin --o
# visual debug
d.mon start=wx0
d.erase
d.vect -c basin
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user
