| map_proj_init [message #53364] |
Sat, 07 April 2007 15:09  |
marit
Messages: 4
Registered: September 1995
|
Junior Member |
|
|
Has anyone run into the problem that subsequent calls to map_proj_init
using a GCTP projection interfere with previously defined maps? Here
is an example:
south_proj=MAP_PROJ_INIT('Polar Stereographic' ,/GCTP ,$
semimajor_axis=6378273.0,semiminor_axis=6356889.4,$
center_lon=0,center_lat=-70.0,false_easting=0,false_northing =0)
print,map_proj_forward([0,0],[-90.0,-89.0],map_structure=sou th_proj)
0.0000000 0.0000000
0.0000000 108332.24
north_proj=MAP_PROJ_INIT('Polar Stereographic' ,/GCTP ,$
semimajor_axis=6378273.0,semiminor_axis=6356889.4,$
center_lon=0,center_lat=70.0,false_easting=0,false_northing= 0)
print,map_proj_forward([0,0],[-90.0,-89.0],map_structure=sou th_proj)
0.0000000 -2.0002841e+23
0.0000000 -1.4035070e+09
I haven't yet seen this occur if the second projection is an IDL
projection; however the IDL projections are not usually useful since
they can't be set up like a normal projections with false easting and
false northing and they are mostly spherical not ellipsoidal.
|
|
|
|
| Re: MAP_PROJ_INIT [message #78186 is a reply to message #53364] |
Sat, 29 October 2011 11:41  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
alx writes:
> In order to process some GPS data, I logically write:
> IDL> map = map_proj_init(101, ELLIPSOID=24, ZONE=31)
> wanting to use WGS84 ellipsoid and UTM projection (zone 31 is for
> Paris), for further use of the "map_proj_forward" function. Then :
> IDL> print, map.A, map.E2, map.PROJECTION
> 6370997.0 0.00000000 20
> shows that IDL rather chooses the SPHERE and the projection n°20.
> Moreover this projection is not referenced in the IDL_help, the
> projection index ranging from 0 to 19.
> Forcing GCTP keyword to 1 does not change anything.
> What does it mean ?
Well, I think it means the IDL UTM map transformation is
screwed. You can ONLY use a spherical datum with this
transformation. The datum I was trying to use the other
day was WGS84 and I was getting the magnitude of errors
you are getting. This is a problem, since ANY data you
get from satellites is probably using a WGS84 datum.
Which means there is no way you can use this data in
IDL. (Well, unless accuracy doesn't matter to you.)
Maybe this is why IDL examples always use global map
projections. A couple hundred meters here and there
is chump change on that scale.
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: MAP_PROJ_INIT [message #78187 is a reply to message #53364] |
Sat, 29 October 2011 11:10 |
lecacheux.alain
Messages: 325
Registered: January 2008
|
Senior Member |
|
|
On 29 oct, 19:14, alx <lecacheux.al...@wanadoo.fr> wrote:
> On 29 oct, 18:48, David Fanning <n...@dfanning.com> wrote:
>
>
>
>
>
>> alx writes:
>>> My problem is actually the following: when applying the
>>> map_proj_forward function to real data,- that is to say when computing
>>> plane coordinates from given GPS longitudes and latitudes, a very
>>> standard operation from GPS data -, I find results wrong by a few
>>> hundred meters...
>>> For instance:
>>> IDL> lon0 = 2.1937863d0 & lat0 = 47.3808737d0
>>> IDL> map = map_proj_init(101, ELLIPSOID=24, ZONE=31)
>>> IDL> print, map_proj_forward(lon0, lat0)
>>> 439142.34 5247587.0
>>> Actual easting and northing (as calculated by using other tools) are
>>> 439144 and 5247806 instead.
>>> I then suspect somme error in map_proj_init, or maybe I am not using
>>> it correctly.
>>> I was wondering if other people had same experience.
>
>> Well, this is *exactly* the problem I was having this
>> week, and the "other tool" was ENVI. Maybe we are going
>> to have to look into this some more. What "other tools"
>> are you using?
>
>> 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.")- Masquer le texte des messages précédents -
>
>> - Afficher le texte des messages précédents -
>> What "other tools" are you using?
>
> Several online tools, like "http://www.rcn.montana.edu/resources/tools/
> coordinates.aspx".
> But my reference is the french one: "http://geodesie.ign.fr/index.php?
> p=53&page=circe".
> All these give same results within 1 meter.
> alx.- Masquer le texte des messages précédents -
>
> - Afficher le texte des messages précédents -
Well, I ran the command "map = map_proj_init(101, ELLIPSOID=24,
ZONE=31)" in debug mode.
I found that the ELLIPSOID index is wrongly changed to -1 at line
1164 (in "map_proj_getellipsoid"). This may explain why the output map
finally contains a SPHERE reference.
Unavoidably, the rest of the computation will be in error or, at
least, inaccurate.
The ITTVIS programmer seems to feel a bit unconfortable, because on
lines 1376-1379 (in map_proj_init) you can get his own warning :
; Internal routines to initialize GCTP forward and inverse
; projections. Use at your own peril.
MAP_PROJ_GCTP_FORINIT, sMap.simple, gctpZone, sMap.p,
ellipsoid
MAP_PROJ_GCTP_REVINIT, sMap.simple, gctpZone, sMap.p,
ellipsoid
Indeed, sMap.p does no longer contains the original ellipsoid
index ...
alx.
|
|
|
|
| Re: MAP_PROJ_INIT [message #78189 is a reply to message #53364] |
Sat, 29 October 2011 10:14 |
lecacheux.alain
Messages: 325
Registered: January 2008
|
Senior Member |
|
|
On 29 oct, 18:48, David Fanning <n...@dfanning.com> wrote:
> alx writes:
>> My problem is actually the following: when applying the
>> map_proj_forward function to real data,- that is to say when computing
>> plane coordinates from given GPS longitudes and latitudes, a very
>> standard operation from GPS data -, I find results wrong by a few
>> hundred meters...
>> For instance:
>> IDL> lon0 = 2.1937863d0 & lat0 = 47.3808737d0
>> IDL> map = map_proj_init(101, ELLIPSOID=24, ZONE=31)
>> IDL> print, map_proj_forward(lon0, lat0)
>> 439142.34 5247587.0
>> Actual easting and northing (as calculated by using other tools) are
>> 439144 and 5247806 instead.
>> I then suspect somme error in map_proj_init, or maybe I am not using
>> it correctly.
>> I was wondering if other people had same experience.
>
> Well, this is *exactly* the problem I was having this
> week, and the "other tool" was ENVI. Maybe we are going
> to have to look into this some more. What "other tools"
> are you using?
>
> 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.")- Masquer le texte des messages précédents -
>
> - Afficher le texte des messages précédents -
> What "other tools" are you using?
Several online tools, like "http://www.rcn.montana.edu/resources/tools/
coordinates.aspx".
But my reference is the french one: "http://geodesie.ign.fr/index.php?
p=53&page=circe".
All these give same results within 1 meter.
alx.
|
|
|
|
| Re: MAP_PROJ_INIT [message #78190 is a reply to message #53364] |
Sat, 29 October 2011 09:48 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
alx writes:
> My problem is actually the following: when applying the
> map_proj_forward function to real data,- that is to say when computing
> plane coordinates from given GPS longitudes and latitudes, a very
> standard operation from GPS data -, I find results wrong by a few
> hundred meters...
> For instance:
> IDL> lon0 = 2.1937863d0 & lat0 = 47.3808737d0
> IDL> map = map_proj_init(101, ELLIPSOID=24, ZONE=31)
> IDL> print, map_proj_forward(lon0, lat0)
> 439142.34 5247587.0
> Actual easting and northing (as calculated by using other tools) are
> 439144 and 5247806 instead.
> I then suspect somme error in map_proj_init, or maybe I am not using
> it correctly.
> I was wondering if other people had same experience.
Well, this is *exactly* the problem I was having this
week, and the "other tool" was ENVI. Maybe we are going
to have to look into this some more. What "other tools"
are you using?
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: MAP_PROJ_INIT [message #78191 is a reply to message #53364] |
Sat, 29 October 2011 08:27 |
lecacheux.alain
Messages: 325
Registered: January 2008
|
Senior Member |
|
|
On 29 oct, 16:18, David Fanning <n...@dfanning.com> wrote:
> alx writes:
>> In order to process some GPS data, I logically write:
>> IDL> map = map_proj_init(101, ELLIPSOID=24, ZONE=31)
>> wanting to use WGS84 ellipsoid and UTM projection (zone 31 is for
>> Paris), for further use of the "map_proj_forward" function. Then :
>> IDL> print, map.A, map.E2, map.PROJECTION
>> 6370997.0 0.00000000 20
>> shows that IDL rather chooses the SPHERE and the projection n°20.
>> Moreover this projection is not referenced in the IDL_help, the
>> projection index ranging from 0 to 19.
>> Forcing GCTP keyword to 1 does not change anything.
>> What does it mean ?
>
> I doubt it means anything. At most it means IDL maintains
> a different indexing scheme internally than they do in
> their public interface. That's not unusual. I don't think
> I would spend any time worrying about it. :-)
>
> 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.")
My problem is actually the following: when applying the
map_proj_forward function to real data,- that is to say when computing
plane coordinates from given GPS longitudes and latitudes, a very
standard operation from GPS data -, I find results wrong by a few
hundred meters...
For instance:
IDL> lon0 = 2.1937863d0 & lat0 = 47.3808737d0
IDL> map = map_proj_init(101, ELLIPSOID=24, ZONE=31)
IDL> print, map_proj_forward(lon0, lat0)
439142.34 5247587.0
Actual easting and northing (as calculated by using other tools) are
439144 and 5247806 instead.
I then suspect somme error in map_proj_init, or maybe I am not using
it correctly.
I was wondering if other people had same experience.
alain.
|
|
|
|
| Re: MAP_PROJ_INIT [message #78192 is a reply to message #53364] |
Sat, 29 October 2011 07:18 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
alx writes:
> In order to process some GPS data, I logically write:
> IDL> map = map_proj_init(101, ELLIPSOID=24, ZONE=31)
> wanting to use WGS84 ellipsoid and UTM projection (zone 31 is for
> Paris), for further use of the "map_proj_forward" function. Then :
> IDL> print, map.A, map.E2, map.PROJECTION
> 6370997.0 0.00000000 20
> shows that IDL rather chooses the SPHERE and the projection n°20.
> Moreover this projection is not referenced in the IDL_help, the
> projection index ranging from 0 to 19.
> Forcing GCTP keyword to 1 does not change anything.
> What does it mean ?
I doubt it means anything. At most it means IDL maintains
a different indexing scheme internally than they do in
their public interface. That's not unusual. I don't think
I would spend any time worrying about it. :-)
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
