| Re: More Map Projection Madness [message #78215] |
Thu, 03 November 2011 11:43  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Chris Torrence writes:
> Please let me know if you see any other issues with the GCTP map projections.
What is wrong with the right side of this Equirectangular
map projection?
IDL> m = map_proj_init(117, datum=19, /gctp)
IDL> map_grid, map_structure=m, linestyle=0
IDL> map_continents, map_structure=m
There seem to be extra lines at odd angles at
the top and bottom of the map on the right-hand
side.
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 #78241 is a reply to message #78215] |
Tue, 01 November 2011 13:20  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> It seems the name "Albers Equal Area" is ambiguous.
> In fact, it seems to choose the old MAP_SET routine
> "Albers Equal Area Conic" as the map projection.
> To make it choose the same projection as projection
> 103, I have to set the CGTP keyword.
>
> alberMap = MAP_PROJ_INIT('albers equal area', /GCTP, $
> DATUM='WGS 84', $
> CENTER_LATITUDE=geotag.PROJNATORIGINLATGEOKEY, $
> CENTER_LONGITUDE=geotag.PROJNATORIGINLONGGEOKEY, $
> STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
> STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY)
>
> Then, this behaves identically to the the projection when
> I use the projection index 103.
This silent substitution of command thing is more
pernicious than I first expected. Not only is the
"albers equal area conic" map projection substituted
for the "albers equal area" projection, but the
sphere radius that is used for the conic map projection
is set to the length of the semi-major axis of the
WGS84 ellipse!
In other words, this command:
mapStruct = Map_Proj_Init('Albers Equal Area', $
ELLIPSOID='WGS 84', $
STANDARD_PAR1=21.0, STANDARD_PAR2=-19.0)
is silently substituted with this command:
mapStruct = Map_Proj_Init('Albers Equal Area Conic', $
SPHERE_RADIUS=6378137.0, $
STANDARD_PAR1=21.0, STANDARD_PAR2=-19.0)
Naturally, you get the wrong results if you don't
realize this. You can find out more about this in
this article:
http://www.idlcoyote.com/map_tips/choosemap.php
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 #78249 is a reply to message #78248] |
Tue, 01 November 2011 10:32 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> But, this contradicts what I learned this morning
> about the WGS84 datum, because here is behaves
> perfectly.
OK. To clarify my clarification. (I guess I'm not
the only person confused by map projections this
morning.)
First, apparently, the UTM projection is the
*only* projection affected by the WGS84 bug,
contrary to what I heard this morning.
Second, there is ambiguity when you set the map
projection to the name "Albers Equal Area". In
fact, Map_Proj_Init will silently select the
"Albers Equal Area Conic" map projection for
you and ignore the value selected with the
DATUM keyword (more below). You get around
this by setting the GCTP keyword to Map_Proj_Init.
The "Albers Equal Area Conic" projection is
the "same as" the "Albers Equal Area" projection,
as far as I can tell from several sources this
morning, but I do note that the old Map_Set
"Albers Equal Area Conic" projection can only
use a spherical datum. Presumably that is the
reason for the different output when these two
different values are used to select the map
projection.
I also note that there would be less confusion
if this selection and the ignoring of the DATUM keyword
were flagged for the user. You don't typically
see the "Albers Equal Area Conic" selection if you
are using on-line help to look for projections to use
with Map_Proj_Init. And you certainly wouldn't
be aware that they produce different results.
In fact, you probably will only learn this when your
paper with the incorrect results has gone to print. :-(
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 #78251 is a reply to message #78248] |
Tue, 01 November 2011 09:35 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> 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. :-(
OK, the situation is more complicated than this.
Apparently, I chose a perverse example. Although,
I have to admit, over the past couple of months
I have the distinct feeling that examples are
choosing me, rather than visa versa. :-(
It seems the name "Albers Equal Area" is ambiguous.
In fact, it seems to choose the old MAP_SET routine
"Albers Equal Area Conic" as the map projection.
To make it choose the same projection as projection
103, I have to set the CGTP keyword.
alberMap = MAP_PROJ_INIT('albers equal area', /GCTP, $
DATUM='WGS 84', $
CENTER_LATITUDE=geotag.PROJNATORIGINLATGEOKEY, $
CENTER_LONGITUDE=geotag.PROJNATORIGINLONGGEOKEY, $
STANDARD_PAR1=geotag.PROJSTDPARALLEL1GEOKEY, $
STANDARD_PAR2=geotag.PROJSTDPARALLEL2GEOKEY)
Then, this behaves identically to the the projection when
I use the projection index 103.
But, this contradicts what I learned this morning
about the WGS84 datum, because here is behaves
perfectly.
I appreciate Chris's attempt to shed light on this subject.
It is greatly appreciated. But, I am still very confused.
I understand this may not be completely ITTVIS's fault. :-)
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 #78283 is a reply to message #78215] |
Mon, 07 November 2011 08:19 |
Paul Van Delst[1]
Messages: 1157
Registered: April 2002
|
Senior Member |
|
|
David Fanning wrote:
> Chris Torrence writes:
>
>> Please let me know if you see any other issues with the GCTP map projections.
>
> What is wrong with the right side of this Equirectangular
> map projection?
>
> IDL> m = map_proj_init(117, datum=19, /gctp)
> IDL> map_grid, map_structure=m, linestyle=0
> IDL> map_continents, map_structure=m
>
> There seem to be extra lines at odd angles at
> the top and bottom of the map on the right-hand
> side.
Huh. Apart from an unlabeled X-Y axis from map_grid, I don't get ANY output from the above.
Weird.
cheers,
paulv
|
|
|
|
| Re: More Map Projection Madness [message #78302 is a reply to message #78215] |
Fri, 04 November 2011 10:34 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> What is wrong with the right side of this Equirectangular
> map projection?
>
> IDL> m = map_proj_init(117, datum=19, /gctp)
> IDL> map_grid, map_structure=m, linestyle=0
> IDL> map_continents, map_structure=m
>
> There seem to be extra lines at odd angles at
> the top and bottom of the map on the right-hand
> side.
As I have received no answer, I thought perhaps my
example wasn't working. Here is code that reliable
shows the problem I am having with map projections.
In particular, when I set up the projection with
Map_Proj_Init and try to draw grid lines, I am
getting extraneous and ill-positioned lines showing
up on the plot. If you look at the right-hand side
of the third window, you can see an example of
what I mean.
;**********************************************************
lons = -180 > Indgen(11)*36-180 < 180
lats = -90 > Indgen(11)*18-90 < 90
window, 0, Title='Normal Map Grid'
map_set, /Cylindrical, Position=[0.1, 0.1, 0.9, 0.9]
map_grid, lats=lats, lons=lons
window, 1, Title='OK GCTP Map Grid'
m = Map_Proj_Init('Miller Cylindrical', Datum='Sphere', /GCTP)
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
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
END
;**********************************************************
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
