representing geology with rasters?

[Malcolm Williamson]

Tom Charnock writes:
> 
> I am trying to link GRASS to a ground-water flow model.  At the
> moment am a considering the problem of representing geology within
> the GRASS database. The model I am using divides the geology up into
> 3d cells with constant properties (eg hydraulic conductivity) within a
> cell.  Rasters are two dimensional, so the only two ways I can see of
> representing geology are as a series of planar slices through the
> geology,  or as a series of "DEMS" following strategraphic surfaces.
> Either way I can see a lot of problems looming, particularly with
> maintaining consistency between raster layers as GRASS has no way
> (that I know) of linking rasters except by giving them sensible
> names; eg layer_1 layer_2 etc
> 
> At the moment I favour using strategraphic surfaces as this fits the
> model I am using most closely, and GRASS is quite rich with surface
> interpolation routines.  I am concern however that when I
> change region, resampling will lead to layers intersecting, a kind of
> 3d analogy to the problem of slivers in vector data.
> my head hurts :-(
> 
> Has anyone addressed this problem before?
> 
> You can probably tell I'm no geologist or database engineer, comments,
> discussion and references welcome :-)
> 
> much cheer    Tom.

I haven't worked with models quite like yours, Tom, but I have given this 
quite a bit of thought. Think of your collection of raster layers as being
a conventional tabular database. Each spatial location, represented by a cell,
would represent a database record. Each raster layer holds a piece of 
information about that record - a field, if you will. The name of that raster
layer would be the heading for the field. You add additional fields by 
creating new raster layers.

You might use one raster layer to define the upper limit of a geologic layer,
and a second raster layer to define the bottom, both in absolute elevations.
Another raster layer would hold one of the constants, such as hydraulic
conductivity. As additional pieces of information are needed for each spatial
location, new raster layers are created. r.mapcalc will be your primary 
analysis tool, enabling you to conditionally mix and match the different
raster layers. As long as you are ultimately interested in outputting 
*numbers* from your model, this approach should work well. If you want
to display data in three dimensions, you will have a tough time with it (unless
you have a Silicon Graphics machine, enabling you to use SG3d). You'll just
have to be creative.

I believe that your concerns about "slivers" are unfounded, since you are
really only dealing with values in a fixed matrix.

Hope that this gives you some ideas!
-- 
Malcolm D. Williamson - Research Assistant       E-mail: malcolm at cast.uark.edu
Center for Advanced Spatial Technologies      Telephone: (501) 575-6159
Ozark Rm. 12                                        Fax: (501) 575-3846 
University of Arkansas              
Fayetteville, AR 72701