Comments on Map Projections in the Vicinity of the SAFOD drill site

This page summarizes the steps that I have been using to convert from Latitudes and Longitudes to a local coordinate system for the PASO dataset in the vicinity of the SAFOD drill site. It has been put together in reponse to a request at the Fall 2003 AGU meeting by members of the SAFOD working group to come up with a uniform coordinate system for investigators working in the area.

The area of interest around SAFOD is small enough (no more than a few 10's of km) that a "flat earth" representation is adequate. One would nevertheless seek to minimize the distortion inherent in mapping an ellipsoidal surface (i.e., latitude, longitude, depth) onto a plane (i.e., cartesian X, Y, Z). Based on its minimal distortion capabilities, I use the Universal Transverse Mercator (UTM) projection for this purpose. There are three reference frames to consider:

Ellipsoidal Reference Frame.

There are at least 24 existing definitions of the ellipsoidal shape of the Earth. The software I used to convert from Lat/Lon to UTM is general enough to accomodate any of these frames, but because we now almost always use GPS coordinates, I default to the the WGS-84 frame.

UTM Reference Frame.

There are a large number of UTM reference frames developed for different parts of the world. The two main ones we need to worry about here are the North American Datum 1927 (NAD27) and North American Datum 1983 (NAD83) frames. The NAD27 frame is used by the USGS in producing their topo maps, and for this reason is my preferred UTM projection.

Local Reference Frame.

It will be useful for many applications to specify a reference frame local to the drill site area, as distance from a given point in a given direction. In the past I have used the 1966 Parkfield epicenter and mean sea level (MSL) as the coordinate center, with a rotation of 138.5 degrees clockwise to align the Y axis more or less with the local strike of the SAF. The Y axis is positive to the SE, the X axis is positive to the NE and Z is positive down.

I employ a few algorithms sequentially to convert from Lat/Lon to Local frames (and reverse). The procedure is obviously not as elegant as it could be, but basically I just stopped worrying about this when I got something to work. Please feel free to do additional programming.

There are 3 steps involved:

1. Conversion from Lat/Lon to NAD83

2. Conversion from NAD83 to NAD27

3. Conversion from NAD27 to the Local Frame

Conversion from Lat/Lon to NAD83 (and back again):lat2utm and utm2lat

lat2utm and utm2lat are simple subroutine drivers that perform the conversion from Lat/Lon to UTM NAD83 (and back again). They are easy to modify to fit your particular formatting needs. The subroutines (in utmsubs.f) are my FORTRAN adaptions of a nice bit of C++ code written by Chuck Gantz of the Defence Mapping Agency.

Here are links to the source code:

lat2utm.f

utm2lat.f

utmsubs.f

Here are the Compilation steps:

f77 -c utmsubs.f
f77 lat2utm.f utmsubs.o -o lat2utm
f77 utm2lat.f utmsubs.o -o utm2lat

Here's an example of running lat2utm. The stuff that you type in is in italics:

gretchen.geo.rpi.edu% lat2utm
Enter a lat, lon pair: 35.970 -120.5225
Starting position(Lat, Long): 35.970000000000 -120.52250000000
Calculated UTM position(Northing, Easting, Zone):
3983458.9545501 723397.49328311 10S
Calculated Lat, Long position(Lat, Long): 35.970000000842 -120.52249999903
gretchen.geo.rpi.edu%


Comments:

I input the Lat/Lon of interest (35.970 -120.5225) and am returned the UTM83 Northing and Easting (3983458.9545501, 723397.49328311). Note that the units are Meters. I also get the UTM Zone (S) and Zone Number (10) for this region (10S will apply everywhere in the vicinity of the SAFOD site). The program then does the inverse operation as a check to see if the procedure is working.

Note that this software is hardwired for WGS-84.

Here's an example of utm2lat. Again, the stuff that you type in is in italics:

gretchen.geo.rpi.edu% utm2lat
Enter a UTM (N,E) pair: 3983458.9545501 723397.49328311
Starting position(UTMNorthing, UTMEasting): 3983458.9545501 723397.49328311
Enter Zone Number (e.g. 18 for New England):
10
Enter UTM Zone (e.g. T for New England):
S
Calculated Lat, Long position(Lat, Long): 35.970000000841 -120.52249999903
Calculated UTM position(Northing, Easting, Zone):
3983458.9546457 723397.49336818 10S
gretchen.geo.rpi.edu%


Comments:

Note that I had to input the appropriate UTM Zone and Zone Number (10S) for this to work. Of course, you can hardwire this in your own application.

nad2nad is part of a larger cartographic projection package called PROJ.4 written by Gerald I. Evenden of the USGS. You can obtain the latest version of PROJ.4 from http://remotesensing.org/proj. I have a older version that I can pass on to you if you wish. After you compile this package (not difficult - just follow the README in the distribution), you can (for example) use the following syntax to convert from NAD83 to NAD27:

nad2nad -i 83,utm=10 -o 27,utm=10 -r conus < park_utm83 > park_utm27 

Translation:

-i 83,utm=10	means the input is NAD83, UTM Zone 10
-o 27,utm=10	means the output is NAD27, UTM Zone 10
-r conus	means use the conversion table appropriate for the conterminous 48 States of the US.
< park_utm83	means redirect the standard input from the keyboard (or whatever) to the file park_utm83
> park_utm27	means redirect the standard output from the screen (or whatever) to the file park_utm27

Note that the order of input and output is (Easting, Northing)

The reverse operation (NAD72 to NAD83) is:

nad2nad -i 27,utm=10 -o 83,utm=10 -r conus < park_utm27 > park_utm83

If you desire instant gratification use a "-" after the "conus" and input the Easting and Northing. For example, to convert the 1966 Parkfield Epicenter directly (you type the italics):

nad2nad -i 83,utm=10 -o 27,utm=10 -r conus -
723397.49328311    3983458.9545501
723492.12    3983265.09
[CNTRL-D] <- end of input


or going backwards:
nad2nad -i 27,utm=10 -o 83,utm=10 -r conus -
723492.12    3983265.09
723397.49    3983458.96
[CNTRL-D]


The redirected input file will be a series of lines of Eastings and Northings.

Once we have decided on a reference location and strike, the conversion from NAD27 to the local frame is simply an affine translation followed by a rotation. Here's an example of a simple FORTRAN code to convert from NAD27 a local frame based on the 1966 Parkfield Epicenter (723492.12, 3983265.09) with a rotation of 138.5 degrees clockwise. Again, remember that the Y axis (utmnr) will be positive to the SE, and the X axis (utmer) positive to the NE (that's the reason for the minus sign near the end). Also, the output units are kilometers.

