[GRASSLIST:3558] Geophysical framework: WHITE PAPER

[Michael Barton]

Benjamin,

Your white paper clearly represents considerable thought into the needs 
for processing geophysical survey data. This is a good beginning for 
looking at GRASS from a geophysical survey perspective. I think it is 
great that you are serious about pursuing this and enhancing GRASS for 
this aspect of archaeological and geological research. The value of 
using GRASS as a platform for processing geophysical data is that it 
provides a rich and complex tool set for the analysis and management of 
spatial data. However, this richness comes at the price of a rather 
steep learning curve (though not any steeper than other full-featured 
GIS systems). Because you are at the early stages of learning GRASS, I 
hope to help by offering ideas about how to best use the existing tools 
  in GRASS to further this objective.

First, I want to be sure that I understand the overarching objectives 
of such work. Please, Benjamin and anyone, feel free to offer 
amendments or corrections here. As I see it...

1. Geophysical equipment in the field record data at a series of points 
within a survey area. These points theoretically sample a continuous 
field of some parameter (earth's magnetic field, soil conductivity, 
radar reflections, etc.). Usually values are recorded at some set 
'depth' below the surface (or depth range) for all points in a given 
survey. But in some cases, GPR or electromagnetic tomography for 
example, values may be recorded for multiple depths at each sample 
point.

2. A model of the continuous field must be constructed from the sample 
points. Usually, this model takes on the format of a continuous surface 
and is digitally represented by a regular grid that is conceptually 
equivalent to a digital elevation model (i.e., a raster grid of 
equal-sized square cells, in which each cell has a value of the field 
at that point).

3. This field model--or surface--needs to be transformed to enhance 
specific signals of potential interest. This transformation commonly 
involves moving window convolving filters, fast fourier (for 
destriping) transformations, and whole image transformations (e.g., 
contrast enhancement or gamma correction, Laplacian edge detection, 
etc.).

4. The transformed field model needs to be georeferenced so that it can 
be integrated with field models from other geophysical surveys and 
different types of archaeological or geological data.

Given these objectives, I offer some suggestions as to how to integrate 
existing GRASS modules  with some of the new ones you are proposing. 
I'm attaching a text file with these suggestions so as not to make the 
text in this email overly long.

Cheers
Michael Barton

