[GRASS-user] patching vectors - some are not patched

[Ken Mankoff]

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 at 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 at 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 at 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 at routing res=30 e=e+1000 w=w-1000 n=n+1000 s=s-1000 --o
    r.water.outlet input=dir at 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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot_20181118_203058.png
Type: image/png
Size: 46122 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20181119/0e058440/attachment-0001.png>