Re: Splitting Up an Area

Moritz Lennert <[email protected]> Tue, 5 Dec 2017 16:05:22 +0100

Newsgroups	gmane.comp.gis.grass.user
Message-ID	<[email protected]>

On 05/12/17 01:09, Joshua Quesenberry wrote:
> All,
> 
> Any suggestions on how best I could split up a region into 40-60 acre 
> segments using linear features (ridges, streams, roads, trails, etc.)? 
> The watershed tool seems to split along ridges and ravines fairly well, 
> however, it seems I can only define the lower area limit not the upper..
> 

Are you speaking about linear features in images, or about vector maps 
of linear features ?

If this is for images, what you are speaking about are cutlines. AFAIK, 
we currently do not have an implementation of cutlines in GRASS GIS. 
Attached you can find a bash script I threw together trying to implement 
a cutlines algorithm using i.zc and r.cost+r.drain. It is in a rough 
state, without any comments, and it works only moderately well. It also 
doesn't scale well as i.zc puts everything into memory. I will probably 
have to work on a better version in a few months, but if you or someone 
else wants to do it before, I won't complain ;-)

You should be able to run it as is in the North Carolina demo dataset to 
see the results.

Moritz

_______________________________________________
grass-user mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/grass-user

cutlines.sh (application/x-shellscript, 3.3 KB)

RASTMAP=ortho_2001_t792_1m
ZC_THRESH=1
ZC_WIDTH=9
MINPIXELS=50
NUMBERLINES=10
EDGEWEIGHT=5
OUTPUTRASTER=cutlines10
OUTPUTVECTOR=cutlines10
MINTILESIZE=40000
MEMORY=10000


date
eval $(g.region rast=$RASTMAP -g)
i.zc $RASTMAP out=zctmp width=$ZC_WIDTH thresh=$ZC_THRESH --o --q
#i.group ${RASTMAP}_group in=${RASTMAP} --o --q
#i.superpixels.slic ${RASTMAP}_group out=tmp_superp minsize=${MINPIXELS} memory=${MEMORY} --o --q
#r.neighbors tmp_superp method=interspersion out=tmp_interspers --o --q

date
NSSTEP=$(python -c "print 1.0*($n-$s-$nsres)/$NUMBERLINES")
python -c "for i in range(0,($NUMBERLINES+1)): print str($w+0.5)+','+ str(($n-i*$NSSTEP)-$nsres/2.0)" > hstartpoints
python -c "for i in range(0,($NUMBERLINES+1)): print str($e-0.5)+','+ str(($n-i*$NSSTEP)-$nsres/2.0)" > hstoppoints
v.in.ascii in=hstartpoints sep=comma out=hstart --o --q
v.in.ascii in=hstoppoints sep=comma out=hstop --o --q
v.distance from=hstart to=hstop out=hborderlines_tmp --o --q
HALFNSSTEP=$(python -c "print $NSSTEP / 2.0")
v.transform hborderlines_tmp yshift=-$HALFNSSTEP out=hborderlines --o --q
v.to.rast hborderlines use=val out=hborderlines type=line --o --q
r.mapcalc "hbase = if(isnull(hborderlines), if(zctmp==0, $EDGEWEIGHT, 1), 9999)" --o --q
#r.mapcalc "hbase = if(isnull(hborderlines), if(tmp_interspers==1, $EDGEWEIGHT, 1), 9999)" --o --q
r.cost hbase startp=hstart stopp=hstop outp=hcumcost --o --q
v.transform hstop out=hdrainstart xshift=-1 --o --q
r.drain in=hcumcost startp=hdrainstart out=hlines --o --q

date
EWSTEP=$(python -c "print 1.0*($e-$w-$ewres)/$NUMBERLINES")
python -c "for i in range(0,($NUMBERLINES+1)): print str(($e-i*$EWSTEP)-$ewres/2.0)+','+str($n-0.5)" > vstartpoints
python -c "for i in range(0,($NUMBERLINES+1)): print str(($e-i*$EWSTEP)-$ewres/2.0)+','+str($s+0.5)" > vstoppoints
v.in.ascii in=vstartpoints sep=comma out=vstart --o --q
v.in.ascii in=vstoppoints sep=comma out=vstop --o --q
v.distance from=vstart to=vstop out=vborderlines_tmp --o --q
HALFEWSTEP=$(python -c "print $EWSTEP / 2.0")
v.transform vborderlines_tmp xshift=+$HALFEWSTEP out=vborderlines --o --q
v.to.rast vborderlines use=val out=vborderlines type=line --o --q
r.mapcalc "vbase = if(isnull(vborderlines), if(zctmp==0, $EDGEWEIGHT, 1), 9999)" --o --q
#r.mapcalc "vbase = if(isnull(vborderlines), if(tmp_interspers==1, $EDGEWEIGHT, 1), 9999)" --o --q
r.cost vbase startp=vstart stopp=vstop outp=vcumcost --o --q
v.transform vstop out=vdrainstart yshift=+1 --o --q
r.drain in=vcumcost startp=vdrainstart out=vlines --o --q

date
r.patch hlines,vlines out=$OUTPUTRASTER --o --q
r.to.vect $OUTPUTRASTER out=$OUTPUTVECTOR type=line --o --q

v.in.region out=region type=line --o --q
v.patch cutlines,region out=tmp1 --o --q
snapdist=$(python -c "print $nsres*2")
echo "snapdist = $snapdist"
v.clean tmp1 tool=snap,break,rmdupl out=tmp2 thresh=$snapdist --o --q
v.type tmp2 out=tmp3 --o --q
v.centroids tmp3 out=tmp4 --o --q
v.clean tmp4 tool=rmarea thresh=$MINTILESIZE out=cutpolygons --o --q

g.remove type=rast,vect name=base,hstart,vstart,hstop,vstop,hdrainstart,vdrainstart,hcumcost,vcumcost,hlines,vlines,zctmp,hbase,hborderlines,vbase,vborderlines,hborderlines_tmp,vborderlines,region,vborderlines_tmp,hborderlines,tmp1,tmp2,tmp3,tmp4,tmp_interspers,tmp_superp -f --q
rm hstartpoints hstoppoints vstartpoints vstoppoints


date