-------------- next part --------------
GEOPHYSICAL PROCESSING AND GRASS
Michael Barton, ASU
(I've emphasized GRASS commands as ***command*** to make them easy to see in plain text)

PART I - Raw Data to Surface Model

1.A Raw Data to Georeferenced GRASS Vector Point File

The first issue is how to get raw geophysical data into a GRASS GIS format for subsequent processing. Hamish recently suggested a way to do this with raster data. However, after thinking about it, I have come to favor initially importing raw data as vector points (using ***v.in.ascii***) for several reasons I'll outline below. To import raw data in this way, it needs to come in a format like this...

x-coordinate, y-coordinate, other data (ID, values, dates, etc)

Some equipment, like the Cesium magnetometers that have been used on my sites recently, automatically generate xy coordinates for each data point; others do not. For data sets that do not come with xy coordinates, it will be necessary to assign xy coordinates based on transect spacing (x-coordinate) and transect length + number of data points/transect (y-coordinate). It would be convenient if a data input module could assign x&y coordinates to the original ascii raw dataset if needed

Unless the instrument incorporates a GPS, the xy coordinates will be geographically arbitrary for a given dataset. Hence, sooner or later it will be necessary to georeference these data. Although this can be done on the initial vector data points or at the end on the analyzed raster surface model, I suggest that it is analytically better to do this on the points prior to transforming them to a raster surface model. Doing it later requires 2 interpolations of the data values (once to create the surface model and a second time to georeference), potentially degrading data quality.

If the data are a set of vector points recorded in an orthogonal grid, they can be georeferenced using ***v.transform*** and a simple first order transformation. This requires only that the real-world coordinates (e.g. UTM's) be known for 3-4 of the original data points (i.e., ground control points or GCP's). The corners of the original survey area are convenient places to do this, but any well-spaced set of points will do in theory. 

The workflow of this data input module would be as follows:

Import ascii raw data into GRASS vector points using ***v.in.ascii*** [optionally, assign arbitrary xy coordinates prior to import]. These data are imported intially into a temporary xy location. Then ***v.transform*** is run on the data, employing user-supplied GCP's to georeference them into the current working, georeferenced GRASS location. The temporary xy data can then be deleted to keep from cluttering up the work. There is some overlap with the functions of r.in.grid module proposed in the white paper. Perhaps it could be modified to make use of exiting GRASS functions in this way.

Import module arguments would include:

input file name, output file name, GCP's (UL x,y,east,north; LL x,y,east,north; UR x,y,east,north; LR x,y,east,north), [transect spacing, transect length, data points per transect], field to use as x-coordinate, field to use as y-coordinate.

The result would be a georeferenced file of vector points, representing the real-world location of each sample spot (each data point from the field), and containing the value of the geophysical parameter sampled at each spot. 

Because it is georeferenced, this geophysical dataset can be easily joined with other geophysical datasets at this time, using ***v.patch***. Any raster surface models created from the dataset can be mosaic'ed using ***r.patch*** (see also **r.mpatch** and i.image.mosaic scripts). Either can easily overlay or be overlayed on any other georeferenced archaeological or geological data.

This does not deal with 3D data. If we don't want to explicitly include G3D data at this time, it is still worth considering how to develop such functionality so as to permit the addition of additional dimensions at a future date.

1.B - Creating a Surface Model

Because conceputally geophysical data represent point samples of continuous field phenomena, interpolation seems the best means of creating a surface model from the original point data set. 

GRASS currently includes interpolation modules for bilinear (***r.bilinear***), inverse weighted distance (***v.surf.idw***, with a choice of several algorithms), and regularized spline tension (***v.surf.rst***). It lacks modules for trend surface creation and krieging (regularized spline tension is computationally is equivalent to krieging for surface creation, but does not include the theoretical/analytical component of computing a semivariagram and fitting it to a model). 

Surface model creation could be part of the input module, be a separate module incorporating various interpolation choices, or simply use existing (or new, e.g., krieging) GRASS interpolation modules. It depends on whether or not there are effective ways to optimize interpolation for surface model creation in geophysical survey contexts. If part of the input module or a new comprehensive interpolation module, arguments would include: input file, output file, interpolation method, [optional interpolation parameters]. 

Note that surface module resolution can be specified on the fly via region settings (***g.region***) or changed permanently in the raster file via ***r.resample*** or the map calculator.

Currrently, there is only regularized spline tension to create a volumetric 3D model from vector points (***v.vol.rst***).

ADDITIONAL COMMENTS ON WHITE PAPER

In the white paper, I was confused by the use of the term grid. It seemed to refer alternatively to a GRASS raster data file, an ascii point file or its derivative, and a parcel of landscape surveyed. I suggest referring to a dataset from a single geophysical data collection activity as a quadrant, block, dataset, or other such term; I suggest referring to GRASS raster files as simply raster maps or raster files; "grid" can be then used to refer to the original organization of the transects and data collection sample points in the field. 

Original raw data can be stored in any convenient locality. I'd argue against storing raw data in the $GISDBASE folder unless there is compelling reason to do so. To keep from mixing things up, this should be reserved for GRASS data only. Also, the structure of raster file storage may change in the near future as GRASS 5.7 continues to be developed. Importing the raw ascii data from geophysical instruments as vector points obviates the need to define a separate GRASS file structure, $GISDBASE directories, or a new GRASS data type. Given the complexity of GRASS currently and the fact that it includes vector, raster, and 3D data models (covering all those currently used in GIS work), I don't think it is a good idea to create new data models unless there is a compelling reason to do so. Vector points match the data model of the original field sampling points and GRASS raster maps serve well as the surface model created from the data points. 

It seems that the functions of proposed m.grid.layout.create and m.grid.layout.edit are either already handled by existing GRASS modules or not necessary.

Combining datasets (or survey block) can be done using ***v.patch*** for the vector  points or ***r.patch*** or ***r.mpatch (a script)*** for raster files. Note that r.composite (proposed in the white paper as a patching module) is already a GRASS command module to produce composite rgb color images/maps from individual images/maps for the red, green, and blue channels. Again, the reason to do this kind of work in GRASS is to take advantage of its already developed data models and management/analysis routines, avoiding starting from scratch as much as possible.

PART II - Transforming Surface Models 

Many tools for transforming surface models exist in GRASS. However, it would be very helpful to identify and perhaps enhance those most useful for processing geophysical data. 

Moving window filters in grass include ***r.neighbors*** (with choice of window size and automatic calculation of average,median,mode,minimum,maximum,stddev,sum,variance,diversity, and interspersion) and ***r.mfilter*** (a basic convolving filter module in which the user can specify iterations, along with filter size and values via a text file). These modules apply to any raster file; their effects can be spatially limited via MASK and region settings. As suggested in the white paper, set of predefined filter files for processing geophysical data would be helpful. Even better would be a new GUI interface in which users could interactively vary the parameters for various filters (e.g., by sliders), also suggested in the white paper. Better yet would be a preview window in which users could see the effects of the filter before applying it. I don't know how possible the latter is given the interaction of GRASS commands with the x-display drivers. Devising a series of sequential filters, as suggested in the white paper, is a useful idea if there is a 'normal' set that can be run without examining the results after each filter. ***r.mfilter*** currently can be set to rerun the same filter a number of times. It should be easy to modify r.mfilter to run a sequence of filters. The r.filter module described in the white paper seems very similar to the existing GRASS ***r.mfilter*** module, especially if enhanced in this way. However, there are also other types of transformations to consider.

Transformations that affect the entire raster can be done by the map calculator (***r.mapcalc***); "canned" transformations that are directed at geophysical data processing could be created. Fast fourier transformations are done by ***i.fft*** and ***i.ifft*** (to reverse the transformation). Destriping or other kinds of processing can be done by masking during the inverse transformation. A new module could include i.fft, i.ifft, and appropriate masks (or better, an interactive, GUI based mask creation function) for tasks such as destriping. GRASS includes a fairly sophisticated edge detection filter--***i.zc***-- that combines FFT, Laplacian, and Gaussian manipulation of images, with considerable user control. Histogram stretching (equalized and log) is done via color table support (***r.colors***). However, there is no interactive (or non-interactive) module for other kinds of histogram manipulation like "contrast" or "brightness" controls, or gamma correction (although they could be done via the map calculator). 

All (or nearly all) transformations available in GRASS can be spatially controled in two ways--via region settings that set the extents of the transformation and via a MASK. The latter is a raster file of any shape. Transformations only occur in parts of a raster file that coincide spatially with the cells of the MASK file.

A masking function is available for G3D data (***r3.mask***); the functions of a separate 3D map calculator have just been incorporated into ***r.mapcalc***. No other analytical functions are yet available for 3D datasets. This is an important area for development highly relevant to archaeology and geophysical survey.


ADDITIONAL COMMENTS ON WHITE PAPER

I think the most important first step is to identify exactly what kinds of transformations/filters are most important to geophysical data processing. Then we could better assess whether these could be best accomplished using existing GRASS modules, enhancement of existing modules, or the creation of new modules. You mention destripe, upward, downward, highpass, lowpass. All of these can be accomplished by r.mfilter and i.fft/i.ifft. Are there others?

PART III - Integration of Processed Geophysical Data with Other Data (including other geophysical data)

If geophysical data are georeferenced at the outset, this task is inherent within the GRASS GIS framework. Geophysical data sets only need to A) be collected within an orthagonal grid for which xy coordinates for each data collection point can be assigned, and B) have real world coordinates known for 3-4 locations on the grid. Datasets do not have to be rectangular, just collected within a rectangular grid. Any other data that can be georeferenced to ANY standard coordinate system can be combined with the analyzed geophysical data using the GRASS reprojection modules ***r.prov*** and ***v.prov***. This leaves a great deal of flexibility for data collection in the field and potential for use subsequently.

