| UTM mapping problems [message #41594] |
Thu, 11 November 2004 06:50 |
sso
Messages: 13
Registered: February 2002
|
Junior Member |
|
|
I have a problem with the UTM map projection which I cant find any
solution for. Perhaps someone could help? I want to plot a UTM map
for an approx 300 km area (in x and y) with a model grid (58x74 cells)
on top, and I compute the lat/lon of the grid cells. However, when
overplotting this grid the grid lines become very uneven and nothing
like a straighth line at all. Is there any way to avoid this? And why
is it behaving like this? The program below shows the problem (at
least on my screen!). I'm using IDL 5.5 (yes I know there are
upgrades) for Unix, and the 'utm_to_ll' function was taken from Ben
Tuppers UTM IDL utilities found on this newsgroup.
many thanks for any help!
Sverre Solberg
-----------------------------------------
PRO Test_utm
;..Switch to black on white:
loadct, 39
!P.color = 0
!P.background = !D.n_colors-1
;..SW corner (WGS84 system):
east = 692089.d
north = 3868229.d
dx = 5000.d
nx = 58 + 1
ny = 74 + 1
zone = 34
;..compute the lon/lat coordinates of the grid cells:
utmgrid = dblarr(nx, ny, 2)
FOR ix = 0, nx-1 DO BEGIN
FOR iy = 0, ny-1 DO BEGIN
utmx = east + ix*dx - dx/2
utmy = north + iy*dx - dx/2
IF KEYWORD_SET(utm) THEN BEGIN
utmgrid(ix, iy, 1) = utmx
utmgrid(ix, iy, 0) = utmy
ENDIF ELSE BEGIN
latlon = utm_to_ll(utmx, utmy, 'WGS84', zone = zone)
utmgrid(ix, iy, 1) = latlon(0)
utmgrid(ix, iy, 0) = latlon(1)
ENDELSE
ENDFOR
ENDFOR
lat0 = utmgrid(0, 0, 0)
lon0 = utmgrid(0, 0, 1)
lat1 = utmgrid(nx-1, ny-1, 0)
lon1 = utmgrid(nx-1, ny-1, 1)
;..Draw UTM MAP WGS84 (Aka NAD83)
; from Chuck Gantz via http://gpsy.com/gpsinfo/geotoutm.htm
; taken from the idl newsgroups
; assign constant values for WGS 84 datum for lat/lon to UTM
A = 6378137.d
eccsq = 0.00669439d
k0 = 0.9996d
map_set, 0, 0, 12.7, /trans, limit = [lat0, lon0, lat1, lon1], $
ellipsoid = [A, eccsq, k0], title = title
map_continents, /coast, /hires, thick = 2
;..Draw the grid lines:
FOR iy = 1, ny DO BEGIN
x = utmgrid(*, iy-1, 1)
y = utmgrid(*, iy-1, 0)
oplot, x, y, color = 0
ENDFOR
FOR ix = 1, nx DO BEGIN
x = utmgrid(ix-1, *, 1)
y = utmgrid(ix-1, *, 0)
oplot, x, y, color = 0
ENDFOR
END
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
