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
