map_set stereographic projection [message #48520] |
Thu, 27 April 2006 05:50  |
dvila
Messages: 13 Registered: April 2006
|
Junior Member |
|
|
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?
On the other hand, can I change the projection to a regular lat-lon
grid in IDL?
Thanks!
Daniel Vila
ESSIC/CICS 2207 Computer & Space Building - Bldg # 224
College Park, MD 20742
|
|
|
Re: map_set stereographic projection [message #48563 is a reply to message #48520] |
Mon, 01 May 2006 08:12  |
dvila
Messages: 13 Registered: April 2006
|
Junior Member |
|
|
I agree with Kuyper, the article is very clear. Thank for the citation
and I hope that this example can be usefull for many others IDL users.
Thanks again, David
Daniel
kuyper@wizard.net wrote:
> David Fanning wrote:
>> Folks,
>>
>> I've written a short article on this topic, which you can
>> find here:
>>
>> http://www.dfanning.com/map_tips/precipmap.html
>
> Thanks for the citation - I'm glad I could help! You explained what I
> was doing more clearly than I could have, even if I'd had the time to
> do so.
|
|
|
|
Re: map_set stereographic projection [message #48573 is a reply to message #48520] |
Fri, 28 April 2006 16:06  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Folks,
I've written a short article on this topic, which you can
find here:
http://www.dfanning.com/map_tips/precipmap.html
I'm especially interested in debugging this code, so
read it quickly. :-)
In addition to the interesting topic of fitting a map
projection to an image (which I'm sure you know is just
the *opposite* of the way IDL normally does things), I
show you how to use my COLORBAR routine to display
a non-linear color scheme. (I knew it could do it,
I just didn't know how!)
Flail away!
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 #48574 is a reply to message #48520] |
Fri, 28 April 2006 15:59  |
James Kuyper
Messages: 425 Registered: March 2000
|
Senior Member |
|
|
David Fanning wrote:
> Mr. Kuyper,
>
> What can you tell us about UV coordinates? The IDL
> documentation is (as usual) silent about this point,
> simply referring to them as "XY Cartesian coordinates".
> Do we know, or care, what they *really* are?
Every map "projection", whether or not it's geometrically describable
as a projection, maps lat-lon positions on the surface of the earth to
a positions on a flat plane in space. It's conventional to use u and v
as the names of an orthogonal coordinate system describing positions on
that flat plane. A second, trivial scale-and offset mapping usually
connects positions on this plane to positions on your actual printed
map. U/V coordinates are used for the flat plane, because x and y are
reserved for the position on the printed map. In IDL terms, x and y are
device coordinates.
The relationship between lat/lon and u/v depends upon the map
projection. If you're familiar with the geometric definition of the
particular map projection you're working with, the relationship is
usually unsurprising. If you're not familiar with the map projection, I
suppose the relationship would probably be very difficult to figure
out, but I wouldn't know about that :-). What I generally do is use
MAP_PROJ to calculate a few key positions, so I can figure out what
choices they've made. For polar projections they generally use a plane
tangent to the earth at the center of the projection, with a scale in
meters and the u-v coordinates oriented toward local East and local
North respectively, at the point of tangency. For cylindrical
projections, they use a plane wrapped around the earth at the equator,
with it's center at the center point of the projection, with a scale in
meters, and U/V oriented toward local East/North at the center point.
However, I generally try to avoid building assumptions about the
scaling and orientation of the axes into my code. The V direction, in
particular, sometimes points north and sometimes south. There's a
couple of projections where they appear to use a plane tangent to a
unit sphere, rather than one tangent to the surface of the Earth, which
means that the scaling is equivalent to radians instead of meters near
the center of the projection. What I do to cope with these issues is
convert a few well-chosen positions to U/V coordinates, and use the
values for those positions to resolve any such ambiguities.
> An hour search on Goggle was equally unenlightening,
> except in referring to UV coordinates in the context
> of texture mapping, which I assume is not *really*
> how they are being used here.
>
> And please, sir, tell us how you stumbled onto using
> MAP_PROJ_* functions to solve this registration problem.
I've always been fascinated by map projections. I had a very clear
understanding of them before I ever ran into the IDL map projection
routines. My current job has required me to become familiar with the
GCTP map projection library. Before the MAP_PROJ_* routines came out, I
used MAP_SET and COORD_CONV for similar purposes, but I immediately
recognised the MAP_PROJ_* as an easier way to do things.
|
|
|
Re: map_set stereographic projection [message #48581 is a reply to message #48520] |
Fri, 28 April 2006 12:52  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Mr. Kuyper,
What can you tell us about UV coordinates? The IDL
documentation is (as usual) silent about this point,
simply referring to them as "XY Cartesian coordinates".
Do we know, or care, what they *really* are?
An hour search on Goggle was equally unenlightening,
except in referring to UV coordinates in the context
of texture mapping, which I assume is not *really*
how they are being used here.
And please, sir, tell us how you stumbled onto using
MAP_PROJ_* functions to solve this registration problem.
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 #48582 is a reply to message #48520] |
Fri, 28 April 2006 12:41  |
mattie
Messages: 4 Registered: January 2006
|
Junior Member |
|
|
kuyper@wizard.net writes:
> mattie wrote:
>> kuyper@wizard.net writes:
>>
>>> dvila wrote:
> ...
>>> ; 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.)
>
> Yes, that's what my comment line says.
Ah, so it does. I read it as the center point of each edge.
>>> ; 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?
>
> Yes, that is precisely why I wrote such things as "leftu-0.5*xscale",
Yup again. When I looked at it, I kept thinking u was upper, rather
than the
uv direction, and it's completely obvious now. I was getting thrown by
the
u[1] factor as (leftu + rightu) * .5.
So I don't have any reason to say "But" when I said:
"But this was a very informative tutorial on how to register images
when
corner points are known."
Thanks for the lesson.
Matt
--
Matthew Savoie - Scientific Programmer
National Snow and Ice Data Center
(303) 735-0785 http://nsidc.org
|
|
|
Re: map_set stereographic projection [message #48588 is a reply to message #48520] |
Fri, 28 April 2006 10:24  |
James Kuyper
Messages: 425 Registered: March 2000
|
Senior Member |
|
|
mattie wrote:
> kuyper@wizard.net writes:
>
>> dvila wrote:
...
>> ; 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.)
Yes, that's what my comment line says.
...
>> ; 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?
Yes, that is precisely why I wrote such things as "leftu-0.5*xscale",
etc.
|
|
|