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
|