Steve Mauget's recent email, along with one from yesterday, revisit
the lon-lat <--> UTM conversion issue. Six months ago I submitted
a similar query to this mailing list, got a number of helpful suggestions,
and then posted the appended message with solutions. Sorry for the
re-run, but it seems like an issue which crops up from time-to-time.
Also, the full reference to the Snyder publication is:
Snyder, J.P., Map Projections -- A Working Manual, USGS Professional
Paper 1395, 1987.
This is a very useful reference, containing formulas (and tables for
checking implementations of these formulas) and example charts.
Dan
...................................................................
Dan Scheirer Email: scheirer@emma.geo.brown.edu
Department of Geological Sciences Office: 401-863-7573
Box 1846 Lab: 401-863-1701
Brown University Fax: 401-863-2058
Providence, RI 02912
################## Originally posted 26 Nov 1997 ############
Hi, GMT-folks
I raised the question of conversions between UTM meters and lon-lat
a week or so ago, got a helpful reply, and had a number of requests to
pass on any solutions which might come my way. Here goes...
1) converting UTM meters <--> lon-lat
Andy Goodliffe pointed out a straightforward way of doing this with
mapproject. Here is an example of some lines to convert a single lat-lon
pair to the equivalent UTM meter position (eastings and northings):
set REF_ELLIPSOID = WGS-84
set UTM_ZONE = 25
set LAT = 37.331 ; set LON = -32.266
gmtset D_FORMAT %12.5f
gmtset ELLIPSOID $REF_ELLIPSOID
echo $LON $LAT $LON $LAT | \
mapproject -R0/1/0/1 -F -C -Ju$UTM_ZONE/1:1
It's probably a good idea to explicitely specify the reference ellipsoid
(GMT defaults to WGS-84 but the UTM projection as defined, is based on
Clarke-1866). You also need to specify the UTM zone. I change the
D_FORMAT gmtdefault so it doesn't truncate some useful digits. I checked
out a couple of these conversions using examples from John Snyder's book
(USGS, 1987?), and my results (UTM meters) agreed to better than 1 ppm
with his. Thus, I think this conversion can be viewed as exact.
The way to go from UTM meters to Lon-lat's is to use exactly the same
command but add the -I (inverse) option.
I'm pretty sure that when using mapproject with -C, the -R region which
is specified really doesn't make any difference; hence, I just used a
generic region.
The extension of these conversions to many lon-lat or XY pairs is
straightforward.
2) Superimposing Lon-lat ticks/lines on UTM XY ticks/lines.
My goal is to plot grids with both lon-lat and UTM XY ticks/lines/annotations
superimposed.
I am able to do this fairly simply when my grid is in UTM meters. Basically,
just plot the grid and whatever UTM meter basemap on top of it. Then,
figure out the lon-lat values of the lower-left and upper-right positions of the
UTM meter basemap. Then, plot a new basemap on top of this, using
the -JU$UTM_ZONE/$XSIZE specification (where $XSIZE is the width of the
UTM meter basemap) and the -R(lower-left/upper-right)r lon-lat specification.
This -R...r makes the outline of the basemap from this call rectangular;
individual lon or lat grid lines are not exactly vertical or horizontal.
Actually, the above solution *should* work but doesn't quite. I think there
may be a bug in the -JU psbasemap specification. When I tried -JU on a
number of platforms, I got sluggish performance and invalid postscript.
When I used -JT with the central meridian appropriate for the $UTM_ZONE,
the lon-lat psbasemap came out great. I've reported this behavior as
a possible bug to the GMT developers.
If one has a lon-lat grid, I think you'll need to plot it with -JU (or -JT)
before you can overlay UTM meter lines.
3) When working with UTM meter annotations (which range in magnitude of
10**5 -- 10**6 meters, I found the perl script which Steven Howell posted
to this mailing list in August of 1996 helpful to remove the ensuing
exponential annotations for these large numbers. Appended at the bottom
of this message is his version which has been useful for me.
Regards,
Dan
##############################################################################
Here is a short (4 lines of code) perl script that might work. If applied
to a postscript file, it searches for text strings containing a single
number in exponential format and replaces them with the same number in
decimal form.
I tried it using GMT 3.0 and Perl 5.003 on a NeXT running system 3.2.
#--------cut here for the perl script. I called it refexp.p--------
#!/usr/local/bin/perl
#You might have to change the path
#
#Replaces postscript strings containing single numbers in exponential
#format with an integer representation.
$expnum = '\(([+-]?\d\.?\d*[eE][+-]?\d{1,2})\)';
#Postscript strings are encased in parentheses. This expression matches a
#'(' followed immediately by an optional + or minus, a single digit
#optionally followed by a period, then zero or more digits, E or e,
#another optional + or -, and one or two digits, then a closing ')'.
while(<>) {
s/$expnum/sprintf("(%d)",$1*1)/eg;
#replace the matched expression (multiplied by 1 to convince perl it's
#a number) with an integer encased in parentheses.
print $_;
}
#-----------------end of the perl script-------------------------------
I tested it with a short GMT script:
#------------Cut here for the csh script-----------------------------
#!/bin/csh
psxy -JX5/5 -R0/5000000/0/10000000 -P -Ba2000000 > test.ps <<END
1 1
200000 7000000
4000000 3000000
END
refexp.p test.ps > newtest.ps
#--------------------------end of csh-----------------------
Good luck,
Steven Howell