comp.lang.idl-pvwave archive
Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado

Home » Public Forums » archive » map_set stereographic projection
Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend 
Switch to threaded view of this topic Create a new topic Submit Reply
map_set stereographic projection [message #48520] Thu, 27 April 2006 05:50 Go to next message
dvila is currently offline  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 Go to previous message
dvila is currently offline  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 #48572 is a reply to message #48520] Fri, 28 April 2006 21:18 Go to previous message
James Kuyper is currently offline  James Kuyper
Messages: 425
Registered: March 2000
Senior Member
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 Go to previous message
David Fanning is currently offline  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 Go to previous message
James Kuyper is currently offline  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 Go to previous message
David Fanning is currently offline  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 Go to previous message
mattie is currently offline  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 Go to previous message
James Kuyper is currently offline  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.
  Switch to threaded view of this topic Create a new topic Submit Reply
Previous Topic: Re: How to rotate a figure in PLOT?
Next Topic: Upgrade to IDL 6.3 breaks IDLtoAVI

-=] Back to Top [=-
[ Syndicate this forum (XML) ] [ RSS ] [ PDF ]

Current Time: Wed Oct 08 13:35:46 PDT 2025

Total time taken to generate the page: 0.00709 seconds