comp.lang.idl-pvwave archive
Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado

Home » Public Forums » archive » create an UTM grid
Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend 
Return to the default flat view Create a new topic Submit Reply
Re: create an UTM grid [message #78556 is a reply to message #78417] Mon, 28 November 2011 17:26 Go to previous message
natha is currently offline  natha
Messages: 482
Registered: October 2007
Senior Member
I think we are doing exactly the same but you are taking in account
half pixel. You said that in your grids, the distances are exactly
1000 meters apart. Not in mines. Take a look :

center_lat= 45
center_lon=-74
xdim=100
ydim=100
resolution=1000.

map_utm=MAP_PROJ_INIT('UTM', /GCTP, CENTER_LON=center_lon,
CENTER_LAT=center_lat, ELLIPSOID=24)

xycenter=MAP_PROJ_FORWARD(center_lon, center_lat,
MAP_STRUCTURE=map_utm)

xstart=xycenter[0] - (xdim/2.-.5)*resolution
ystart=xycenter[1] - (ydim/2.-.5)*resolution

xgrid=FINDGEN(xdim)*resolution + xstart
ygrid=FINDGEN(ydim)*resolution + ystart

xgrid=REBIN(xgrid, xdim, ydim)
ygrid=REBIN(REFORM(ygrid, 1, ydim), xdim, ydim)

result=MAP_PROJ_INVERSE(xgrid, ygrid, MAP_STRUCTURE=map_utm)

res_lon=REFORM(result[0,*],xdim,ydim)
res_lat=REFORM(result[1,*],xdim,ydim)

sz=SIZE(res_lon,/DIM)

xres=FLTARR(sz)
yres=FLTARR(sz)

FOR i=0, sz[0]-2 DO FOR j=0, sz[1]-1 DO $
xres[i,j]=MAP_2POINTS(res_lon[i,j],res_lat[i,j],res_lon[i
+1,j],res_lat[i,j],/METERS,RADIUS=6378137.)
FOR i=0, sz[0]-1 DO FOR j=0, sz[1]-2 DO $

yres[i,j]=MAP_2POINTS(res_lon[i,j],res_lat[i,j],res_lon[i,j] ,res_lat[i,j
+1],/METERS,RADIUS=6378137.)

xres[sz[0]-1,*]=xres[sz[0]-2,*]
yres[*,sz[1]-1]=yres[*,sz[1]-2]

mmin_x=MIN(xres,MAX=mmax_x)
mmin_y=MIN(yres,MAX=mmax_y)

PRINT, mmin_x, mmax_x
PRINT, mmin_y, mmax_y


IDL prints:
998.289 998.727
1001.60 1002.14

Maybe this is because MAP_2POINTS do not use the same semimajor and
semiminor axis of the ellipsoid but I tried the same code using the
Sphere (6370997.0,6370997.0) and giving the same number to the RADIUS
keyword on MAP_2POINTS. The result I get is:
999.991 1000.38
999.989 1000.38

And if my grid is 1000x1000, the results are:
990.751 1000.40
990.745 1000.40

I can consider this results correct even if I was expecting all
distances to be exactly 1000m. Do you think that this is due to the
same error you explained in http://www.idlcoyote.com/map_tips/utmwrong.php
?

Anyway, 1% is not a big error but I am missing some precision here.

nata
[Message index]
 
Read Message
Read Message
Read Message
Read Message
Read Message
Previous Topic: IDL's thread pool
Next Topic: Re: IDL's thread pool

-=] Back to Top [=-
[ Syndicate this forum (XML) ] [ RSS ] [ PDF ]

Current Time: Fri Oct 10 07:35:39 PDT 2025

Total time taken to generate the page: 1.27868 seconds