Datasets can be patched together at two different stages (the georeference vector points or the raster surface model). GRASS script i.image.mosaic provides an idea of how this can be accomplished for raster files. A new function that would be handy would be the ability to develop a color table using the combined values of multiple datasets that have NOT been patched together. This would allow them to remain separate but be displayed as if they were a single, combind data set.

The abilities of r.colors should be reassessed; enhancements could be proposed if needed. A GUI version is highly desirable here, and possible combination with information derived from r.univar or r.statistics. 

Finally, I can't think of any needed enhancements to current GRASS display capabilities (d.rast and NVIZ) that are specific to geophysical survey data processing, but perhaps others can. 

-------------- next part --------------

____________________
C. Michael Barton, Professor
School of Human Origins, Cultures, & Societies
PO Box 872402
Arizona State University
Tempe, AZ  85287-2402
USA

Phone: 480-965-6262
Fax: 480-965-7671
www: <www.public.asu.edu/~cmbarton>
On May 29, 2004, at 10:01 PM, Multiple recipients of list wrote:

> as promised, this is the white paper in which I attempt to
> outline a unified framework for implementing geophysical
> functionality in GRASS.
>
> What is missing:
>
> - anything relating to 3D analysis
> - color ramp modification tools
> - a full list of the actual filter functionality we want to
>   implement (see section at bottom of document),
>   should also include a short description of how each filter
>   operates and what it is useful for
>
>
> I am attaching a PDF version for viewing and a M$Word version
> for your edits.
> I would like everyone interested to add comments, make changes
> and additions, look for inconsistencies (terminology), spelling,
> grammar (I am not a native English speaker), clarifications
> - please mark your changes with another colour in the DOC file
> and send it back to me, so I can synchronise all the edits.
> Craig would then be able to forward the white paper for review.