| create an UTM grid [message #78417] |
Mon, 28 November 2011 06:29  |
natha
Messages: 482
Registered: October 2007
|
Senior Member |
|
|
Hi folks,
I am trying to create an UTM grid. I never did something similar and I
got the idea from the last IDL webinar...
What I want is a grid centered in a center point and I would like to
have the lat/lon values associated to every pixel. Lets say that I
want a grid of 500x500 km centered on -74lon 45lat, is there a way to
create that grid ?
What I tried until now is something like :
map_utm = MAP_PROJ_INIT('UTM', CENTER_LONGITUDE=-74
CENTER_LATITUDE=45, ELLIPSOID=24)
llrange = MAP_PROJ_INVERSE([-500000.,500000], [-500000.,500000.],
MAP_STRUCTURE=map_utm)
Maybe it's not that simple so that's why I call your wisdom to help me
with this. Thank you in advance,
nata
|
|
|
|
| Re: create an UTM grid [message #78515 is a reply to message #78417] |
Wed, 30 November 2011 15:01  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
nata writes:
> Well, it seems that nobody knows the answer... I hope that my code is
> correct and this is another IDL weird issue.
Humm. I wrote an answer, but I don't think it has been delivered,
for some reason.
I think the values are essentially correct. Your XY projected
grid is, essentially, a piece of graph paper put down on a flat
map projection. Your great circle distance is the distance
around a sphere. I would *hope* they were different values,
in fact.
That said, plotting a row of your longitude grid shows
a disturbing regularity to the way the data varies. As
if there were a systematic error of some sort. I don't
really know what to make of it, and haven't had time
the past couple of days to investigate further.
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
|
|
| Re: create an UTM grid [message #78556 is a reply to message #78417] |
Mon, 28 November 2011 17:26 |
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
|
|
|
|
| Re: create an UTM grid [message #78656 is a reply to message #78515] |
Thu, 01 December 2011 07:36 |
natha
Messages: 482
Registered: October 2007
|
Senior Member |
|
|
Ok, thank you David. I was thinking that maybe this was something
similar to the problems you explained in "UTM Map Projection Produces
Incorrect Results".
Anyway, I will deal with my results. Thank you again,
nata
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
