Tools and programs- please acknowledge as appropriate:

Plot shaded gOcad t-surfs using GMT for easy manipulation in graphics packages.

gOcad is used by our research group for structural modeling of triangulated surfaces. Often, we make images of these surfaces by either the "old screenshot" method or a direct export of the camera view to a raster image (TIFF or JPEG). I often am not happy with the results, as the screenshot is lossy and unscalable. And to be honest, since I mostly work on my Macintosh, I hate to have to go back to the lab, open up the project again, get the right view, with all the right elements to shoot another screenshot if I find errors in the image of the surface or can't scale it.

So I developed a quick routine to plot shaded t-surfs as Postscript using GMT (Generic Mapping Tools). I use a flat shading technique that gives the illumination based upon the angle between the illumination direction and the normal of the surface (use your favorite search engine to find further info on "flat shading"). This technique was developed to make maps of gently undulating surfaces (geologic horizons), but I show an example of it for a steeply dipping surface. There are many caveats like which way the normals are facing in your calculation, etc that may not make the results consistant from surface to surface, but this is a start.

This script fault_plot.awk (control - click) generates the multi-line file for use in GMT. It reads through the .ts file and takes the triangles and calculates the illumination (see above). It then spits out a GMT ready file with the color values (angle/illumination) in each header of each triangle.

Once you do this, you have a file that can then be plotted. HERE is an example file from the Hollywood fault as per the CFM. In this example the light comes from 0,-1,0 (north direction) and it is plotted in perspective (the -E in the GMT command). You will also need a color file HERE that scales the illumination to some user defined RGB values. This is the trickiest part and often first requires inspection of the angle values to find your min and max illumination.
This can be done with your t-surf by using the awk script and the GMT command minmax:

./fault_plot.awk yourtsurf.ts | awk '{if ($1==">") print $3}' | minmax


The command in GMT I then use to plot the t-surf in 3D space:

psxyz 3d_plot_examp.xyz -JX10/4.1698495 -JZ8.6273 -R367906/390161/3772140/3781420/-19000/200 -E330/15 -M'>' -L -P -W1 -Cplot.cpt -B10000/10000/5000 > plot.ps



You also need to know the X-Y-Z scale/limits of your t-surf for the -R command. Notice my -J commands scale 1 to 1 to the extent given in -R. This postscript and the surface elements are now directly editable in a graphics package of your chosing. That gives you this:
gmt

Surface plot from GMT script with shading on the triangular elements from a northerly light source. Note that some of the triangular elements seem to have 0 area, which may lead to problems in the calculation step.		 gocad

Same surface in gOcad. Note this was only shot to illustrate the suface extent, not the resolution which is variable and can be changed.

5/2006

How to restore seismic data using a structural model following inclined shear restoration techniques in Adobe Illustrator is HERE.

1/5/2006

How to convert an ERDAS Imagine image to a geospatially referenced NETCDF .grd file for use in GMT (Generic Mapping Tools)

We commonly use images (LANDSAT, geologic maps, etc) in our research and often have to change the geographic projection of the image or actually take a TIFF or JPEG and geospatially reference it. I use the Generic Mapping Tools to make many of my figures and wanted a way to incorporate geospatially referenced images with other datasets (earthquakes, shotpoint maps, etc).

Here is how you do it:

  1. First convert the .img file to Sun raster files and export it as something like blah.rs in ERDAS
  2. Use gmt2rgb to make individual rgb files

    gmt2rgb blah.rs -Gcolor%c.grd -R0/3144/0/2580 -F -I1/1


    Where 3144 are the columns and 2580 are the rows
    This makes 3 R,G,B color files that can be used in grdimage
  3. Plot a test image

    grdimage colorr.grd colorg.grd colorb.grd -R0/3144/0/2580 -B100 -JX10 > plot.ps

  4. Now reconvert the .grd files to a geographic coordinate system. To do this you need to know the corners of the image (from ERDAS- image viewer of the original image, get layer info).

    grdedit -R-121.664/-118.52/35.061/37.64/35.06 colorr.grd


    where the X distance at 3144 is equal to 0.001 degrees spacing and the Y distance at 2580 is very close to 0.001 spacing.
  5. Do this for all 3 color .grd files
  6. Plot image with other geographic datasets.
8/9/2005

gOcad script (in awk) to generate a group of surfaces from user-defined locations with a user defined strike and dip
How to use:
This program, make_plane.awk, is written in awk so that it can be read on multiple platforms. I have no plans to port it to other platforms.
The program reads a tabulated ascii file that has the following information about a point:
X, Y, Z, strike, dip
where X,Y,Z are in user defined units and strike and dip are in degrees following the Right Hand Rule (RHRule).
It then makes a surface (in a gOcad group) centered at that point with the given strike and dip value, with the dip of the surface being a surface property.
The original program was generated for gOcad v. 1.5 but should work for later versions as well.

get make_plane.awk (control - click)
get test data (control - click)
get the group of surfaces made by the test data (control - click)

coalinga_planes

Subset of fracture data measured from Coalinga anticline (2002 CAG)
Contact me: guzofski@fas.harvard.edu with any questions or problems (which I am sure there will be!). However, use at your own risk, I claim no responsibility for problems or errors caused by the software or materials they generate. And ALWAYS SAVE YOUR PROJECT!!!!!
Useage:
  1. In a text editor go to line 22 of the file make_plane.awk and change the dimensions you want the surface to be. The variable DIST is the half length of each side of the plane to be generated. Currently it is set to 5 user defined units.
  2. In a text editor go to lines 44 and 45 of the file make_plane.awk and change the amount subtracted from your X and Y values. This is done to make the computations more robust. In the default settings the program subtracts 730000 from the X values (so to make, say, 734956 a value of 4956 for the computations) and subtracts 4010000 from the Y values (so to make, say, 4010020 a value of 20 for the computations). If you are dealing with small numbers in the first place you can remove these subtractions from the awk file.
  3. You may (but don't have to) change the name of the group (line 27- default is calc_plane) or the name of each plane (line 180 - default is plane_#). You may also need to change the first line from #!/usr/bin/awk -f to the path where you have your executable awk file.
  4. In some shell where you can execute awk files type:
make_plane.awk datapoints.txt > datapoints.gp
where the file datapoints.txt is the XYZ data file and the file datapoints.gp is the gOcad Group file
From here you just have to load the Group file into your gOcad project and then go to the pulldown window in gOcad for : Compute -> Apply Script -> On Object and do the following:
{X = X + 730000}
{Y = Y + 4010000;}
or whatever the value you subtracted in step 2.
to return the X and Y locations back to their original position.

12/20/2002

Program (in gawk) and spreadsheet (in Excel) to calculate the strike and dip of a surface defined by three points (following Bilotti et al., AAPG Bulletin, V. 84, No. 6 (June 2000), p. 727 - 740.)
How to use:
This program, mapping.gawk, is written in gawk so that it can be read on multiple platforms. The spreadsheet sandd.xls is an Excel spreadsheet tested on a Macintosh and Windows. I have no plans to port these programs to other platforms.
The gawk program reads a tabulated ascii file that has the following information:

X, Y, Z for three points of a surface (say a geologic contact on a map)

and computes the strike and dip of the surface defined by those three points.

The Excel spreadsheet (screenshot below) has two sheets and allows the user to input 3 points on a surface and calculates the strike and dip of that surface.

get mapping.gawk (control - click)
get the Excel spreadsheet (control - click)
excel_sheet
Contact me: guzofski@fas.harvard.edu with any questions or problems. However, use at your own risk- I claim no responsibility for problems or errors caused by the software or the data they generate.
12/20/2002