kuyper@wizard.net writes:
> Sorry - I didn't notice your question because I lost track of it in
> the middle of a discussion of someone else's quesiton..
>
> What precisely are you trying to do? More to the point, what's going
> wrong?
>
> What map projection do you want to use? If you're having trouble
> choosing a map projection, what characteristics are most important to
> you: equal area, true angles, true distance in radial, horizontal, or
> vertical directions, etc.
>
> What do you want to do with the projection once you've initialized it?
Here is the problem. I want to use the GSHHS data base on
a map:
http://www.dfanning.com/map_tips/gshhs.html
Liam has suggested a way to do this with MAP_PROJ_INIT. I
used his example at the bottom the page above. It works
well. Except...I wanted not a map of the area around Cuba,
but a map of the area around the Great Lakes.
Here is the code I used to get the Cuba map.
datafile = 'gshhs_h.b'
Window, XSize=500, YSize=350
Erase, color=fsc_color('sky blue')
map = Map_Proj_Init('UTM', $
Limit=[10.15, -95.29, 29.87, -78.77], Zone=16)
; Create a data coordinate space based on map positions.
Plot, map.uv_box[[0, 2]], map.uv_box[[1, 3]], $
Position=[0.0, 0.0, 1.0, 1.0], $
/NoData, /Isotropic, XStyle=5, YStyle=5, /NoErase
Map_GSHHS_Shoreline, datafile, /Fill, Level=3, Map_Projection=map, $
Water='dodger blue'
A couple of days ago I fooled with this for a couple of hours
and couldn't get anything that even resembled the Great Lakes.
I played around with it for a few more hours this morning, I can
now get the Great Lakes displayed (they are in the same Zone 16,
after all!), but there is other weirdness. (By the way, I think
Liam's original code had the LONMIN and LONMAX values reversed,
but *oddly* is doesn't matter! That's weird, right there.)
For example, if I do nothing but change the latitude values
in the LIMIT keyword:
map = Map_Proj_Init('UTM', $
Limit=[40.15, -95.29, 59.87, -78.77], Zone=16)
I get a representation of the Great Lakes. Fine. It doesn't
look like it is even using the limits, but fine. I won't
quibble just yet.
If I remove the ZONE=16, I get something that looks like
it is upside down and zoomed.
map = Map_Proj_Init('UTM', $
Limit=[40.15, -95.29, 59.87, -78.77)
If I add CENTER_LAT and CENTER_LON, then I get back to
something I understand.
map = Map_Proj_Init('UTM', $
Limit=[40.15, -95.29, 59.87, -78.77], $
Center_Lat=50, Center_Lon=-83)
OK. So now I think I understand a little bit of what's going on.
(It must be using the ZONE keyword to center things. Maybe that's
why it ignores LIMITS. It only centers. Maybe that's why Liam
had those numbers reversed. I don't know. I want to press ahead.)
But the map fills up my window. And I want to restrict
it to just a certain portion of my window. In Liam's code,
he has a PLOT command with POSITION=[0,0,1,1].
Plot, map.uv_box[[0, 2]], map.uv_box[[1, 3]], $
Position=[0.0, 0.0, 1.0, 1.0], $
/NoData, /Isotropic, XStyle=5, YStyle=5, /NoErase
I presume this is to set up a data coordinate space so
he can transform lat/lon values to normalized coordinates
later on to see if these polygons are in his window.
However, if I remove the ISOTROPIC keyword, my map is
distorted (which I expect), but appears zoomed, too. Why?
It seems like LIMITS *are* being used now. :-(
So, OK, I put the ISOTROPIC keyword back in. But now
I don't really want to plot into the entire window. I
want to plot into a portion of the window. I change the
POSITION coordinates to these [0.1, 0.1, 0.9, 0.9].
It doesn't make any difference. The map still fills up
the whole window. This is true even if I switch the code
in my MAP_GSHHS_SHORELINE code to be identical to the code
I used from MAP_SET. Why? It works with MAP_SET. There is
nothing about the *polygons* that prevents clipping.
So, I can ONLY get this particular map (with only God knows
what limits) to fill up my window. I can't get it in just
a portion of my window.
So, bottom line, this whole business just confuses me.
I don't readily see how to get just that portion of the
map I really want in exactly the portion of the window
I want it in. It is possible with MAP_SET, clearly.
What don't I know yet about MAP_PROJ_* to make this happen?
The documentation is NO HELP. I've been all over the
web reading up on things. I still don't understand it.
It may be "robust", but it sure ain't "straightforward". :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
