| More Map Projection Madness [message #78253] |
Tue, 01 November 2011 08:58  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Folks,
I have really bad news today.
I am still trying to get to the bottom of problems
I am having with map projections producing incorrect
results. Today I turned my attention to an image
that uses an Albers Equal Area projection, rather than
a UTM projection.
I used a data file, AF03sep15b.n16-VIg.tif, that you
can download from here, if you want to check my work:
http://www.idlcoyote.com/data/
Here is the code I ran:
geofile = 'AF03sep15b.n16-VIg.tif'
geoImage= Read_Tiff(geoFile, GeoTIFF=geotag)
geoImage= Reverse(geoImage, 2)
xscale = geotag.ModelPixelScaleTag[0]
yscale = geotag.ModelPixelScaleTag[1]
tp = geotag.ModelTiePointTag[3:4]
s = Size(geoImage, /Dimensions)
xrange = [tp[0], tp[0] + (xscale * s[0])]
yrange = [tp[1] - (yscale * s[1]), tp[1]]
alberMap = MAP_PROJ_INIT('albers equal area', $
DATUM='WGS 84', $
CENTER_LATITUDE=geotag.PROJNATORIGINLATGEOKEY, $
CENTER_LONGITUDE=geotag.PROJNATORIGINLONGGEOKEY, $
STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY)
Print, map_proj_inverse(tp[0], tp[1], map_structure=alberMap)
The results I get are these:
-24.538705 43.358419
I know these values to be wrong. The correct values are:
-24.521589 43.618949
I checked with ITTVIS technical support to see if I had
misunderstood what I was told yesterday about the UTM
projection and the WGS84 datum. I did. They claim that
the WGS84 datum is not working for *any* map projection.
So, okay. I substituted the WALBECK datum for the WGS84
datum and ran the program again. I got the same
incorrect results! That's weird!
So, I was preparing a note to send to technical support,
when I noticed that in my example to them, I was suddenly
getting the correct result. What was different in the code
I was using now?
Well, instead of using the name of the map projection,
"Albers Equal Area", as I usually do for readability and
pedological reasons, I was taking a short-cut and substituting
the map projection index number, 103. That was the
*only* difference!
That's right, this code will produce accurate results:
alberMap = MAP_PROJ_INIT(103, $
DATUM='WGS 84', $
CENTER_LATITUDE=geotag.PROJNATORIGINLATGEOKEY, $
CENTER_LONGITUDE=geotag.PROJNATORIGINLONGGEOKEY, $
STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY)
This goes a long way in explaining why I have seriously
thought I was going crazy in the past week or so! I keep
getting different results from what I think is exactly
the same code!
Unfortunately, this does NOT apply to the UTM projection
problem from yesterday, where the WGS84 datum is, in
fact, totally screwed.
In fact, I'm really afraid map projections in general
in IDL are totally screwed. How can we ever be sure
we are getting correct results!?
So, bottom line. Don't use the UTM projection with
the WGS84 datum, don't use map projection names at
all, and get VERY familiar with proj4 so you can
check to see if anything at all that comes out of
an IDL map projection is accurate. :-(
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