c---Origin and rotation
        utmno = 3983265.09
        utmeo = 723492.12
        az = 138.5

        azr = az*3.14159/180.
        cosaz = cos(azr)
        sinaz = sin(azr)
        call openold(' Enter name of NAD27 file: ',4)
        open(10,file='temp.local')

2       read(4,*,end=60) utme, utmn

        de = utme - utmeo
        dn = utmn - utmno

        utmer = de*cosaz - dn*sinaz
        utmnr = de*sinaz + dn*cosaz

        utmer = -utmer/1000.
        utmnr =  utmnr/1000.

        write(10,301) utmer, utmnr
301     format(2f12.5)

        go to 2
60      stop
        end

The reverse steps would substitute the following lines in the above:

        dy = dy*1000.0
        dx = -dx*1000.0

        de =  dx*cosaz + dy*sinaz
        dn = -dx*sinaz + dy*cosaz

        utme = de + utmeo
        utmn = dn + utmno

where dx and dy are the locally defined X and Y corrdinates (= utmer and utmnr above).

The SAFOD Pilot Hole is located at:

Latitude:

35 deg 58 min 27.3286 sec

35.97425794 deg

Longitude:

-120 deg 33 min 07.5857 sec

-120.55210714 deg

NAD83 coordinates from lat2utm:

3983863.9008305    720715.42853055    10S

NAD27 coordinates from nad2nad:

3983669.99    720810.15

Position relative to 1966 Parkfield Epicenter frame (in km):

x = -1.74045

y = -2.07929

z = -0.6675

There are two ways of specifying an elevation that we need to discuss: ellipsoidal elevation and orthometric elevation.

Elliposidal elevation is defined as the displacement normal to the surface of a reference ellipsoid that approximates the shape of the Earth. This is the elevation provided by GPS, and typically the reference ellipsoid used is WGS84.

Orthometric elevation is defined as the vertical displacement from a reference gravitational equipotential surface known as the Geoid. This is what we usually mean by elevation w.r.t. Mean Sea Level (MSL) since the oceans in some sense are used to define this reference surface. This is the elevation that appears on USGS topo maps.

