| Overlay images from WMS servers (web mapping servers) on map projections [message #47435] |
Wed, 08 February 2006 06:01  |
Jan Kristian Jensen
Messages: 10
Registered: February 2006
|
Junior Member |
|
|
I thought I had a brilliant idea of how to combine images from WMS
servers (OGC Web Mapping Services compliant servers) on map projections
using map_set or map_proj_init functions. The result is not by far good
enough - and I can't figure what's wrong.
First I use my favourite http client to send a OGC WMS getMap-request to
my favourite WMS server. Part of the URL is a boundingbox (BBox)
parameter containing coordinates for lower left and upper right corner,
and how many x-y pixels there should be. Further the projection must be
specified. An easy test case would be the equirectangular projection, or
EPSG:4623 in OGC speak.
Then the easiest thing would be to open a window with the exact same
size, call map_set with the appropriate projection and display the
image. Something like this:
--- Example - you have to provide the wms file yourself ... -
; Limit = boundingBox coordinates to WMS server
limit = [ 59.5, 3.5, 61.0, 5.2]
map_set, 0, 0, /cylindrical, xmargin = 0, ymargin = 0, limit = limit
wmsimage = read_png( filename, r, g, b)
tvlct, r, g, b
window, xsize = 1200, ysize = 1000
tv, wmsimage
; Add a nice grid
device, decomposed = 0
tek_color
map_Grid, latdel = 0.25, londel = 0.25, color = 2
-------------------------------
My wmsimage example displays the coast contour and a grid line for every
0.25 degree lat/lon. This to compare with the lines from map_grid. The
two set of grid lines are fairly close in the center of the image, but
off by up to 1.9 km at the egde of the image. It is _/*very*\_ apparent
that the wms image and the map_grid lines do not match at all exept for
a small area near the center of the image.
I may have misunderstand something about projections and/or how to set
them correctly: I think equirectangular = EPSG:4326 = cylindrical
projection with standard parallell at 0 degrees, which should be the
default projection for map_set.
The symmetry about the image center makes me suspect there is some
border issue here. If I don't use the xmargin or ymargin keywords (or
/noborder) I get an even larger error. I can thank this newsgroup for
figuring that out :)
I'm quite frustrated at the moment: The result of the above method is
NOT good enough. I can't bring my mind to accept that this method
shouldn't work - but it doesn't, and I can't figure out why....
I know you can use map_projection_init and map_projection forward to get
more fine granular control - but that really shouldn't be nescessary
with this feeble task!!! If so, Ben Tupper posted an example (Nov 15
2004) of how to do this with direct graphics. I may start there...
Cheers, Jan
--
Jan Kristian Jensen
Remove the obvious from the email adress to email me.
|
|
|
|
| Re: Overlay images from WMS servers (web mapping servers) on map projections [message #47447 is a reply to message #47435] |
Sat, 11 February 2006 00:26  |
Jan Kristian Jensen
Messages: 10
Registered: February 2006
|
Junior Member |
|
|
liamgumley@gmail.com wrote:
> I like MAP_PROJ_INIT and the related functions precisely because they
> are decoupled from the graphics device. I think the potentially weak
> link comes when you have to transfer the projection to the graphics
> device (i.e., the PLOT call).
Is this perhaps why _*ALL*_ map_proj* examples in the documentation are
using object graphics?
Extracting something useful from the documentation sometimes has more in
common with theology or astrology than science and engineering. Most of
the time, I search this group, Davids exellent book (sorry Liam, haven't
got around to buy one of yours yet!) and some triendly IDL websites to
find out how things work.
Jan
--
Jan Kristian Jensen
Remove the obvious from the email adress to email me.
|
|
|
|
|
|
|
|
| Re: Overlay images from WMS servers (web mapping servers) on map projections [message #47463 is a reply to message #47435] |
Fri, 10 February 2006 10:35 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Liam writes:
> Sorry for the delay; I got busy with some other tasks.
What!? Oh, ok. :-)
> I was initially equally frustrated with MAP_PROJ_INIT, however that
> changed when Ben Tupper posted some very good instructions on how to
> use it.
I should have known Tupper was behind this. Where is he,
anyway? He hasn't exactly been weighing in on this discussion.
He probably wants to *sell* all that knowledge!
> Since then, I've had success in using MAP_PROJ_INIT to
> establish map projections which match up with GeoTIFF satellite images.
> I've created a brief example which is available at
>
> ftp://ftp.ssec.wisc.edu/pub/gumley/IDL/geotiff/
Uh, well.... Is it my imagination, or does that shoreline not
quite line up with the shore?
I think I am beginning to see the problem. My expectations
are WAY too high! :^)
> I like MAP_PROJ_INIT and the related functions precisely because they
> are decoupled from the graphics device. I think the potentially weak
> link comes when you have to transfer the projection to the graphics
> device (i.e., the PLOT call).
I agree with this. They work GREAT in theory!
> That said, I'd love to hear other examples of how people are using
> MAP_PROJ_INIT, particularly in conjunction with georegistered satellite
> images.
Yes, me, too. I'm sure I'm going to learn a lot. :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
|
|
|
|
| Re: Overlay images from WMS servers (web mapping servers) on map projections [message #47464 is a reply to message #47435] |
Fri, 10 February 2006 10:28 |
liamgumley
Messages: 74
Registered: June 2005
|
Member |
|
|
David pointed out an error in my LIMIT definition, i.e.
limit=[46.13, -89.37, 41.52, -84.42]
which should be
limit=[41.52, -89.37, 46.13, -84.42]
If you look at the source code for MAP_PROJ_INIT, it checks the
ordering of latmin and latmax, and actually switches them if latmin >
latmax. So that explains why the LIMIT keyword was interpreted
correctly in the original code.
|
|
|
|
| Map_Proj_* Problems Re: Overlay images from WMS servers [message #47466 is a reply to message #47435] |
Fri, 10 February 2006 10:05 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
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/
|
|
|
|
| Re: Overlay images from WMS servers (web mapping servers) on map projections [message #47471 is a reply to message #47435] |
Fri, 10 February 2006 09:18 |
liamgumley
Messages: 74
Registered: June 2005
|
Member |
|
|
Sorry for the delay; I got busy with some other tasks.
I was initially equally frustrated with MAP_PROJ_INIT, however that
changed when Ben Tupper posted some very good instructions on how to
use it. Since then, I've had success in using MAP_PROJ_INIT to
establish map projections which match up with GeoTIFF satellite images.
I've created a brief example which is available at
ftp://ftp.ssec.wisc.edu/pub/gumley/IDL/geotiff/
In brief, here's how I use MAP_PROJ_INIT to overlay GSHHS coastlines on
a MODIS true color GeoTIFF image:
image = read_tiff('t1.06026.1655.LakeMichigan.143.250m.tif')
help, image
swindow, xsize=1652, ysize=2061
imdisp, image, order=1, /noresize, /noscale
map = map_proj_init('UTM', limit=[46.13, -89.37, 41.52, -84.42], $
zone=16)
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
gshhs_plot, 'gshhs_h.b', map=map, level=2
I like MAP_PROJ_INIT and the related functions precisely because they
are decoupled from the graphics device. I think the potentially weak
link comes when you have to transfer the projection to the graphics
device (i.e., the PLOT call). I have also observed some strange
behavior when the POSITION keyword is set to a subset of the graphics
window. However for cases where the projection spans the entire window,
MAP_PROJ_INIT seems to work quite well. It also ties nicely into
projection libraries such as GCTP and PROJ.4, and allows you to define
a projection the same way in IDL as you would in the external library.
That said, I'd love to hear other examples of how people are using
MAP_PROJ_INIT, particularly in conjunction with georegistered satellite
images.
Cheers,
Liam.
Practical IDL Programming
http://www.gumley.com/
|
|
|
|
| Re: Overlay images from WMS servers (web mapping servers) on map projections [message #47472 is a reply to message #47435] |
Fri, 10 February 2006 09:07 |
James Kuyper
Messages: 425
Registered: March 2000
|
Senior Member |
|
|
David Fanning wrote:
> David Fanning writes:
>
>> I absolutely cannot figure out how this projection space works.
>> The documentation, as far as I can tell, is hopeless. I presume
>> you have to look elsewhere to figure out how it works. But
>> I haven't found the right place yet.
>>
>> I wanted a map projection of the Great Lakes Region of the US,
>> positioned in the window according to a POSITION keyword. It
>> seems straightforward, but a couple of hours of work got me
>> exactly nowhere. :-(
>>
>> Can you show me how you would do this robustly? :-)
>
> Nothing!? No takers?
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?
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
