| 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.")
|
|
|
|
| Re: More Map Projection Madness [message #78281 is a reply to message #78253] |
Mon, 07 November 2011 09:20  |
Paul Van Delst[1]
Messages: 1157
Registered: April 2002
|
Senior Member |
|
|
Hello,
David Fanning wrote:
> Paul van Delst writes:
>
>> Huh. Apart from an unlabeled X-Y axis from map_grid, I don't get ANY output from the above.
>
> Say what!?
>
>> Weird.
>
> Indeed.
>
> Just to let you know, I've spent the past week working
> on Coyote Graphics map projection routines. I've fixed
> all of the problems I know about in MAP_GRID (and there
> were many!) and I have turned these into routines that
> work in the way of other Coyote Graphics routines.
>
> I'll be making an announcement shortly. I'm working
> on polishing these routines up and writing the
> documentation now.
>
> Of course, if MAP_GRID and MAP_CONTINENTS don't work
> on UNIX machines anymore, I'm screwed, so please fix
> this problem ASAP. :-)
I already debug several other licensed, and not inexpensive, software products thanks very much. :o)
If it makes a difference, your other example with the three windows produces output consistent with your text explaining
the problem.
cheers,
paulv
|
|
|
|
| Re: More Map Projection Madness [message #78282 is a reply to message #78253] |
Mon, 07 November 2011 08:26 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Paul van Delst writes:
> Huh. Apart from an unlabeled X-Y axis from map_grid, I don't get ANY output from the above.
Say what!?
> Weird.
Indeed.
Just to let you know, I've spent the past week working
on Coyote Graphics map projection routines. I've fixed
all of the problems I know about in MAP_GRID (and there
were many!) and I have turned these into routines that
work in the way of other Coyote Graphics routines.
I'll be making an announcement shortly. I'm working
on polishing these routines up and writing the
documentation now.
Of course, if MAP_GRID and MAP_CONTINENTS don't work
on UNIX machines anymore, I'm screwed, so please fix
this problem ASAP. :-)
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.")
|
|
|
|
| Re: More Map Projection Madness [message #78294 is a reply to message #78253] |
Fri, 04 November 2011 14:23 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> window, 2, Title='Screwed Up GCTP Map Grid'
> m = Map_Proj_Init('Cylindrical')
> plot, [0,1], xrange=m.uv_box[[0,2]], yrange=m.uv_box[[1,3]], $
> xstyle=5, ystyle=5, /noerase, /nodata, $
> Position=[0.1, 0.1, 0.9, 0.9]
> map_grid, lats=lats, lons=lons, map_structure=m
The problem here appears to be a bug I first reported
to ITTVIS back in December of 2009.
Map_Proj_Forward has a POLYLINES keyword that returns correct
values itself, but causes the projected XY latitude values to
be incorrect. Thus, they get drawn in the wrong locations.
The only work-around is to call Map_Proj_Forward *without*
using this keyword. For example, in Map_Grid, you can call
Map_Proj_Forward twice, like this:
uv = MAP_PROJ_FORWARD(REPLICATE(lon, N_ELEMENTS(lati)), lati, $
MAP_STRUCTURE=mapStruct, POLYLINES=polylines)
uv = MAP_PROJ_FORWARD(REPLICATE(lon, N_ELEMENTS(lati)), lati, $
MAP_STRUCTURE=mapStruct)
And this problem goes away.
Here is an example of what I mean:
IDL> lons = [-180, -90, 0, 90, 180]
IDL> lats = Replicate(90, 5)
IDL> m = Map_Proj_Init('Cylindrical')
IDL> uv = Map_Proj_Forward(lons, lats, map_structure=m, polylines=p)
IDL> print, uv[1,*]
10007523.
10007523.
10007539.
10007539.
10007539.
10007539.
10007523.
IDL> uv = Map_Proj_Forward(lons, lats, map_structure=m)
IDL> print, uv[1,*]
10007539.
10007539.
10007539.
10007539.
10007539.
You can see there are incorrect values (10007523.) in the
first Map_Proj_Forward call that are not in the second call.
Plus there are fewer values!!
Very strange!
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
