[GRASSLIST:5041] Re: Shaded relief + geology ???

[cheg01 at attbi.com]

----- Original Message -----
From: "John Doucette" <doucettj at gbis.com>
To: <GRASSLIST at baylor.edu>
Sent: Friday, November 22, 2002 5:24 PM
Subject: [GRASSLIST:5039] Shaded relief + geology ???


> Hi all:
>
> I am new to GRASS.  I have played a bit with version 4.3 and have version
> 5.0 running on my computer.
>
> What I would like to do is make a map incorporating geologic raster area
and
> shaded relief to come out with a 3-d geologic map.  I have brought vector
> files and dem's into GRASS and have done a bit of fiddling around.  In
> addition I have searched the web and looked at tutorial files and man
pages.
>
> But I really am not sure how to even get started.  Can anyone point me in
> the direction of a method or set of procedures?
>
> John
>

I've done this for every 1:100000 quad in the State of Washington. First
you have to import the data into GRASS in various mapsets, one for each
projection in use. GRASS requires discipline in creating the directory
structure. You also have to be careful with the boundaries and resolution
defined for the region when you do imports or projections of data. Make
sure that the Window (region) files are correct for what you are doing.,
The data I started with was:

Data set                                                        Projection
Format
Washington geology 1:100000  quads          Albers Equal Area
E00
Western Washington 10 meter DEMs            UTM Zone 10               SDTS
Eastern Washington 10 meter DEMS             UTM Zone 11               SDTS
Western Washington 1:100000 vector files     UTM Zone 10               SDTS
Eastern Washington 1:100000 vector files       UTM Zone 11
SDTS

The first hurdle was getting the geology data imported into the AEA
directories. There are some quirks in the E00 format that can cause errors
when the you try to build raster files from them. I wrote some AWK scripts
to fix the errors, you can find them on
http://op.gfz-potsdam.de/GRASS-List/ by searching for PNT_TO_AREA.

The next problem was that GRASS does not automatically maintain the
association between area labels and area properties that you might expect.
Each area is given a unique identification number and the labels are
associated with those numbers in the dig_cats files. If you have 2000 areas
of 200 different types you get 2000 different categories, some of which
share the same label. I used a Microsoft Excel macro to sort the labels and
create reclass files. I used r.reclass and r.calc to create new layers with
all layers with the same label classed in the same category. There is
another story about how to get the colors to be the same between different
quads, also involving an Excel macro. You can probably do the same things
with StarOffice or another Unix spreadsheet.

There are 32 ten meter DEM files in each 1:100000 quad (takes a loooong
time to download them for the whole state). To create a single layer, use
r.in.gdal to import the files and r.patch to combine them for each quad.
Once the DEM files were made, I used r.slope.aspect to create the aspect
file for terrain shading.  I think the default color range for aspect, 0 to
255, is too dark so I adjusted the range to 120 to 255 by editing the
aspect file in the colr directory.

I chose to use UTM projection for the combined quads, so I imported the
geology layers using r.proj. and brought in the vector layers with
v.in.sdts. I was VERY lucky that all of the data used the same datum
references, since GRASS does not do datam shifts. Once everything is in, I
used d.his with h=geology, i=aspect, and s=nothing to create the terrain
shaded geology maps. You can overlay the roads, rivers, whatever you like
from vector files. If you want to save high resolution images for plotting,
I recommend that you use the PNG driver to create the image files since it
will do 24 bit color images at any resolution.

I ended up writing some shell scripts to automate different parts of the
process. They won't work for  anyone else as-is because they depend on
directory names and GRASS environment variables for the particular system I
was on. With the scripts, it took about 18 hours of processing time to
create a quad image on a dual 150 MHz Sparc 20. Probably much faster on a 2
GHz Linux box. Anyway, here is a scripts which may be adaptable to your
project. Be aware that the released version of GRASS 5 treats some of the
environment variables differently from the version I used (pre3), so the
variable references may not work as written. You don't need any scripts to
do any of this, they just reduce the keystroke count..

This one uncompresses and loads the DEM files:

g.region -d
v.proj -s input=bnd location=wa_utm mapset=$MAPSET output=bnd
g.region vect=bnd
d.erase
d.vect bnd
cd $LOCATION
mkdir load
cd 10meter
for file in *.GZ
do
cp $file ../load
cd ../load
gunzip $file
tar -xvf *.TAR
r.in.gdal -e -o input=`ls | grep CATD` output=`echo $file | sed
s/.DEM.SDTS.TAR.GZ//`
rm README
rm *.DDF
rm *.TAR
cd ../10meter
mv $file `echo $file | sed s/GZ/gz/`
d.rast -o `echo $file | sed s/.DEM.SDTS.TAR.GZ//`
done


This one creates the DEM patch file, runs r.patch, runs r.slope.aspect,
brings in the geology data, and plots the d.his image to the screen:

#load_final
    cd $GISDBASE/$LOCATION_NAME
    g.region vect=bnd
    DDIR=`echo $MAPSET | sed -e 's/q//'`
    export DDIR
    ln -s /usr4/wadems/$MAPSET $DDIR
    g.mapsets $MAPSET,$DDIR
    echo "r.patch input=" > patch
    ls -m /gd/wadems/$MAPSET/cats >> patch
    echo " output=dem" >> patch
    cat patch | sed -e 's/,./,/g' > patch2
    sed -n -f /gd/join.sed patch2 > patch
    chmod 777 ./patch
    ./patch
    cp /gd/wageol_aea/geology/colr/dted $LOCATION/colr/dem
    d.erase
    d.rast dem
    r.slope.aspect elevation=dem aspect=aspect min_slp_allowed=0.1
    cd $LOCATION/colr
    cat aspect | sed s/:0/:120/ > asp
    rm aspect
    mv asp aspect
    d.erase
    d.rast aspect
    r.proj input=geol location=wa_100k_lcc mapset=$MAPSET output=geol
    mkdir /gd/wa_utm/$MAPSET/colr
    cp /gd/wa_100k_lcc/$MAPSET/colr/geol $LOCATION/colr
    cp /gd/wa_100k_lcc/$MAPSET/cats/geol $LOCATION/cats
    d.erase
    d.his h_map=geol i_map=aspect

The previous script calls the following script, "join.sed" during the
creation of the patch script to concatenate all the lines into one line for
execution:
#join.sed
H
$ g
$ s/\n//g
$ p