| Re: map_set stereographic projection [message #48493] |
Fri, 28 April 2006 07:24  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
kuyper@wizard.net writes:
> stereo = MAP_PROJ_INIT('Stereographic', CENTER_LONGITUDE=-105, $
> CENTER_LATITUDE=90)
>
> etc.
Ah,*very* helpful!! I had to add a POSITION=[0,0,1,1] to your
MAP_SET command, but now I have the image in the correct position
with respect to the map. This whole business with the UV coordinates
was *exactly* what I wanted to know about. Thanks for this!
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
|
|
|
|
| Re: map_set stereographic projection [message #48494 is a reply to message #48493] |
Fri, 28 April 2006 06:40  |
James Kuyper
Messages: 425
Registered: March 2000
|
Senior Member |
|
|
dvila wrote:
> Hi all,
>
> I'm trying to deal with a pre-projected polar stereographic image with
> this geometrical characteristics:
>
> # It is a 1121x881 polar stereographic grid.
> # Point (1,1) is at 23.117N 119.023W.
> # Point (1,881) is at 53.509N 134.039W.
> # Point (1121,1) is at 19.805N 80.750W.
> # Point (1121,881) is at 45.619N 59.959W.
> # The y-axis is parallel to 105W.
> # The resolution is 4.7625km at 60N.
> # The pole point is (I,J) = (400.5,1600.5)
>
> I don't know how may I set the map_set routine to fit the map with the
> data. Is it possible to work with pre-projected data?
stereo = MAP_PROJ_INIT('Stereographic', CENTER_LONGITUDE=-105, $
CENTER_LATITUDE=90)
; The repeat at the end will serve to close the box in oplot command
longitude = [-119.023D, -134.039, -59.959, -80.750, -119.023]
latitude = [23.117D, 53.509, 45.619, 19.8057, 23.117]
uv = MAP_PROJ_FORWARD(longitude,latitude,MAP_STRUCTURE=stereo)
; These are the u-v values corresponding to pixel centers along each
edge.
topv = (uv[1,1]+uv[1,2])*0.5
botv = (uv[1,0]+uv[1,3])*0.5
leftu = (uv[0,0]+uv[0,1])*0.5
rightu = (uv[0,2]+uv[0,3])*0.5
xscale = (rightu-leftu)/(1121-1)
yscale = (topv-botv)/(881-1)
; U-V coordinates of midpoints of outer edges
u = [leftu-0.5*xscale, 0.5*(leftu+rightu), rightu+0.5*xscale, $
0.5*(leftu+rightu)]
v = [0.5*(botv+topv), topv+0.5*yscale, 0.5*(botv+topv),
botv-0.5*yscale]
lonlat = MAP_PROJ_INVERSE(u, v, MAP_STRUCTURE=stereo)
limit = [lonlat[1,*],lonlat[0,*]]
WINDOW,XSIZE=1120,YSIZE=880
MAP_SET, 90, -105, /STEREOGRAPHIC,
LIMIT=limit,/NOBORDER,XMARGIN=0,YMARGIN=0
; To verify that we've set up the map projection correctly:
pixels = CONVERT_COORD(longitude,latitude,/DATA,/TO_DEVICE)
print,pixels
data = FLTARR(1120,880)
; Fill in data
TVSCL,data
MAP_GRID,/LABEL
MAP_CONTINENTS,/HIRES
; Draw closed box through corner pixels
OPLOT,longitude,latitude
|
|
|
|
|
|
|
|
| Re: map_set stereographic projection [message #48497 is a reply to message #48496] |
Fri, 28 April 2006 04:23 |
dvila
Messages: 13
Registered: April 2006
|
Junior Member |
|
|
Thank you, Peter. The data provider offers a little Fortran program to
calculate the lat, lon of each pixel, so it's very easy to perform an
array with the lat,lon of each point of the data array. I'll download
the routines to know which are the input data and the data format.
Thanks!!
Daniel
Peter Albert wrote:
> Hi Daniel,
>
> I found it often useful to use the 8-element limit keyword with
> map_set. You are not resticted to the corners of the mapped region then
> and have more freedom to choose the "anchor points". Moreover,
> sometimes the corners aren't actually on earth, in those cases you
> _must_ use the 8-element vector. It's a bit tricky to get used to, but
> it's worth trying. Basically, you are specifying 4 points (as lat/lon
> pairs) _anywhere_ on the left, upper, right and lower boundary of the
> area. "Anywhere, how's that supposed to work?" Well, it just does.
> Funny enough.
>
> But then this is probably not the right way for your problem, as you'd
> like to get the regular lat-lon projection. What you actually need are
> the lon/lat values for each pixel. With those, you can either project
> your data in any projection (using e.g. Liam Gumleys IMAGEMAP() routine
> or our MAP_IMAGE__DEFINE object found at
> http://wew.met.fu-berlin.de/idl).
>
> If you actually want to transform the data itself into regular lon /
> lat grid, e.g. for pixelwise comparison with other datasets, you might
> want to use LONLAT2REG(), found on the same website. This routine
> averages irregularly gridded lon / lat data into any regular lon / lat
> array.
>
> But now, how do you get the lon/lat values for each pixel?
>
> First try: Ask the data provider. They should have the data. Somewhere.
> Second try: Download the "proj" software and calculate them yourself.
> (http://proj.maptools.org/). proj is _not_ easy to get started with,
> but if you have to do map transformations more than once (including,
> possibly, different ellipsoids), it's well woth every minute spent
> reading the manual.
>
> I would recommend calculating all map coordinates in map space first
> (i.e. coordinates being given in x = meters east of greenwich meridian;
> y = meters from the equator) and use proj to transform those into
> lat/lon. I have to admit that I only hat to deal with sinusoidal
> projection so far, which was pretty easy, and do not know how to
> actually approach polar stereographic. But in case you don't get the
> lat/lon data from the data providers, I would recommend to give proj a
> good try. And the proj people have a quite helpful mailing list, too...
>
>
> Best regards,
>
> Peter
|
|
|
|
| Re: map_set stereographic projection [message #48501 is a reply to message #48497] |
Fri, 28 April 2006 01:33 |
peter.albert@gmx.de
Messages: 108
Registered: July 2005
|
Senior Member |
|
|
Hi Daniel,
I found it often useful to use the 8-element limit keyword with
map_set. You are not resticted to the corners of the mapped region then
and have more freedom to choose the "anchor points". Moreover,
sometimes the corners aren't actually on earth, in those cases you
_must_ use the 8-element vector. It's a bit tricky to get used to, but
it's worth trying. Basically, you are specifying 4 points (as lat/lon
pairs) _anywhere_ on the left, upper, right and lower boundary of the
area. "Anywhere, how's that supposed to work?" Well, it just does.
Funny enough.
But then this is probably not the right way for your problem, as you'd
like to get the regular lat-lon projection. What you actually need are
the lon/lat values for each pixel. With those, you can either project
your data in any projection (using e.g. Liam Gumleys IMAGEMAP() routine
or our MAP_IMAGE__DEFINE object found at
http://wew.met.fu-berlin.de/idl).
If you actually want to transform the data itself into regular lon /
lat grid, e.g. for pixelwise comparison with other datasets, you might
want to use LONLAT2REG(), found on the same website. This routine
averages irregularly gridded lon / lat data into any regular lon / lat
array.
But now, how do you get the lon/lat values for each pixel?
First try: Ask the data provider. They should have the data. Somewhere.
Second try: Download the "proj" software and calculate them yourself.
(http://proj.maptools.org/). proj is _not_ easy to get started with,
but if you have to do map transformations more than once (including,
possibly, different ellipsoids), it's well woth every minute spent
reading the manual.
I would recommend calculating all map coordinates in map space first
(i.e. coordinates being given in x = meters east of greenwich meridian;
y = meters from the equator) and use proj to transform those into
lat/lon. I have to admit that I only hat to deal with sinusoidal
projection so far, which was pretty easy, and do not know how to
actually approach polar stereographic. But in case you don't get the
lat/lon data from the data providers, I would recommend to give proj a
good try. And the proj people have a quite helpful mailing list, too...
Best regards,
Peter
|
|
|
|
|
|
|
|
| Re: map_set stereographic projection [message #48507 is a reply to message #48506] |
Thu, 27 April 2006 11:51 |
dvila
Messages: 13
Registered: April 2006
|
Junior Member |
|
|
David,
Certainly I started with that link but I couldn't deal with the limit
keyword on map_set routine because the four corners of the image have
different latitude and longitude (due to the projection, I supposed).
Otherwise, I have all other data (to compare with these ones) in a
lat-lon regular grid projection...
You can find an image in
http://essic.umd.edu/~dvila/ftp/
It's daily precipitation data (in mm) with 9.999E20 as missing data.
Thanks again!
Daniel
David Fanning wrote:
> dvila writes:
>
>> I'm trying to deal with a pre-projected polar stereographic image with
>> this geometrical characteristics:
>>
>> # It is a 1121x881 polar stereographic grid.
>> # Point (1,1) is at 23.117N 119.023W.
>> # Point (1,881) is at 53.509N 134.039W.
>> # Point (1121,1) is at 19.805N 80.750W.
>> # Point (1121,881) is at 45.619N 59.959W.
>> # The y-axis is parallel to 105W.
>> # The resolution is 4.7625km at 60N.
>> # The pole point is (I,J) = (400.5,1600.5)
>>
>> I don't know how may I set the map_set routine to fit the map with the
>> data. Is it possible to work with pre-projected data?
>
> On occasion. See, for example:
>
> http://www.dfanning.com/map_tips/georeference.html
>
>> On the other hand, can I change the projection to a regular lat-lon
>> grid in IDL?
>
> Probably not. :-)
>
> I don't know the definitive answers to these questions, but
> I have a reason to find out. Where can I find one of these
> images you are trying to fit a map projection to?
>
> Cheers,
>
> David
> --
> David Fanning, Ph.D.
> Fanning Software Consulting, Inc.
> Coyote's Guide to IDL Programming: http://www.dfanning.com/
|
|
|
|
| Re: map_set stereographic projection [message #48516 is a reply to message #48507] |
Thu, 27 April 2006 07:23 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
dvila writes:
> I'm trying to deal with a pre-projected polar stereographic image with
> this geometrical characteristics:
>
> # It is a 1121x881 polar stereographic grid.
> # Point (1,1) is at 23.117N 119.023W.
> # Point (1,881) is at 53.509N 134.039W.
> # Point (1121,1) is at 19.805N 80.750W.
> # Point (1121,881) is at 45.619N 59.959W.
> # The y-axis is parallel to 105W.
> # The resolution is 4.7625km at 60N.
> # The pole point is (I,J) = (400.5,1600.5)
>
> I don't know how may I set the map_set routine to fit the map with the
> data. Is it possible to work with pre-projected data?
On occasion. See, for example:
http://www.dfanning.com/map_tips/georeference.html
> On the other hand, can I change the projection to a regular lat-lon
> grid in IDL?
Probably not. :-)
I don't know the definitive answers to these questions, but
I have a reason to find out. Where can I find one of these
images you are trying to fit a map projection to?
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
|
|
|
|
| Re: map_set stereographic projection [message #48553 is a reply to message #48495] |
Tue, 02 May 2006 00:56 |
peter.albert@gmx.de
Messages: 108
Registered: July 2005
|
Senior Member |
|
|
Humm, well, I guess those !quiet=0 lines are remnants from the times
when I did not know about the SILENT keyword to LOADCT and at some
point did not want to read those messages any more ... :-) I have to
admit that I completely left IDL programming last December, but a very
quick check on my old PC seems to show that no real message gets
surpressed. I'll see whether I can delete those lines on the web page.
Cheers,
Peter
|
|
|
|
| Re: map_set stereographic projection [message #48589 is a reply to message #48494] |
Fri, 28 April 2006 09:58 |
mattie
Messages: 4
Registered: January 2006
|
Junior Member |
|
|
kuyper@wizard.net writes:
> dvila wrote:
>> Hi all,
>>
>> I'm trying to deal with a pre-projected polar stereographic image with
>> this geometrical characteristics:
>>
>> # It is a 1121x881 polar stereographic grid.
>> # Point (1,1) is at 23.117N 119.023W.
>> # Point (1,881) is at 53.509N 134.039W.
>> # Point (1121,1) is at 19.805N 80.750W.
>> # Point (1121,881) is at 45.619N 59.959W.
>> # The y-axis is parallel to 105W.
>> # The resolution is 4.7625km at 60N.
>> # The pole point is (I,J) = (400.5,1600.5)
>>
> ; These are the u-v values corresponding to pixel centers along each
> edge.
> topv = (uv[1,1]+uv[1,2])*0.5
> botv = (uv[1,0]+uv[1,3])*0.5
> leftu = (uv[0,0]+uv[0,1])*0.5
> rightu = (uv[0,2]+uv[0,3])*0.5
Aren't these the u-v values corresponding to the _center_ of the
gridcell along each edge? (assuming the initial data gave the
centerpoints of the gridcell.)
>
> xscale = (rightu-leftu)/(1121-1)
> yscale = (topv-botv)/(881-1)
>
> ; U-V coordinates of midpoints of outer edges
> u = [leftu-0.5*xscale, 0.5*(leftu+rightu), rightu+0.5*xscale, $
> 0.5*(leftu+rightu)]
> v = [0.5*(botv+topv), topv+0.5*yscale, 0.5*(botv+topv),
> botv-0.5*yscale]
> lonlat = MAP_PROJ_INVERSE(u, v, MAP_STRUCTURE=stereo)
> limit = [lonlat[1,*],lonlat[0,*]]
Again, I'm not sure, but don't you have to add half a gridcell to each
direction to get the outter limit of each grid cell?
I'm asking these questions honestly because I've never quite convinced
myself which is the actual limit of the map_set command. I think it's
the outside edge of the grid cell that you want to set. So for this
example, it will be quite close (half a pixel over 800. But if you're
working with zoomed regions. 10x10 grid cells or whatever, it becomes
noticable.
But this was a very informative tutorial on how to register images when
corner points are known.
Thanks
Matt
btw. Does anyone know of a free news server that will let me post
(besides google)? The university just turned ours off. Thanks.
--
Matthew Savoie - Scientific Programmer
National Snow and Ice Data Center
(303) 735-0785 http://nsidc.org
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
