On Thu, Jan 9, 2014 at 9:57 AM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:
On Thu, Jan 9, 2014 at 1:46 AM, John Ciolek <jciolek@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@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