[GRASSLIST:3589] Re: Geophysical framework: WHITE PAPER

[Michael Barton]

Craig and Benjamin

Your thoughts are much appreciated. This project is really beginning to 
take shape. I have a couple more suggestions from the GRASS 
perspective.

It is easiest to include them below.

>
> Thanks Michael and Benjamin, I have several comments/questions:
>
> Acquisition/Raw data:
> <quote>
> 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.
> </quote>
>
> In my experience, the data is rarely uniform. For mineral exploration, 
> data is often collected in less then hospitable terrain which means 
> that certain points on the proposed acquisition grid cannot be 
> accessed in the real world. For the surveys I have been involved in 
> where the instrument does not have a GPS, a coordinate was obtained 
> using a technique appropriate for the accuracy desired ie GPS 
> instrument or using laser Theodolites. Even in this case though, the 
> data could still be imported as vectors as you described after some 
> manipulation:
>
> x-coordinate, y-coordinate, other data (ID, values, dates, etc)

Thinking it over, there are at least 3 possible kinds of data sets:

1. Data with georeferenced xy coordinates, such as those derived from 
intruments with integrated GPS (or post processed to include GPS data).

2. Data with xy coordinates that are not georeferenced, such as those 
from an arbitrary grid laid out in the field

3. Data lacking xy coordinates that need to have them assigned.

If the data flow involves input points into vector format, optionally 
assign coordinates if needed, optionally georeference is needed, it 
would accomodate all three.



>
> Regarding corrections that need to be performed on the data (also 
> often referred to as reductions):
>
> 1) Drift correction.
> 2) Latitude correction
> 3) Eötvos correction.
> 4) Reduction to the pole.
> 5) Free air and Bouger correction.
> 5) Terrain corrections.

These all seem conceptually be be errors in measurement that need 
correction. If so, they should be corrected in the original 
measurements (i.e., vector points) prior to surface model generation. 
Georeferencing does not affect these, but would be necessary at least 
to apply terrain correction (i.e., because DEM's are usually 
georeferenced).


> 6) Trend surface removal. This applies to grav and mag data. Usually 
> involves fitting a 1st or 2nd order surface to the data and then 
> differencing that surface with the data set to get the residual data. 
> There are many ways to do this. I recently came across a paper from 
> Paolo Zatelli and Andrea Antonello titled "New GRASS modules for 
> multiresolution analysis with wavelets". Multiresolution analysis 
> would be a great way to do this reduction, although a more traditional 
> approach would be to use least squares to fit a 1st or 2nd order 
> surface.

I **think** that this is better done on the surface model than the 
original points because it is a type of signal processing (see below).

>
> To summarize: All the standard corrections (except maybe terrain 
> correction) could be applied using existing GRASS modules. These 
> corrections could be implemented either in shell scripts or in a new 
> GRASS module. My inclination is to write a GRASS module to do this. 
> Any comments?

There is no map calculator for vector points. However, they can be 
altered using SQL update queries within GRASS. The calculations and 
data manipulations necessary to do the corrections might be better 
implemented in a GRASS module than a shell script.

> NOTE: The last message I sent got garbled up,
> so please disregard and read this one instead:

Thanks for resending this. I thought it was my email.

>
> (a) can very easily be converted to (b) by adding coordinates, so it 
> would
> indeed be easy to have a unified data model.

I agree. That's why I came to the same conclusion about using vector 
points after some thought.

>
> For data type (a), rasterisation can be done without interpolation and 
> the
> filters should be run on this, un-interpolated, data.
> Interpolation will only be necessary in the final step, if the axes 
> and/or
> resolution of the working location are different than those of the 
> data.
> So -- just one interpolation in this case.
>
> For type (b), interpolation will have to be done right after import, 
> because
> ** I think ** filters will work better on data which is still properly 
> aligned
> to the region's X/Y axes. I have a hard time, e.g. imagining a 
> destripe filter
> work well on stripes running diagonally, which might happen,
> if the filter is to run on a rectified image (correct me, if this is 
> wrong!).
> In this case, we would indeed need two interpolations for the final 
> result.

This is somewhat complicated. With respect to interpolation, 
conceptually all of the rasterization should be by interpolation. All 
geophysical measurements are point samples of a continuous field. 
Hence, you need to model the continuous field from the samples. While 
there are a number of ways to do this, interpolation is the most 
common.

If a set of points are orthagonally distributed and equidistant from 
each other, you can certainly directly convert them to raster grids 
(***v.to.rast***). However, this assumes that the field in the space 
between points is constant up to the edge of a grid square (half-way 
between points), where it abruptly changes to the value of the next 
point. In some cases, where orthagonal, equidistant sample points are 
so closely spaced that inter-point variation is irrelevant, 
interpolation is unnecessary. However, I suspect that this is the 
exception rather than the rule. Direct rasterization should probably be 
an option, but interpolation should probably be the norm.

The second issue is what order to interpolate, filter, and 
georeference. In my experience with satellite imagery, georeferencing a 
raster image causes it to lose visual clarity a bit--that is, signals 
lose definition. Hence, signals enhanced by filtering will lose some of 
their visual impact (and possibly quantitative definition and 
separation) if a raster is georeferenced after enhancement. This may or 
may not be a problem, depending on what one is doing with the raster. 
This does not happen, however, if the raster surface model is created 
after the vector points are georeferenced. That is, georeferencing a 
raster will change the values of a raster grid, but georeferencing a 
set of vector points will only change the xy coordinates, not the z 
value.

While most filters should work equally well on georeferenced and 
non-georeferenced data sets, destriping may be another problem. I 
thought that a fft filter destriped by focusing on signal frequency 
rather than orientation, but I have very little experience with this 
and may well be wrong here. I suggest some experiments with filtering 
and georeferencing to see what order produces the best results.

The other reason to georeference early is that georeferened point sets 
(e.g., from multiple sample grids) can be easily combined visually 
simply by displaying and combined into a single analytical file using 
v.patch. The same is true of georeferenced rasters. They can be 
mosaiced by using r.patch. There is no need for the 'r.composite' 
module mentioned.

One potentially useful enhancement of r.mfilter is the ability to run 
the filter across multiple rasters in a single pass. Currently, a 
filter only alters a single raster. This could be enhanced so that 
r.mfilter could be run across a set of rasters specified in the input. 
Also note that r.mfilter allows the specification of multiple filters 
to be run across a raster.

>
> I would like to keep the grid layout information intact, because I deem
> it very handy.

This is a good idea. I'd suggest keeping it in a simple text file with 
structured variables much like a GMT or ESRI header file. Thus, each 
data set (i.e., from a grid sample) would be stored in a text file 
named <filename>.grd (or something like that) and a corresponding grid 
definition file would be stored in <filename>.gdf. A *.gdf file might 
have a series of variables like

TITLE
DATE
ROWS
COLUMNS
ULX
ULY
ULXCOORD
ULYCOORD
LLX
LLY
LLXCOORD
LLYCOORD
...
and so on.

Such a file could be easily parsed by a module that would read it and 
determine if it was type 1,2, or 3 above and process it accordingly. It 
would be kept associated with the data file by having the same name but 
a different extension.

I hope this is helpful.

MIchael Barton

____________________
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 Jun 1, 2004, at 10:01 PM, Multiple recipients of list wrote:
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: text/enriched
Size: 9269 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-user/attachments/20040603/cd3abe6e/attachment.bin