sso@nilu.no (Sverre Solberg) wrote in message news:<55f39a3a.0411110650.63c260ce@posting.google.com>...
> 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
>
Hi Sverre,
I don't think this is your fault.
If you comment out your plotting of the grid, and instead use IDL's
own Map_grid wih the Latdel and Londel keywords, then the same effect
occurs. (See amendment to code below). So it's not just your coding at
fault.
And if your remove the /trans, then the grid lines are wonderfully
straight,
though rotated. I suspect that there is an inherent problem in IDL's
mapping routines in the way they handle Transverse Mercator and
rotation.
Send your code to support@rsinc.com and say "What gives?" Or better
still,
search RSI's Tech Tips first.
Cheers,
Andrew
dot
Cool
at
dsto
dot
defence
dot
gov
dot
au
Adelaide, South Oz
>
> -----------------------------------------
>
> 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
>
lat_range = lat1 - lat0
lon_range = lon1 - lon0
latdel = lat_range/float(ny)
londel = lon_range/float(nx)
map_grid,londel=londel,latdel=latdel,glinestyle=0
> ;..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
