v.import creates a GRASS vector map from the ASCII DLG data with the same name. (Note that v.import does not import the spot elevation points in the DLG file. When creating a DEM, unless these denote peaks or bottoms of depressions, they are best omitted anyway, since they have unduly large influence on the interpolation algorithm.) v.import throws a benign error and core file as it terminates; just rm core afterward.
Display each vector map to make sure it imported okay, e.g.:
g.region vect=newark_e.hyp
d.erase
d.vect newark_e.hyp
Run v.support on each of these vector maps to create the topology for them (this creates the file in the dig_plus directory that matches label points in the dig_att directory with the contour lines in the dig directory). At this point the contour lines are labeled with the ID numbers 10001, 10002, etc. You will reclass these ID's to elevation values later on.
Now create raster versions of the two vector maps. The current
region settings and MASK (if any) don't matter for vector maps, but they
do control the resolution and extent of raster maps you create. So
set a 30-meter resolution region (the standard for 1:24,000-scale DEM's)
for each quad before rasterizing your hypsography files, e.g.:
g.region vect=newark_e.hyp res=30
v.to.rast i=newark_e.hyp o=newark_e.hyp
If you d.rast the raster maps produced by v.to.rast, you will see random line colors reflecting the randomly-assigned ID numbers.
10026 20 200 22 70These are 10-foot-interval elevation contour lines. On near-vertical slopes several contours may converge to a single line which will be assigned multiple elevation values. For example, line 10027 might be a convergence of 90-, 100- and 110-foot contours:
10027 20 201 22 90 22 100 22 110If you are good with AWK or have a text editor like textedit which can do global find-and-replace of text strings that include line-feeds, edit the two attribute files into reclass files to be redirected into r.reclass, e.g.:
10026 = 70 10027 = 90 100 110Note that the second line here is assigned a reclass value of 90 feet; the extra elevation values are retained as a label.
If you have much trouble reformatting the attribute files into reclass rules files, you can simply use the files newark_e.att and newark_w.att which have already been reformatted. These are also located in the "dlg" subdirectory of PERMANENT's mapset. These will work fine as long as your west quad contours are labeled starting with 10001 and your east quad contours are labeled starting with 20001.
Now r.reclass the raster contour maps to elevation values,
redirecting in the reclass rules files reformatted from the attribute files
created by v.import, e.g.:
r.reclass i=newark_e.hyp o=newark_e.hyp.recl
< newark_e.att
The (reclassed) raster hypsography maps will display the familiar GRASS default color ramp.
Write down the bounding coordinates of both raster hypsography maps' default regions. Then manually re-set your region (g.region option 1) so it just encompasses both maps, with a resolution of exactly 30 x 30 meters. Since g.region modifies the cell resolution slightly to fit any set of bounding coordinates, you may have to tweak the bounding coordinates as necessary to get an exact 30 x 30 cell resolution. Now r.patch the two reclassed raster hypsography maps together. Display the patched map to check it.
(Note that the preceding three steps--rasterize, reclass and patch-- can be done in pretty much any order, since GRASS also has v.reclass and v.patch utilities. For example, you could v.patch the vector hypsography maps together, then cat your two reclass rules files into a single file and v.reclass the combined vector map, then v.to.rast that.)
Creating the DEM
Once you have a patched map of both quads with a consistent color scheme, create the DEM with r.surf.contour. This module takes a minute or two.
The DEM is created without a specified color table in the colr directory of your mapset, so it gets the default ("rainbow") color table as shown here. Try designing a better color table with r.colors by specifying rules. You can get a "purple mountains majesty" effect by specifying the sequence of standard GRASS colors: aqua, green, yellow, orange, brown, purple, magenta, white.
1 thru 3 = 1 0-2% slope
4 thru 6 = 2 3-5% slope
7 thru 11 = 3 6-10% slope
12 thru 16 = 4 11-15% slope
17 thru 26 = 5 16-25% slope
...
Use r.neighbors "mode" method to smooth this slope map using a 3x3 window..
Create a hill-shaded DEM with d.his, using the DEM for hue and the aspect map for intensity (nothing for saturation). You can get a brighter version (on which vector overlays will stand out more clearly) if you edit the color table of the aspect map (located in the colr subdirectory of your mapset) with a text editor. The default aspect color table has two ramps: 0 degrees aspect = 0 (black) to 180 degrees aspect = 255 (white), and 180 degrees aspect = 255 to 360 degrees aspect = 0:
% 1 360
0:0 180:255
180:255 360:0
Use a text editor to change the minimum values to some lighter gray, say 100 or 128, and then re-run d.his so the hillshading isn't so dark, outputting the image to a new map.
Drape the hill-shaded DEM image over the source DEM with d.3d. Try to get a nice viewing angle with enough vertical exaggeration so the separage drainages are visually intelligible but the hills aren't ridiculously spiky. (Keep in mind that while the X and Y units are meters, Z (elevation) is in feet, so a vertical exaggeration of 0.3048 would display the real proportions of the terrain.
r.watershed takes some time, so after you get it launched, you can background it. Maybe check out the module's options in the g.manual entry on r.watershed).
The basin map output by r.watershed identifies each basin by an even integer ID value. r.watershed gives the basin map a "random" color table. If you use r.colors to assign a "rainbow" color table instead, you will see how basins are numbered numbered sequentially from the outlet of each major drainage. Create a hill-shaded version of your basin map by compositing it with the aspect map using d.his. Try a 3D drape of this.
The accumulation map calculates the number of cells draining through
each cell. Negative values indicate the cell receives some off-map flow
also. When you d.rast this map, it shows the White Clay Creek, Christina
River and Elk Creek channels, but not much else. Use r.mapcalc to
calculate the log of the absolute value of the accumulation map, e.g.:
r.mapcalc log_accum='log(abs(accum)+1)'
and d.rast this to get a clearer view of these drainage systems.
Use r.mapcalc to extract inferred streams from this log accumulation map: a good threshold value is 6: r.mapcalc inf_streams='if(log_accum>6)'
Create a vector layer of inferred streams. First, use r.thin to thin the binary raster streams map to 1-pixel thickness. Then use r.line to vectorize the streams. Then run v.support to build the topology for this stream map.
Overlay your vector map of inferred streams on the basins map and identify the basins associated with the Elk, White Clay and Christina watersheds. Reclass the basin map to distinguish just these three main watersheds. Calculate the total map area of each of these three watersheds. Use r.poly to create a vector map of these watershed boundaries. Run v.support to build topology for this file. Create some nice final display maps.
Extra Credit: Zoom to a hilly portion of the DEM, perhaps 4KM by 4KM, and save this as a separate DEM. Using the watershed script in the Shapiro & Westervelt tutorial r.mapcalc: an algebra for GIS and image processing as an example, write a shell script to simulate the iterative down-gradient flow of a chemical spill from a randomly-located point on this DEM. Have the script capture a user's mouse-click from d.where locating the spill point. The script should output a map for each iteration or two. If you convert each iteration map "flow.n" to a .GIF format file with the command r.out.ppm i=flow.n o=- | ppmtogif > flow_n.gif you can animate these with a GIF animator.