| map_proj_* help [message #66805] |
Sun, 07 June 2009 09:55  |
mankoff
Messages: 131
Registered: March 2004
|
Senior Member |
|
|
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.
I thought I could do this easily with MAP_SET like so:
limit=[min(lat),min(lon),max(lat),max(lon)]
map_set, mean(limit[[0,2]]), mean(limit[[1,3]]), $
limit=limit, /CYL, position=[0,0,1,1], /noborder
plots, lon, lat, data
img =tvrd()
write_png, 'bathy_for_ge.png', img
My logic being that since this is /CYL, ROT=0, and the center of the
map is the center of the limit area, the gridlines should be evenly N/
S and E/W. I can then grab the entire window, save to PNG, and use it
in Google Earth with the same limits. However, my grid lines are not
vertical or horizontal on the screen. Almost, but not quite.
So I think maybe map_proj_* can help, but I'm having trouble getting
them to work. My first attempt so far is this:
proj = map_proj_init( 'Equirectangular', $
limit=limit, $
;; sphere_radius =
;; center_longitude=mean(limit[[1,3]]), $
;;true_scale_latitude, false_easting,
false_northing
;;center_latitude=mean(limit[[0,2]]), $
_EXTRA=e)
result = map_proj_forward( lon, lat, map_structure=proj )
But "result" doesn't make sense to me, I'm not sure what to do with
it.
============================================================ ========
2) I have MODIS images with the following attributes:
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
Number of Rows (AVHRR): 800
Number of Columns (AVHRR): 1000
Meters per Pixel (MODIS): 250
Meters per Pixel (AVHRR): 500
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 can post the code I've used to try to get this working, but it is so
clearly incorrect it might waste your time to read it. I'm sure MODIS
is popular on this group. Any help you can provide will be much
appreciated.
Thanks,
-k.
|
|
|
|
|
|
|
|
|
|
| Re: map_proj_* help [message #66848 is a reply to message #66805] |
Fri, 12 June 2009 12:49  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> I probably should tell you how I solved this, although
> God knows I could probably get rich if I kept it to
> myself. Goes against my nature though, and, anyway,
> I'm storing up karma for the next life. Maybe I can
> be a Matlab expert. ;-)
By the way, if this kind of detective work appeals
to the younger programmers lurking about, I HIGHLY
recommend you go get the book Desert Notes/River Notes
by Barry Lopez and read the story Directions. That
is what it is *always* about with me!
http://tinyurl.com/n23xbb
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 #66849 is a reply to message #66805] |
Fri, 12 June 2009 12:40 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> Oh, never mind. I figured it out. :-(
I probably should tell you how I solved this, although
God knows I could probably get rich if I kept it to
myself. Goes against my nature though, and, anyway,
I'm storing up karma for the next life. Maybe I can
be a Matlab expert. ;-)
What I did was I opened the TIFF file in Frank
Warmerdam's FW_TOOLS application and I looked to see
how he was interpreting the GeoTIFF structure.
If I hadn't insisted on IDL all the damn day, I
probably could have figured it out a bit sooner. :-(
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 #66850 is a reply to message #66805] |
Fri, 12 June 2009 12:32 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> I've come to the tentative conclusion that it is NOT
> possible to navigate this image with IDL map projections.
> I do note that ENVI (which uses its own map projection
> software) *is* able to navigate the image correctly.
>
> This is MOST disconcerting. But not, unfortunately,
> all that surprising to me. I could still be wrong
> (I *hope* I am), but the current seems to be running
> in the other direction. :-(
Oh, never mind. I figured it out. :-(
Here you go. The good news is I probably got another
web page article out of it. :-)
;*****************************************************
file = 'G:\data\pine_2009084_1545_modis_ch02.tif'
image = Read_Tiff(file, Geo=geotag)
image = BytScl(Reverse(image,2))
; 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=8, CENTER_LAT=-70, CENTER_LON=-100)
; Set the limits of the image.
s = Size(image, /Dimensions)
xscale = geotag.ModelPixelScaleTag[0]
yscale = geotag.ModelPixelScaleTag[1]
tp = geotag.ModelTiePointTag[3:4]
xOrigin = tp[0]
yOrigin = tp[1] - (yscale * s[1])
xEnd = xOrigin + (xscale * s[0])
yEnd = tp[1]
; Set up the map coordinate space.
Plot, [xorigin, xend], [yorigin, yend], $
XStyle=5, YStyle=5, /NoData, POSITION=pos, /NOERASE
; Draw map grids.
Map_Grid, MAP_STRUCTURE=map, LONDEL=5, LATDEL=1, $
LATLAB=-110, LONLAB=-75, /LABEL
;*****************************************************
Didn't I tell you we would probably have it whipped by Friday!
This is pretty much how map projections go. At least for me.
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 #66851 is a reply to message #66805] |
Fri, 12 June 2009 11:51 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Ken Mankoff writes:
> Thanks for trying so hard on this one David. I'm sure you'll crack it
> open either as soon as you really give up and not just claim to, or
> they provide the correct numbers.
This is one of my very biggest weaknesses. I need to
learn to move on in the face of insurmountable
technical hurdles. :-(
Partly, I think, I am just reluctant to admit that
IDL--my life for the last 20 years--is becoming
increasingly irrelevant to me. Geez, if I can't
do map projections with it, what's the point?
For what it's worth, I do think the numbers are correct.
I just don't think it is possible to set up this map
projection in IDL. I'm certainly to the point now where
I would be willing to bet some money on it. :-)
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 #66853 is a reply to message #66805] |
Fri, 12 June 2009 11:43 |
Bruce Bowler
Messages: 128
Registered: September 1998
|
Senior Member |
|
|
On Fri, 12 Jun 2009 12:16:50 -0600, David Fanning wrote:
> David Fanning writes:
>
>> I can get closer if I navigate the GeoTiff image you can download from
>> the same place, but I still can't line my grid up with the illustration
>> on the web page:
>
> Actually, this only *looks* better. In fact, it is much, much worse!
>
> I've come to the tentative conclusion that it is NOT possible to
> navigate this image with IDL map projections. I do note that ENVI (which
> uses its own map projection software) *is* able to navigate the image
> correctly.
>
> This is MOST disconcerting. But not, unfortunately, all that surprising
> to me. I could still be wrong (I *hope* I am), but the current seems to
> be running in the other direction. :-(
>
> Cheers,
>
> David
>
> P.S. I'll buy the beer for anyone who can successfully set up a map
> coordinate system in IDL for this image:
>
> ftp://sidads.colorado.edu/pub/DATASETS/ICESHELVES/pine/pine_ 2009084_
> 1545_modis_ch02.tif
Bear in mind that I work primarily with modis ocean data, not land data
so the following may be completely irrelevant.
position to the 6th decimal place is claiming accuracy of less than 4
inches. The smallest pixel size that I'm aware of (this is where the
"oceans" bit comes in) is 250 meters (well, OK, looking at it, it applies
to this image as well).
Given that, I wouldn't, for half a heart beat, believe their coordinates.
Bruce
PS - Let's just say I understand the mapping system even less than David
so I ain't gonna try for the free beer
|
|
|
|
|
|
| Re: map_proj_* help [message #66855 is a reply to message #66805] |
Fri, 12 June 2009 11:16 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> I can get closer if I navigate the GeoTiff image you
> can download from the same place, but I still can't
> line my grid up with the illustration on the web page:
Actually, this only *looks* better. In fact, it is
much, much worse!
I've come to the tentative conclusion that it is NOT
possible to navigate this image with IDL map projections.
I do note that ENVI (which uses its own map projection
software) *is* able to navigate the image correctly.
This is MOST disconcerting. But not, unfortunately,
all that surprising to me. I could still be wrong
(I *hope* I am), but the current seems to be running
in the other direction. :-(
Cheers,
David
P.S. I'll buy the beer for anyone who can successfully
set up a map coordinate system in IDL for this image:
ftp://sidads.colorado.edu/pub/DATASETS/ICESHELVES/pine/pine_ 2009084_
1545_modis_ch02.tif
--
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 #66856 is a reply to message #66805] |
Fri, 12 June 2009 08:47 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> 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.
I can get closer if I navigate the GeoTiff image you
can download from the same place, but I still can't
line my grid up with the illustration on the web page:
;*********************************************************** *****
file = 'G:\data\pine_2009084_1545_modis_ch02.tif'
image = Read_Tiff(file, Geo=geotag)
image = BytScl(Reverse(image,2))
; 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=8, CENTER_LAT=-90, CENTER_LON=-70)
s = Size(image, /Dimensions)
xscale = geotag.ModelPixelScaleTag[0]
yscale = geotag.ModelPixelScaleTag[1]
tp = geotag.ModelTiePointTag[3:4]
xOrigin = tp[0]
yOrigin = tp[1] - (yscale * s[1])
xEnd = xOrigin + (xscale * s[0])
yEnd = tp[1]
Plot, [xorigin, xend], [yorigin, yend], $
XStyle=5, YStyle=5, /NoData, POSITION=pos, /NOERASE
; Draw map grids.
Map_Grid, MAP_STRUCTURE=map, LONDEL=5, LATDEL=1, $
LATLAB=-105, LONLAB=-72, /LABEL
END
;*********************************************************** *****
I do notice in the output from listgeo that there is the following
tag in the file:
ProjStraightVertPoleLongGeoKey (Double 1): -100
I wonder if this has something to do with the problem. If so,
there is nothing I can do about it in IDL. :-(
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 #66962 is a reply to message #66805] |
Mon, 15 June 2009 14:46 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> OK, I've written this up for general consumption.
> As always (probably the main reason I write these
> articles) I learned even more about the problem as
> I wrote it down. Anyway, the article shows you how to
> navigate both the GeoTiff and the PNG version of this
> image. They have to be done in slightly different ways
> due to unspoken "conventions".
>
> You can find the article here:
>
> http://www.dfanning.com/map_tips/iceshelf.html
Just to wrap this whole sorted tale up, I have asked
the good folks at NSIDC to please add the latitude
which is down from the map origin to the information
provided on their map image page. Without this piece
of information, it would be nearly impossible to navigate
the image successfully.
And I wanted to point out that IDL itself is a bit
confusing in this regard, in that the CENTER_LATITUDE
keyword, which normally identifies the latitude of
the map origin (at least it does in every map projection
I have ever worked with) does no such thing with polar
stereographic map projections. Instead, this keyword
identifies the true-scale latitude or standard latitude
(the terms are often used interchangeably). I've added
a bit more detail on this state of affairs in this article:
http://www.dfanning.com/map_tips/polarstereo.html
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 #66967 is a reply to message #66850] |
Sat, 13 June 2009 14:09 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> Oh, never mind. I figured it out. :-(
>
> Here you go. The good news is I probably got another
> web page article out of it. :-)
OK, I've written this up for general consumption. (Left
all the requisite swear words and references to the
Dalai Lama out of it. )
As always (probably the main reason I write these
articles) I learned even more about the problem as
I wrote it down. Anyway, the article shows you how to
navigate both the GeoTiff and the PNG version of this
image. They have to be done in slightly different ways
due to unspoken "conventions".
You can find the article here:
http://www.dfanning.com/map_tips/iceshelf.html
Cheers,
David
P.S. Oh, if you are going to run the code, get a fresh
copy of the Catalyst Library. Things *always* change when
you probe below the surface of things. ;-)
--
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