An exact transformation between Orthometric and Ellipsoidal elevations is a bit tricky because of the difference in how "vertical" is defined (the former uses the local gravity field and the latter uses a simple geometric figure). However, we can get pretty close (to within a meter or so) by simply adding the local displacement of the Geoid w.r.t. the Ellipsoid.

For example, if my GPS reciever tells me that my elevation is 600 meters, and I know that at the particular latitude and longitude of my measurement the Geoid is 30 meters below the ellipsoid, then my Orthometric elevation is (to a reasonable estimation) 630 meters, or in other words I can say that I am 630 meters above MSL.

Recommendations for SAFOD elevations.

Because the area around Middle Mountain has appreciable relief, getting an accurate elevation from a USGS topo map can be dicey. Not only do you have to know where you are exactly on the side of that mountain, but the contour intervals on the topo maps are on the order of 20-40 feet. Thus, I would recommend using the USGS topo maps only for corroboration (i.e., is your reported elevation reasonable?).

A single GPS receiver observing for a short amount of time will not be able to determine elevations to accuracies better than several 10's of meters. Don't use elevations determined this way.

Two or more GPS receivers can be used to determine differential elevations (with some post-processing) and these can be quite good - cm range accuracies are normal. Note that they still must be tied to some reference (i.e. a measurement must be taken where the elevation is known - but remember if you occupy a bench mark that the elevation on the mark is MSL, not ellipsoidal). You might also get down to a meter of accuracy if you occupy the site for a long period of time with a single receiver (long enough for climate to average itself out - probably several days). Another possibility is that your GPS receiver is capable of receiving and processing corrections from regional radio transmitters (this is what the survey community often calls "Differential GPS" - conflicting with the UNAVCO usage of this term). Elevations thus determined are often accurate at the meter level.

Thus, if you have differential, or otherwise accurate GPS estimates of elevation as described above, you most likely have a better idea about your 3D position than you can get from other sources. In this case you should look up the elevation of the Geoid with respect to the ellipsoid and add this to your GPS elevation to get elevation above MSL.

If you want to be precise, use one of the recent NGS geoid models like GEOID03 or USGG2003. NGS has an online utility to look up the Geoid elevation or you can download the program "intg" from their site if you want to do it yourself. Start at the web site:

http://www.ngs.noaa.gov/GEOID

and click on the Geoid model you want to use. There will be a link to an "Interactive Computation" from these pages.

Below is a plot of the Geoid elevation created with intg.f using the GEOID03 and USGG2003 geoids. Two things to note: first, the difference between geoids amounts to not much more than 10 cm anywhere in our region, and second, the geoid elevation varies by about 2 meters over our region, so if you add, say, 33.2 m to your GPS elevations you'll be within a meter of the appropriate Geoidal correction, and that probably will be good enough (considering all other sources of error).

If you don't have accurate GPS elevations, then I recommend that you use a DEM (Digital Elevation Model). The USGS now has a very convenient setup for obtaining elevation models from the National Elevation Dataset (NED). Start at the web site:

http://seamless.usgs.gov/

There is a link from this page to one that explains what NED is all about. One very nice thing about this new site is that they have transitioned to WGS84 (all their old DEMs were in NAD27), so you can look up your elevations directly using your GPS derived coordinates.

Click on the bar that says "View and Order Data Sets - United States Viewer" to put in a data request. There are a number of ways to make a request - one that I've used successfully was to select the "Define Area by Coordinates" link and put in Lat and Long delimiters for the SAFOD region.

Some insights:

1. Unless you have ArcInfo (I don't), be sure to request GridFloat data. This will give you a binary file with 4 byte reals for elevation. You will also get a header file (*HDR) that tells what's in the binary file, plus a few other files that I haven't bothered with yet.

2. The binary file is Little-Endian, which is fine if you are on a PC but not fine if you are on a Sun or Mac. You will need to swap bytes at some stage before using this data.

3. Use the NED data and NOT the SRTM 30 meter data. The fuzzing of the SRTM dataset to this level of resolution makes it less useful.

In playing with this data, I wrote a program (interp_ned.f) that does a bicubic interpolation over the data in the binary file and allows the user to estimate an elevation at a specific latitude and longitude. This program will also do byte swapping for you if required. You are welcome to use it but as with any freebie software like this there is no warranty.

This page was last updated on 10 February, 2003 by Steven Roecker
roecks@rpi.edu