| Re: map_proj_* help [message #66773] |
Tue, 09 June 2009 07:12  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> If you want to let me know where I can find these images,
> I'd be happy to have a quick look. One of my hobbies, I guess
> you would call it, is aligning map boundaries on images. ;-)
Ken,
You are probably aware by now that the code I sent you
doesn't work. :-(
I continue to work on this in my odd moments. I've
gotten to the point where I believe it is "friggin
impossible!". Which is good news, I guess. It means
we might expect a breakthrough in the next day or
two. :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
| Re: map_proj_* help [message #66791 is a reply to message #66773] |
Mon, 08 June 2009 10:28  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> You find it in the bowels of Map_Set_Init if you go though the code
> step by step and watch as it changes the parameter list it sends to
> the GCTP software.
BTW, the GCTP parameter list can be hard to find itself.
I've never seen it in any IDL documentation. Rather, I
go down the hall and take a binder off the shelf of the
resident map expert. It is in Table 6.5.4, Projection
Parameters in the HDF-EOS Library User's Guide for the
ECS Project, Volume 1: Overview and Examples, published
June 1999.
You can probably find the information elsewhere, too,
I just don't know where. :-)
Cheers,
David
--
David Fanning, Ph.D.
Coyote's Guide to IDL Programming (www.dfanning.com)
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
|
|
| Re: map_proj_* help [message #66794 is a reply to message #66793] |
Mon, 08 June 2009 10:03 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> Well, here is how I would orient this image. Trying
> to put continental boundaries on Antarctica is sketchy.
> But this will certainly give you a map data coordinate
> system that you can draw your own map features on.
>
> ;**********************************************************
> file = 'G:\data\pine_2009084_1545_modis_ch02.png'
> image = Read_PNG(file)
> image = Reverse(image, 2) ; Put the (0,0) point in UL.
>
> ; Display map in window.
> pos = [0.1, 0.1, 0.9, 0.9]
> TVImage, image, POSITION=pos, /KEEP_ASPECT, /NOINTERP, /ERASE
>
> ; Set up Map Projection. Polar Stereographic, WGS-84 Datum.
> map = Map_Proj_Init(106, DATUM='WGS 84', CENTER_LATITUDE=-90)
>
> ; Convert corners to XY coordinates. Clockwise, from UL.
> lons = [-113.594373, -98.868027, -98.568614, -117.004071]
> lats = [-71.998860, -72.492200, -76.115095, -75.491196]
> xy = Map_Proj_Forward(lons, lats, MAP_STRUCTURE=map)
>
> ; Set up map coordinate space for drawing on the map.
> Plot, [xy[0,0], xy[0,1]], [xy[1,0], xy[1,2]], $
> XStyle=5, YStyle=5, /NoData, POSITION=pos, /NOERASE
>
> ; Draw continental outlines.
> ; Map_Continents, MAP_STRUCTURE=map, /COASTS, /HIRES, /FILL
>
> ; Draw map grids.
> Map_Grid, MAP_STRUCTURE=map, LONDEL=2.5, LATDEL=0.25, $
> LATLAB=-105, LONLAB=-72, /LABEL
> ;*********************************************************** ****
Did I mention that learning about map projections was,
as far as I can tell, a never-ending proposition?
Here is another little tidbit I learned today. I was a little
uncomfortable with the NSIDC web page information Ken pointed
me to yesterday. It claimed a Polar Stereographic projection
with a True-Scale latitude of -70 degrees. The MAP_PROJ_INIT
function has a TRUE_SCALE_LATITUDE keyword, but it is not
allowed with the Polar Stereographic map projection.
When I checked with the experts at work today, however, I found
out that if you set the CENTER_LATITUDE keyword to -70 for the
Polar Stereographic projection (106), then instead of setting
the center latitude of the projection, it sets the true-scale
latitude. Isn't that neat? And you won't find that in any IDL
documentation, either! :-)
You find it in the bowels of Map_Set_Init if you go though the code
step by step and watch as it changes the parameter list it sends to
the GCTP software.
So now all we have to worry about is whether the reported image
"corners" are located in the center of the pixel, or at the outside
edge of the pixel. No one is saying. I'll have to run some tests
later today to find out (part of the mapx software I mentioned
yesterday). For IDL users, the only way to find out for sure is
good ol' trial and error.
We are almost there. By Friday we should have things lining up
pretty good. :-)
Cheers,
David
--
David Fanning, Ph.D.
Coyote's Guide to IDL Programming (www.dfanning.com)
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
| Re: map_proj_* help [message #66796 is a reply to message #66794] |
Mon, 08 June 2009 05:29 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Ken Mankoff writes:
> Images from the NSIDC. I noticed you used "we" when referring to them.
> Has coyote moved south from Fort Collins and has an affiliation?
Too many children in college. :-(
> Also, I am familiar with a map datum and other terms, although
> obviously not familiar enough to solve this particular problem. But
> are you sure IDL won't convert datums? If you define a "from" and "to"
> on different datums (datii?) and projections in MAP_PROJ_INIT, and
> then in MAP_PROJ_IMAGE use IMAGE_STRUCTURE=from, MAP_STRUCTURE=to,
> this wouldn't work? Never tried it but it seems like it should.
> Anyway, WGS84 is common enough, and my errors are currently large
> enough, that I don't think datum is (yet) the issue.
I don't think IDL does datum shifts. ENVI does a three-point
datum shift, but not the more comprehensive 8-point (?)
datum shifts you can get with more sophisticated map
projection software. (This is still a little blurry for
me, so maybe someone who knows can elaborate.)
> http://nsidc.org/data/iceshelves_images/pine.html
Well, here is how I would orient this image. Trying
to put continental boundaries on Antarctica is sketchy.
But this will certainly give you a map data coordinate
system that you can draw your own map features on.
;**********************************************************
file = 'G:\data\pine_2009084_1545_modis_ch02.png'
image = Read_PNG(file)
image = Reverse(image, 2) ; Put the (0,0) point in UL.
; Display map in window.
pos = [0.1, 0.1, 0.9, 0.9]
TVImage, image, POSITION=pos, /KEEP_ASPECT, /NOINTERP, /ERASE
; Set up Map Projection. Polar Stereographic, WGS-84 Datum.
map = Map_Proj_Init(106, DATUM='WGS 84', CENTER_LATITUDE=-90)
; Convert corners to XY coordinates. Clockwise, from UL.
lons = [-113.594373, -98.868027, -98.568614, -117.004071]
lats = [-71.998860, -72.492200, -76.115095, -75.491196]
xy = Map_Proj_Forward(lons, lats, MAP_STRUCTURE=map)
; Set up map coordinate space for drawing on the map.
Plot, [xy[0,0], xy[0,1]], [xy[1,0], xy[1,2]], $
XStyle=5, YStyle=5, /NoData, POSITION=pos, /NOERASE
; Draw continental outlines.
; Map_Continents, MAP_STRUCTURE=map, /COASTS, /HIRES, /FILL
; Draw map grids.
Map_Grid, MAP_STRUCTURE=map, LONDEL=2.5, LATDEL=0.25, $
LATLAB=-105, LONLAB=-72, /LABEL
;*********************************************************** ****
Of course, in the case of NSIDC image, I *always* believe the
published map corners. ;-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
| Re: map_proj_* help [message #66799 is a reply to message #66796] |
Sun, 07 June 2009 21:13 |
mankoff
Messages: 131
Registered: March 2004
|
Senior Member |
|
|
On Jun 7, 11:51 pm, David Fanning <n...@dfanning.com> wrote:
> mankoff writes:
>> Anyway, thanks for the advice. I'll see if I can get any help from the
>> data publisher, look at gdal, mapx, etc. And I'll be sure to post back
>> here with anything I find out that might be useful to others.
>
> If you want to let me know where I can find these images,
> I'd be happy to have a quick look. One of my hobbies, I guess
> you would call it, is aligning map boundaries on images. ;-)
Images from the NSIDC. I noticed you used "we" when referring to them.
Has coyote moved south from Fort Collins and has an affiliation?
Also, I am familiar with a map datum and other terms, although
obviously not familiar enough to solve this particular problem. But
are you sure IDL won't convert datums? If you define a "from" and "to"
on different datums (datii?) and projections in MAP_PROJ_INIT, and
then in MAP_PROJ_IMAGE use IMAGE_STRUCTURE=from, MAP_STRUCTURE=to,
this wouldn't work? Never tried it but it seems like it should.
Anyway, WGS84 is common enough, and my errors are currently large
enough, that I don't think datum is (yet) the issue.
http://nsidc.org/data/iceshelves_images/pine.html
-k.
|
|
|
|
| Re: map_proj_* help [message #66800 is a reply to message #66799] |
Sun, 07 June 2009 20:51 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
mankoff writes:
> Anyway, thanks for the advice. I'll see if I can get any help from the
> data publisher, look at gdal, mapx, etc. And I'll be sure to post back
> here with anything I find out that might be useful to others.
If you want to let me know where I can find these images,
I'd be happy to have a quick look. One of my hobbies, I guess
you would call it, is aligning map boundaries on images. ;-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
|
|
| Re: map_proj_* help [message #66802 is a reply to message #66801] |
Sun, 07 June 2009 19:45 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Ken Mankoff writes:
> Google Earth 5.x now supports defining the corners of an image
> overlay, not just the edges. Still, it does not line up, presumably
> because this is Polar Stereographic and Google Earth expects
> Equirectangular. My attempts to warp this using map_proj_* got
> nowhere. I've successfully warped a similar image when the 4 middle-of-
> edge points were defined, along with a scale of 250 m/pixel.
I have never been particularly successful warping images from one
projection to another with IDL. At NSIDC we use a program called
mapx to do the warping. It is a UNIX-only piece of software,
hard to install (although freely available), poorly documented,
and it works great after the resident expert shows you how it
works a half-dozen times or so.
But I recommend you go look up Frank Warmerdam and FWTools.
Frank and others have written a suite of tools that do just
about everything it is possible to do with maps and map
projections. GDAL is what you are looking for to do image
warping.
http://fwtools.maptools.org/
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
| Re: map_proj_* help [message #66803 is a reply to message #66802] |
Sun, 07 June 2009 19:38 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Ken Mankoff writes:
> I began converting from using map_set to map_proj_* about a year ago.
> I've used it successfully in a few situations, but still frequently
> find myself wrestling with it and not getting small details to line
> up. Not that small details always line up with MAP_SET...
Well, I would say about 25% of the time, details not lining
up probably has to do with limitations in IDL's map projection
code. But probably 75% of the time, it has to do with not
completely understanding map projections. I will say this,
they are a LOT harder to understand than they appear to be.
I've been doing map projection stuff of well over a year
now, and just when I begin to think I really understand what
I am doing, I find evidence to the contrary. (Maybe you have
read some of my map articles. :-)
One reason your data may not line up is that the lat/lon
values were produced with one datum, and you are doing
the projection with another datum. (If you don't know what
a datum is, the probability of getting things to line up
goes down by another factor of five.) IDL does not really
do datum shifts, and without them points can definitely be
located in different places. (A "latitude" and a "longitude"
are not the same thing as a "meter" or an "acre" and the point
on the Earth they are meant to pinpoint depends on what datum
you are using, and how that datum is oriented with respect to
the center of the Earth.
But, all that said, that is usually not the problem. More often
one of two things has happened. First, and considerably more
likely in my experience, the published "corner" points are
simply wrong. I don't know if it is because of typos or because
undergraduates who know even less about map projections than you
do are doing the calculations, but whatever it is, those numbers
are frequently figments of someone's imagination.
The second thing that happens frequently is that you (perhaps
naively, if you are not that familiar with map projections)
simply assume that the center of the map projection must be
located half way between the maximum and minimum latitude and
longitude. Nothing is ever that simple on a map projection,
believe me. :-)
So, typically, you have to spend several hours (or not infrequently,
several days) tracking down information and confirming every single
bit of it so you know *exactly* what you are doing. Then--and only
if you have been living your life in an exemplary fashion--things
line up like they are suppose to.
Am I discouraging you? I don't mean to. I'm just trying to
point out that you haven't yet tried hard enough to figure
it all out. There are lots of variables. Code, some of it
written by the good folks at ITTVIS, is written incorrectly
(and not fixed in updates of the software, but that is a topic
better left for another time). There are no good maps (a joke)
of the problem terrain.
Perhaps one good place to start is to read all of the map
articles on my web page. There is a LOT of hard earned
experience in those articles, believe me. If you understand
everything you find there, I would say your are probably at
least half way to your destination. Just grit your teeth and
keep going. :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
| Re: map_proj_* help [message #66804 is a reply to message #66803] |
Sun, 07 June 2009 17:55 |
R.G. Stockwell
Messages: 363
Registered: July 1999
|
Senior Member |
|
|
"mankoff" <mankoff@gmail.com> wrote in message
news:639e6c95-085f-4405-b19e-0b8867508c55@e20g2000vbc.google groups.com...
> Dear Mapping Experts,
>
> I began converting from using map_set to map_proj_* about a year ago.
> I've used it successfully in a few situations, but still frequently
> find myself wrestling with it and not getting small details to line
> up. Not that small details always line up with MAP_SET... I'm writing
> today about two specific problems.
>
> ============================================================ ========
>
> 1) I have 3 vectors (lat, lon, data) that I would like to map.
> Specifically, I'd like to export a PNG of the data for Google Earth,
> which means it needs to be in /CYLINDRICAL projection, and as far as I
> know /ISO is not required. For Google Maps I'd need to use /MERCATOR
> but for now I'm just starting with Google Earth.
for google earth, try:
mapStruct = MAP_PROJ_INIT(117, LIMIT=limit, $
CENTER_LONGITUDE=180)
PLOT, mapStruct.uv_box[[0,2]],mapStruct.uv_box[[1,3]] , $
/NODATA, ISOTROPIC=0, XSTYLE=5, YSTYLE=5,$
xtickname = blank_string(10),ytickname=blank_string(10),$
position = [0,0,1,1]
(just as a start, that may help with the "line up" problem.
cheers,
bob
|
|
|
|
| Re: map_proj_* help [message #66862 is a reply to message #66799] |
Fri, 12 June 2009 07:29 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Ken Mankoff writes:
> Images from the NSIDC.
>
> http://nsidc.org/data/iceshelves_images/pine.html
OK, I guess I am giving up on this as a lost cause. I
have tried enough variations now that I am absolutely
convinced it is NOT possible to navigate this image
with the information provided. The numbers don't make
sense.
Here is the information we are given:
Projection: Polar Stereographic
Datum: WGS84
Standard Parallel: -70.0
Corner Coordinates:
UL: 71.998860S 113.594373W
LL: 75.491196S 117.004071W
UR: 72.492200S 98.868027W
LR: 76.115095S 98.568614W
Number of Rows (MODIS): 1600
Number of Columns (MODIS): 2000
Meters per Pixel (MODIS): 250
It seems pretty straight forward to set this map projection up.
Map projection 106 is Polar Stereographic. Map datum 8 is WGS84.
Recall that for this particular map projection, we learned that
setting the CENTER_LONGITUDE actually sets the longitude of true
scale.
map = Map_Proj_Init(106, DATUM=8, CENTER_LAT=-90, CENTER_LON=-70)
Now, let's transform the corner points into XY coordinates. I
order the corner points so that I start in the upper left and
proceed clockwise.
lons = [-113.594373, -98.868027, -98.568614, -117.004071]
lats = [-71.998860, -72.492200, -76.115095, -75.491196]
xy = Map_Proj_Forward(lons, lats, MAP_STRUCTURE=map)
In XY space, each pixel now represents 250 meters. That is to
say, if we start in the UL corner and I multiply 250 times
the number of pixels in the image, I should find an XY number
on the other end of the image. Do we? Let's see.
What do we calculate the range to be in X and Y?
xorigin = xy[0,0]
xend = xy[0,1]
yorigin = xy[1,3]
yend = xy[1,0]
IDL> Print, xorigin, xend
-1397476.0 -951229.17
IDL> Print, yorigin, yend
1110830.3 1467782.8
xrange = Abs(xend - xorigin)
yrange = Abs(yend - yorigin)
If we divide the ranges (in meters nows) by 250 meters/pixel, we
should get the number of pixels in the image (2000x1600). Let's see:
IDL> Print, xrange / 250
1784.987
IDL> Print, yrange / 250
1427.81
Huh!? No comprende!
OK, what if I choose to trust the lat/lon values in the UL corner
of the image, and I'll compute my own XY end points (after all,
I DO know the image really is 2000x1600).
xorigin = xy[0,0]
xend = xorigin + (250*2000)
yend = xy[1,0]
yorigin = yend - (250*1600)
IDL> Print, xorigin, xend
-1397476.0 -897475.98
IDL> Print, yorigin, yend
1067782.8 1467782.8
These results are quite a bit different from what I found before.
Displaying the image with a coordinate system around it and
grids drawn on it, with the XY coordinates I find most believable
(although I can't reproduce the image on the web page with ANY
combination of coordinates), I do this:
pos = [0.1, 0.1, 0.9, 0.9]
TVImage, image, POSITION=pos, /KEEP_ASPECT, /NOINTERP, /ERASE
Plot, [xorigin, xend], [yorigin, yend], POSITION=pos, $
XStyle=5, YStyle=5, /NoData, /NOERASE
Map_Grid, MAP_STRUCTURE=map, LONDEL=5, LATDEL=1, $
LATLAB=-105, LONLAB=-72, /LABEL
Humm. Close, but no cigar.
Conclusion: Bogus information on the web page or ...
I guess the only other conclusion is that I don't really
understand map projections and what I am doing. It used to
be I could readily jump to that conclusion, but anymore, it's
getting to be a hard sell.
I'll see if I can get the real map experts at NSIDC interested
in the problem.
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
