[GRASSLIST:854] contour->dem - nice new tool + questions

[William Kyngesburye]

I found this mentioned in a bug report for r.contour (#4248).  I  
expanded a bit on it and added more details.  Until (if) the nnbathy  
tool gets added to GRASS, the process is a bit lengthy.  But the  
results are better for contour data than r.surf.contour or  
v.surf.rst.  Thanks to Maciek for the original outline of this process.

The main idea behind my expansion of his outline is that I wanted to  
include extra help for the plain contours:

- any spot elevations - either peaks or just extra points in flat areas

- lake polygons of known elevation - the contour can help.  And if  
you like, you can patch in the area afterwards for a flat fill.

It uses the Natural Neighbors algorythm.  Software is here:

http://www.marine.csiro.au/~sakov/

You must get the nnbathy program - binaries for Windows and some  
flavor of linux, or build it from source.

Also, I worked this out using OGR tools.  This should work with GRASS  
6.0.  The current 6.1 CVS has an improved v.patch that could speed  
things up a bit.


Start with 1-3 vectors, lines minimum: all should have the same  
attribute for elevation (elev in this example), and just that  
attribute, though the OGR method will take care of extra attributes  
for you, but the elevation field still must be the same.  I didn't  
take the time to figure out how do this in GRASS.

1. hypsol - At least your contours - all lines, and they don't have  
to be continuous, you can have supplemental contour segments in flat  
areas.

2. hypsox - Spot elevations are helpful - a point vector, can be  
peaks, low points, or anywhere in between

3. hypsop - Lake polygons - area vector.  the boundary will be used  
as another contour, and the area will be used as a mask and a flat-fill


### merge all the layer together as points
# skip steps referring to lake polys and/or spot points if you don't  
have those
# if you have neither, just lines, there is no need to export, merge  
and import
#   and you can skip to the prepaing for nnbathy steps, and use  
hypsox as the point
#   vector as if it were the merged point vector

# spot elevs are ready to go as is
v.out.ogr in=hypsox dsn=. olay=hypsox typ=point

# turn contour lines into lots of points and export them
#   use a dmax = desired dem cell size, but there may be a better choice
v.to.points in=hypsol out=hypsolx dmax=0.000833
v.out.ogr in=hypsolx dsn=. olay=hypsolx typ=point

# rasterize poly fills with elevation to use as mask for  
interpolation grid
v.to.rast in=hypsop out=hypsofills col=elev
# invert nulls/values - this is the grid for the values we want in  
the end,
#   and we don't need to waste time interpolating the lakes
r.mapcalc hypsointerp="if(isnull(hypsofills),1,null())"
r.stats -g -n hypsointerp > hypsointerp1.asc
# remove extra dummy elev column
cat hypsointerp1.asc | sed 's/ 1$//g' > hypsointerp.asc

# also export boundaries of polys as extra contours
#   don't want the whole fill since that would ruin the slope of the  
shore into the lake
v.category in=hypsop out=hypsopc op=add typ=boundary
# same as for lines, pick dmax = cell size
v.to.points in=hypsopc out=hypsopx type=boundary dmax=0.000833
v.out.ogr in=hypsopx dsn=. olay=hypsopx typ=point

# cleanup and merge points, lines and boundaries
#   if there are extra attributes, this takes care of ignoring those
ogr2ogr -select elev . . hypsolx -nln hypsopt
ogr2ogr -update -append -select elev . . hypsopx -nln hypsopt
ogr2ogr -update -append -select elev . . hypsox -nln hypsopt

v.in.ogr dsn=. lay=hypsopt out=hypsopt

## alt merging in GRASS, instead of all the above: CVS only, all  
columns must be exact same
v.patch -e out=hypsopt in=hypsolx,hypsopx,hypsox
##

### prepare for nnbathy

# need coords in attributes
v.db.addcol map=hypso col="x double,y double"
# ! v.to.db really slow !  there must be a better way (or maybe it's  
the dbf)
v.to.db map=hypso op=coor col=x,y
v.db.select -c hypso column=x,y,elev fs=" " > hypsopoints.asc

### finally, interpolate the grid.  at this point, this step goes  
much faster than r.surf.contour of v.surf.rst
nnbathy -i hypsopoints.asc -o hypsointerp.asc > dem.grd

# a few random cells it just can't interpolate, tidy it up for GRASS
cat dem.grd | sed 's/nan/999999.0/g' > dem.grd.nonan

### bring it back in and cleanup

v.in.ascii -t -z in=dem.grd.nonan out=dem_vect fs=space z=3
v.to.rast in=dem_vect out=dem1 use=z
r.null dem1 set=999999.0

# merge in lake polys
r.patch out=dem in=dem1,hypsofills

# you could run r.fillnulls on this to fill in those last stray cells


You can delete all the inbetween stuff now.  The merging step would  
probably be much faster with the v.patch method, but if the vector  
attributes don't match it won't work.  You can't drop or rename  
columns in a dbf, but with sqlite attrib storage it might be  
possible, certainly with postgres.  Just a couple db.execute's to  
alter tables.


If some industrious programmer adds nnbathy to GRASS as a C module  
(or script I guess, it would just be much slower, and require the  
nnbathy binary), it would be nice to have all threee line, point and  
poly source options, or at least line and point.  To make it easier  
to deal with attribute disparities, all sources should have a way to  
set the elevation column (or to use the z value).


QUESTION:

The output looks great - a coworker said it looks a lot like Arc's  
topogridtool output without any stream input.  So, is there some way  
to enforce stream drainage patterns (if one has the digitized stream  
vectors) in generated DEMs with GRASS commands?  Could it be done  
after the above DEM interpolation, or would it only work as a part of  
the interpolation (assuming someone could add that to the program code)?

-----
William Kyngesburye <kyngchaos at kyngchaos.com>
http://www.kyngchaos.com/

All generalizations are dangerous, even this one.