[GRASS-user] r.contour fails to close a contour at the region's border

[Markus Metz]

On Thu, Jan 9, 2014 at 9:57 AM, Markus Metz
<markus.metz.giswork at gmail.com> wrote:
> On Thu, Jan 9, 2014 at 1:46 AM, John Ciolek <jciolek at alphatrac.com> wrote:
>>
>> On Jan 8, 2014, at 12:35 PM, Markus Metz wrote:
>>
>>> On Thu, Dec 12, 2013 at 11:13 PM, John Ciolek <jciolek at alphatrac.com> wrote:
>>>> I have automatically generated raster data from which I create a contour
>>>> using r.contour.  Sometimes the "feature" to be contoured is not completely
>>>> contained within the specified region.  When this happens, r.contour
>>>> generates an open contour - basically the beginning line segment point does
>>>> not equal the ending line segment point.
>>>>
>>>> Is there a way to get r.contour to close the contour at the region's border?
>>>
>>> No. This is a feature of r.contour because it is impossible to tell
>>> where a contour would continue if it hits the region's border.
>>>>
>>>> If not, is there an easy way to close the contour?
>>>
>>> Why do you want to close the contour?
>>
>> I am trying to create areas from the contours, where the areas clip to the set region.
>>
>>>
>>> If you want to convert the contour lines to areas, you might instead
>>> modify the raster with r.mapcalc, e.g. for contours at every 100 unit
>>> step
>>>
>>> r.mapcalc "surface.classified = int(surface / 100) * 100"
>>>
>>> but here a value of 99 would be converted to 0, thus as an alternative
>>>
>>> r.mapcalc "surface.classified = int((surface + 50) / 100) * 100"
>>>
>>> which would convert e.g. all values >= 50 and < 150 to 100.
>>>
>>> After that the classified raster could be converted to vector areas
>>> with r.to.vect type=area
>>>
>>> HTH,
>>>
>>> Markus M
>>
>>
>> Wouldn't this create a very blocky contour?
>
> The contours are as blocky as the current region settings, also with
> r.contour.

Sorry, you are right, the contours generated by r.contour are not
blocky. You can convert the contour lines to areas with something like

--->
DEM="elev_state_500m"
step=100

halfstep=`echo $step | awk '{printf "%.15g\n", $1 / 2.0}'`

g.region -p rast=$DEM

r.contour in=$DEM out=${DEM}_contour_3d step=$step
r.mapcalc "${DEM}_bin = if(isnull($DEM), null(), 1)"
# double the resolution
eval `g.region -g`
rows2=`echo $rows | awk '{printf "%d\n", $1 * 2}'`
cols2=`echo $cols | awk '{printf "%d\n", $1 * 2}'`
eval `g.region -g rows=$rows2 cols=$cols2`
# grow region by one cell
g.region -g n=n+$nsres s=s-$nsres e=e+$ewres w=w-$ewres
r.mapcalc "${DEM}_bin_inv = if(isnull(${DEM}_bin), 1, null())"
# shrink DEM by one cell
r.grow in=${DEM}_bin_inv out=${DEM}_bin_inv_grow old=1 new=2
metric=euclidean radius=1.5
r.mapcalc "${DEM}_bin2 = if(isnull(${DEM}_bin_inv_grow), 1, null())"
# area matching the r.contour contours
r.to.vect in=${DEM}_bin2 out=${DEM}_bin_area type=area
# patch countours with outer limit and clean up
v.extract in=${DEM}_bin_area out=${DEM}_bin_bounds type=boundary layer=-1
v.type in=${DEM}_bin_bounds out=${DEM}_bin_lines from_type=boundary to_type=line
# flatten contours for patching, needed to build areas
v.transform in=${DEM}_contour_3d out=${DEM}_contour_2d zscale=0
v.patch in=${DEM}_contour_2d,${DEM}_bin_lines out=${DEM}_patch
# clean up the patched vector
v.clean in=${DEM}_patch out=${DEM}_clean tool=snap,break thresh=1.0e-7
type=line -c
v.type in=${DEM}_clean out=${DEM}_bounds from_type=line to_type=boundary
v.clean in=${DEM}_bounds out=${DEM}_bounds_clean
tool=rmdangle,rmbridge type=boundary -c --v
v.centroids in=${DEM}_bounds_clean out=${DEM}_contour_areas

g.remove rast=${DEM}_bin,${DEM}_bin_inv,${DEM}_bin_inv_grow,${DEM}_bin2
g.remove vect=${DEM}_bin_area,${DEM}_bin_bounds,${DEM}_contour_2d,${DEM}_bin_lines
g.remove vect=${DEM}_patch,${DEM}_clean,${DEM}_bounds,${DEM}_bounds_clean
<---

but I did not find a way to assign a reasonable elevation value to the
areas. The DEM used in this example is in the North Carolina sample
dataset.

The r.mapcalc + r.to.vect approach would be

--->
DEM="elev_state_500m"
step=100

halfstep=`echo $step | awk '{printf "%.15g\n", $1 / 2.0}'`

g.region -p rast=$DEM

r.mapcalc "${DEM}_contours = int(round(($DEM - $halfstep) / $step) * $step)"
r.to.vect in=${DEM}_contours out=${DEM}_contours_rtv type=area column=ele -s
<---

Markus M