| Re: create an UTM grid [message #78409] |
Mon, 28 November 2011 14:35  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
nata writes:
> UTM is a cartesian coordinate system. Theoretically, the distance
> between two points on the grid could be performed using the
> Pythagorean theorem (because all distances should be equal). What I
> realize is that the distances are different and never equal to 1km. I
> checked that using the routine MAP_2POINTS. Do you know why ?
I think mostly because you are not following my advice
exactly. ;-)
In my grids, the distances are *exactly* 1000 meters apart.
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 #78410 is a reply to message #78409] |
Mon, 28 November 2011 14:33  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
nata writes:
> Thank you David. Just to be sure, if I want a grid of 1000x1000, 1km
> resolution and centered on 45.42419lat -73.9372lon, is the following
> code right ?
>
> center_lat=45.42419
> center_lon=-73.9372
> xdim=1000
> ydim=1000
> 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)
>
> x_grid=xycenter[0] + FINDGEN(xdim)*resolution - xdim/2.*resolution
> y_grid=xycenter[1] + FINDGEN(ydim)*resolution - ydim/2.*resolution
>
> result = MAP_PROJ_INVERSE(x_grid, y_grid, MAP_STRUCTURE=map_utm)
>
> lon=REFORM(result[0,*])
> lat=REFORM(result[1,*])
Well, something like that. :-)
More like this, I think. (Maybe this is what you had.)
mapStruct = Map_Proj_Init('UTM', /GCTP, $
Center_Lon=-74, Center_Lat=45, Ellipsoid=24)
xycenter = Map_Proj_Forward(-74, 45, Map_Structure=mapStruct)
xstart = [xycenter[0]-49995
ystart = [xycenter[1]-49995
xgrid = Findgen(1000)*1000 + xstart
ygrid = Findgen(1000)*1000 + ystart
xdim = 1000
ydim = 1000
xgrid = Rebin(xgrid, xdim, ydim)
ygrid = Rebin(Reform(ygrid, 1, ydim), xdim, ydim)
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 #78411 is a reply to message #78410] |
Mon, 28 November 2011 12:02 |
natha
Messages: 482
Registered: October 2007
|
Senior Member |
|
|
I think the following code is better :
x_grid=REBIN(x_grid,xdim,ydim,/SAMPLE)
y_grid=TRANSPOSE(REBIN(y_grid,ydim,xdim,/SAMPLE))
result=MAP_PROJ_INVERSE(x_grid, y_grid, MAP_STRUCTURE=map_utm)
res_lon=REFORM(result[0,*],xdim,ydim)
res_lat=REFORM(result[1,*],xdim,ydim)
By the way, I am missing something here...
UTM is a cartesian coordinate system. Theoretically, the distance
between two points on the grid could be performed using the
Pythagorean theorem (because all distances should be equal). What I
realize is that the distances are different and never equal to 1km. I
checked that using the routine MAP_2POINTS. Do you know why ?
|
|
|
|
| Re: create an UTM grid [message #78413 is a reply to message #78411] |
Mon, 28 November 2011 11:25 |
natha
Messages: 482
Registered: October 2007
|
Senior Member |
|
|
Thank you David. Just to be sure, if I want a grid of 1000x1000, 1km
resolution and centered on 45.42419lat -73.9372lon, is the following
code right ?
center_lat=45.42419
center_lon=-73.9372
xdim=1000
ydim=1000
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)
x_grid=xycenter[0] + FINDGEN(xdim)*resolution - xdim/2.*resolution
y_grid=xycenter[1] + FINDGEN(ydim)*resolution - ydim/2.*resolution
result = MAP_PROJ_INVERSE(x_grid, y_grid, MAP_STRUCTURE=map_utm)
lon=REFORM(result[0,*])
lat=REFORM(result[1,*])
Thank you again ! I really appreciate your support !
nata
|
|
|
|
| Re: create an UTM grid [message #78414 is a reply to message #78413] |
Mon, 28 November 2011 10:12 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
nata writes:
> 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,
This is almost right. When you set up a UTM projection,
the CENTER_LONGITUDE and CENTER_LATITUDE keywords are
actually used to choose the correct UTM zone, not really
to center the map projection in that zone. To center the
grid, you will have to project your center latitude and
longitude into XY coordinates, then create your grid:
mapStruct = Map_Proj_Init('UTM', /GCTP, $
Center_Lon=-74, Center_Lat=45, Ellipsoid=24)
xycenter = Map_Proj_Forward(-74, 45, Map_Structure=mapStruct)
xrange = [xycenter[0]-250000, xycenter[0]+250000]
yrange = [xycenter[1]-250000, xycenter[1]+250000]
ll = Map_Proj_Inverse(xrange, yrange, Map_Structure=mapStruct)
Print, 'Longitude Range: ', ll[0,0], ll[0,1]
Print, 'Latitude Range: ', ll[1,0], ll[1,1]
The actual grid will be determined by how many "cells"
you want to have in your data range. In other words,
how big do you want each cell of your grid to be? Or,
what is the resolution of your grid? If you
want each cell to be, say 5000 meters on a side, then
your grid would be a 100x100 array.
IDL> Print, (xrange[1] - xrange[0]) / 100
5000.000
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 #78558 is a reply to message #78410] |
Mon, 28 November 2011 14:38 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> xstart = [xycenter[0]-49995
> ystart = [xycenter[1]-49995
Sorry, this should have been 50000-500=49500:
xstart = [xycenter[0]-49500
ystart = [xycenter[1]-49500
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.")
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
