[GRASS-user] gps import: definition of WGS84

GRASS GRASS user
1

Hello,

I'm trying to figure out why we are seeing different results when importing GPS data into different software (GRASS is doing great, it's the others that are having the trouble). It must be linked to the definition of the projection the GPS data comes in.

v.in.gpsbabel defines the default projection as

IN_PROJ="+proj=longlat +towgs84=0.000,0.000,0.000"

Now my question: is this equivalent to EPSG 4326, which is defined as

<4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs <>

?

And if not, which EPSG code would be the equivalent of the v.in.gpsbabel line ?

Moritz

2

On Fri, 2008-05-23 at 18:39 +0200, Moritz Lennert wrote:

Hello,

I'm trying to figure out why we are seeing different results when
importing GPS data into different software (GRASS is doing great, it's
the others that are having the trouble). It must be linked to the
definition of the projection the GPS data comes in.

v.in.gpsbabel defines the default projection as

IN_PROJ="+proj=longlat +towgs84=0.000,0.000,0.000"

Now my question: is this equivalent to EPSG 4326, which is defined as

<4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs <>

?

Dear Moritz,

I think the first definition is not complete. If the "+datum=WGS84"
parameter was at least present it should be ok since a datum includes
(or should include) all required information.

And if not, which EPSG code would be the equivalent of the v.in.gpsbabel
line ?

Not sure about that. I don't think that an incomplete definition has an
EPSG code.

There are of course EPSG codes for
.the Spheroid

.transformation parameters from ??? to wgs84

.the prime meridian

.the unit of measure

.specific datums (which should include all above normally).

Maybe I am wrong,
please correct me.

Cheers,

Nikos

Moritz
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

3

On 23/05/08 18:45, Nikos Alexandris wrote:

On Fri, 2008-05-23 at 18:39 +0200, Moritz Lennert wrote:

Hello,

I'm trying to figure out why we are seeing different results when importing GPS data into different software (GRASS is doing great, it's the others that are having the trouble). It must be linked to the definition of the projection the GPS data comes in.

v.in.gpsbabel defines the default projection as

IN_PROJ="+proj=longlat +towgs84=0.000,0.000,0.000"

Now my question: is this equivalent to EPSG 4326, which is defined as

<4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs <>

?

Dear Moritz,

I think the first definition is not complete. If the "+datum=WGS84"
parameter was at least present it should be ok since a datum includes
(or should include) all required information.

The "problem" is that it's the first that works.

To be more precise:

I have data in a GPS.

I have a raster topographic map of the area where the GPS data was taken. This map is in ESPG 31370.

I create my EPSG 31370 location and use, import my raster and use v.in.gpsbabel to import the gps data, not setting the proj= parameter.

The gps points and tracks are perfectly positioned compared to the topo map.

However, when I try to use QGIS (or gvSIG) with on-the-fly projection, defining the projection of my tracks/points as EPSG 4326 (and I've tried a few others), or as a custom projection defined as

+proj=longlat +towgs84=0.000,0.000,0.000 +ellps=WGS84

my topo map projection as 31370, and my project/view projection as 31370, I always see my GPS data several tens of meters off...

Moritz

4

On Fri, 2008-05-23 at 19:12 +0200, Moritz Lennert wrote:

However, when I try to use QGIS (or gvSIG) with on-the-fly
projection,
defining the projection of my tracks/points as EPSG 4326 (and I've
tried
a few others), or as a custom projection defined as

+proj=longlat +towgs84=0.000,0.000,0.000 +ellps=WGS84

my topo map projection as 31370, and my project/view projection as
31370, I always see my GPS data several tens of meters off...

Moritz,

could it be a QGIS problem?

I assume you have set up the projection in GGIS to 31370.

What happens if you open your GRASS mapset in QGIS and load the GPS data
as a GRASS vector (with on-the-fly)?

5

I'm trying to figure out why we are seeing different results when
importing GPS data into different software (GRASS is doing
great, it's the others that are having the trouble). It must be
linked to the definition of the projection the GPS data comes in.

v.in.gpsbabel defines the default projection as

IN_PROJ="+proj=longlat +towgs84=0.000,0.000,0.000"

Now my question: is this equivalent to EPSG 4326, which is
defined as

<4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs <>

?

It is intended to be the same as GPSs use WGS84. With the changes to the latest version of PROJ.4 I am not 100% sure that the semantics are not subtly changed for that. (issue of datum transform only happens if datum is/is not present in both "to" and "from" projection definitions; see the proj release notes + news for more on that)
I had /thought/ that the +towgs84= was enough to trigger a datum shift.

And if not, which EPSG code would be the equivalent of the
v.in.gpsbabel line ?

you can try modifying the v.in.gpsbabel script to use WGS84's epsg code if you like and see if you get a difference. what version of PROJ.4 do you use?

another test is to run v.in.gpsbabel into a wgs84 location, then use v.proj to project to your target projection, and see if it lines up exactly with the direct v.in.gpsbabel (with reprojection) way.

If you are only seeing error of a few meters, it may be that the maths of the reverse-projection used by the on-the-fly software is not the same in both directions. e.g. see the forward and reverse transform error given by the gis.m georeferencing tool.

If it is an error of approx 100m, I would expect that to be an incorrect datum issue.

I would note that with current ArcGIS 9 on-the-fly reprojection that* datum shift is totally ignored and only reprojection is done, leading to badly incorrect results. Couple that with "user friendly" rubber sheet drag-and-drop georectifying to match the incorrect reprojection and you get a few wasted afternoons trying to figure out what the hell is going on.

[*] at least for our tests with the NZGD49 datum, no idea if it works with better tested datums like NAD27

Hamish

6

Hamish wrote:

I'm trying to figure out why we are seeing different results when importing GPS data into different software (GRASS is doing
great, it's the others that are having the trouble). It must be
linked to the definition of the projection the GPS data comes in.

v.in.gpsbabel defines the default projection as

IN_PROJ="+proj=longlat +towgs84=0.000,0.000,0.000"

Now my question: is this equivalent to EPSG 4326, which is
defined as

<4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs <>

?

It is intended to be the same as GPSs use WGS84. With the changes to the latest version of PROJ.4 I am not 100% sure that the semantics are not subtly changed for that. (issue of datum transform only happens if datum is/is not present in both "to" and "from" projection definitions; see the proj release notes + news for more on that)
I had /thought/ that the +towgs84= was enough to trigger a datum shift.

Folks,

I'm not sure what the implications are in this context, but any valid
PROJ.4 coordinate system ought to include the ellipsoid definition - either
directly or via the datum parameter. So if the coordinate system is
intended to be WGS84 I'd encourage it to be declared as
"+proj=latlong +datum=WGS84" rather than "+proj=latlong +towgs84=0,0,0".

In normal practice a default ellipsoid of WGS84 is used, and in combination
with +towgs84=0,0,0 this would be equivelent to +datum=WGS84 - but why
risk the defaults file not being found?

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | President OSGeo, http://osgeo.org

7

Frank:

I'm not sure what the implications are in this context,
but any valid PROJ.4 coordinate system ought to include
the ellipsoid definition - either directly or via the
datum parameter. So if the coordinate system is
intended to be WGS84 I'd encourage it to be declared as
"+proj=latlong +datum=WGS84" rather than
"+proj=latlong +towgs84=0,0,0".

In normal practice a default ellipsoid of WGS84 is used,
and in combination with +towgs84=0,0,0 this would be
equivelent to +datum=WGS84 - but why risk the defaults
file not being found?

ok, proj strings updated in SVN for scripts that call cs2cs.

(I take it there is no difference between latlong and longlat? ie it always assumes input is x,y unless you use proj/cs2cs -r ?)

$ proj -l
....
lonlat : Lat/long (Geodetic)
latlon : Lat/long (Geodetic alias)

Hamish

8

I am experienced some trouble with v.distance in my win box. I am running the latest version 6.3.0
my tables are in PG, my db.connect gives:

db.connect -p
driver:pg
database:host=localhost,dbname=thesisWA,port=5432
schema:public
group:

but when I run v.distance I get:
$ v.distance from=cent_g to=road_inv_d6_b dmax=5 upload=cat,dist column=cat_ce
nt,distance
DBMI-Postgres driver error:
Cannot connect to Postgres: could not connect to server: Connection refused (0x0000274D/10061)
        Is the server running on host "localhost" and accepting
        TCP/IP connections on port 6543?

GRASS_INFO_WARNING(4116,1): Cannot open database 'host=localhost,dbname=thesisWA,port=6543'
GRASS_INFO_END(4116,1)

GRASS_INFO_ERROR(4116,2): Unable to open database <host=localhost,dbname=thesisWA,port=6543> by driver <pg>
GRASS_INFO_END(4116,2)

It looks like somehow the port is hardcoded to 6543 in v.distance. If I tried db.tables -p it gives me the list of my tables so it looks like my connection settings are fine.

Jonathan Aguero-Valverde
PhD. Candidate
The Thomas D. Larson
Pennsylvania Transportation Institute
The Pennsylvania State University

9

Thanks to everyone for their answers. I've tried all the different
suggestions and the answer seems to be clear: on-the-fly projection does
not work properly in QGIS and gvSIG... Whatever I do in GRASS, I always
get correct placement of the GPS data.

Using on-the-fly projection in the other programs, I get a difference of
about 100m, so I guess it's a datum issue as Hamish suggests.

Moritz

On Sat, May 24, 2008 05:05, Hamish wrote:

I'm trying to figure out why we are seeing different results when
importing GPS data into different software (GRASS is doing
great, it's the others that are having the trouble). It must be
linked to the definition of the projection the GPS data comes in.

v.in.gpsbabel defines the default projection as

IN_PROJ="+proj=longlat +towgs84=0.000,0.000,0.000"

Now my question: is this equivalent to EPSG 4326, which is
defined as

<4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs <>

?

It is intended to be the same as GPSs use WGS84. With the changes to the
latest version of PROJ.4 I am not 100% sure that the semantics are not
subtly changed for that. (issue of datum transform only happens if datum
is/is not present in both "to" and "from" projection definitions; see the
proj release notes + news for more on that)
I had /thought/ that the +towgs84= was enough to trigger a datum shift.

And if not, which EPSG code would be the equivalent of the
v.in.gpsbabel line ?

you can try modifying the v.in.gpsbabel script to use WGS84's epsg code if
you like and see if you get a difference. what version of PROJ.4 do you
use?

another test is to run v.in.gpsbabel into a wgs84 location, then use
v.proj to project to your target projection, and see if it lines up
exactly with the direct v.in.gpsbabel (with reprojection) way.

If you are only seeing error of a few meters, it may be that the maths of
the reverse-projection used by the on-the-fly software is not the same in
both directions. e.g. see the forward and reverse transform error given by
the gis.m georeferencing tool.

If it is an error of approx 100m, I would expect that to be an incorrect
datum issue.

I would note that with current ArcGIS 9 on-the-fly reprojection that*
datum shift is totally ignored and only reprojection is done, leading to
badly incorrect results. Couple that with "user friendly" rubber sheet
drag-and-drop georectifying to match the incorrect reprojection and you
get a few wasted afternoons trying to figure out what the hell is going
on.

[*] at least for our tests with the NZGD49 datum, no idea if it works with
better tested datums like NAD27

Hamish

10

On Sun, 2008-05-25 at 12:43 +0200, Moritz Lennert wrote:

Thanks to everyone for their answers. I've tried all the different
suggestions and the answer seems to be clear: on-the-fly projection does
not work properly in QGIS and gvSIG... Whatever I do in GRASS, I always
get correct placement of the GPS data.

Sort of off-topic:

I like that GRASS is strict with the projection of data. You are 100%
sure that all data are there where you see them. I don't think it's
really necessary to add an on-the-fly function in GRASS. You can always
re-project and know what happens.

Nikos

11

On Sun, May 25, 2008 12:47, Nikos Alexandris wrote:

On Sun, 2008-05-25 at 12:43 +0200, Moritz Lennert wrote:

Thanks to everyone for their answers. I've tried all the different
suggestions and the answer seems to be clear: on-the-fly projection does
not work properly in QGIS and gvSIG... Whatever I do in GRASS, I always
get correct placement of the GPS data.

Sort of off-topic:

I like that GRASS is strict with the projection of data. You are 100%
sure that all data are there where you see them. I don't think it's
really necessary to add an on-the-fly function in GRASS. You can always
re-project and know what happens.

I agree 100%. I'm not pleading for on-the-fly projection at all. However,
since I'm currently teaching a training course which also includes QGIS
and gvSIG, I also use that feature, but apparently I'd better not !

Moritz

12

Moritz Lennert

Thanks to everyone for their answers. I've tried all the
different suggestions and the answer seems to be clear:
on-the-fly projection does not work properly in QGIS and gvSIG...
Whatever I do in GRASS, I always get correct placement of the GPS data.

This seems to be QGIS bug # 1079
  http://trac.osgeo.org/qgis/ticket/1079

(usual suspects reporting in :slight_smile:

Using on-the-fly projection in the other programs, I get a
difference of about 100m, so I guess it's a datum issue as
Hamish suggests.

It would seem so. I can reproduce it here (see bug report). In my case looking in the layer's "P"roperies projection tab I see the the +datum= and/or +towgs84= terms are missing* versus what is found in the EPSG file.

[*] file: /usr/share/qgis/resources/srs.db (SQLite DB)

Root cause seems to be that QGIS lacks a "pick appropriate datum transform parameters to use" GUI in the case when there are multiple possibilities, and so it doesn't try any datum transform at all:
http://www.nabble.com/why-is-EPSG%3A2180-not-recognised--td16505409.html#a16592790

Hamish

13

On Sun, May 25, 2008 15:16, Hamish wrote:

Moritz Lennert

Thanks to everyone for their answers. I've tried all the
different suggestions and the answer seems to be clear:
on-the-fly projection does not work properly in QGIS and gvSIG...
Whatever I do in GRASS, I always get correct placement of the GPS data.

This seems to be QGIS bug # 1079
  http://trac.osgeo.org/qgis/ticket/1079

(usual suspects reporting in :slight_smile:

Using on-the-fly projection in the other programs, I get a
difference of about 100m, so I guess it's a datum issue as
Hamish suggests.

It would seem so. I can reproduce it here (see bug report). In my case
looking in the layer's "P"roperies projection tab I see the the +datum=
and/or +towgs84= terms are missing* versus what is found in the EPSG
file.

[*] file: /usr/share/qgis/resources/srs.db (SQLite DB)

Root cause seems to be that QGIS lacks a "pick appropriate datum transform
parameters to use" GUI in the case when there are multiple possibilities,
and so it doesn't try any datum transform at all:
http://www.nabble.com/why-is-EPSG%3A2180-not-recognised--td16505409.html#a16592790

Thanks for following up on this.

I don't have the feeling that either 4326 or 31370 offer multiple
possibilities for transform parameters, but don't know for 100% sure. But
the way I read the above nabble thread (and especially Frank's quote) the
issue is that some projections are missing in the proj and thus the QGIS
database because of the multiple transform possibilities (so when you use
the default proj epsg file in GRASS the same problem probably applies).
However both of the projections _are_ in the database and QGIS recognizes
them.

So, I'm not sure whether this is an issue with the epsg codes, or rather a
problem with the transformation code.

But I think it's time to take this over to the relevant QGIS and gvSIG
(which has the same problem) mailing lists...

Moritz

14

Jonathan Aguero-Valverde wrote:

DBMI-Postgres driver error:
Cannot connect to Postgres: could not connect to server:
Connection refused
(0x0000274D/10061)
        Is the server running on host "localhost"
and accepting
        TCP/IP connections on port 6543?

GRASS_INFO_WARNING(4116,1): Cannot open database
'host=localhost,dbname=thesisWA,port=6543'
GRASS_INFO_END(4116,1)

GRASS_INFO_ERROR(4116,2): Unable to open database
<host=localhost,dbname=thesisWA,port=6543> by driver
<pg>
GRASS_INFO_END(4116,2)

can you use any other vector processing modules which use the DB?
disconnect from the internet and try turning off your firewall temporarily?

It looks like somehow the port is hardcoded to 6543 in
v.distance.

this is not the case.

If I tried db.tables -p it gives me the list of my tables so it
looks like my connection settings are fine.

to say that for sure you would need to see that other processing modules work, not just deescriptive ones. v.db.addcol?

Hamish

15

Moritz Lennert:

I don't have the feeling that either 4326 or 31370
offer multiple possibilities for transform parameters,
but don't know for 100% sure.

no, they don't. LL WGS84 is just that (as long as we consider it the root coordinate system), and 31370 explicitly gives 7 term +towgs84= params AFAICT. AFAIU the ambiguity of which transform params to use comes only when +datum= is given without the hint of which path to take to get there.

But the way I read the above nabble thread (and especially
Frank's quote) the issue is that some projections are missing
in the proj and thus the QGIS database because of the multiple
transform possibilities

strictly speaking the "projection" (lcc) is there and fine, but (AFAIU) the datum transform params for a given datum are left out if there are more than one option in the upstream EPSG database. But the +datum= name should still be given in that case.

but I think that is a misdirection: in this case I think QGIS's srs.db is broken, missing *all* datum info not just the ambiguos transform params. The issue Frank raises still applies, but I do not think it is the main problem here.

(so when you use the default proj epsg file in GRASS the same
problem probably applies).

No, in those cases GRASS prompts the user to pick the appropriate method from a list, with the default at the top of the list.

the PROJ epsg file is not broken in this sense, only limited in the information it provides. the QGIS srs.db file is broken. FWIW GRASS uses its own $GISBASE/etc/datum*, etc/proj*, and etc/ogr_csv/* files, the PROJ epsg file is only used in the startup GUI location creation wizard (we should change that so they all use the same). The ogr_csv files are sync'd (typically by Markus) with the latest GDAL versions prior to each release or after each new upstream EPSG DB release. Note AFAIU the current PROJ epsg file awaits regeneration for minor FP tweaks arising from the latest atof() GDAL issue commits. string comparison of the scaling factor etc used by QGIS are no longer equal...

However both of the projections _are_ in the database and
QGIS recognizes them.

but they are incomplete! compare the entry for epsg:31370 in /usr/share/proj/epsg and in /usr/share/qgis/resources/srs.sb (using sqlitebrowser or such). you will see the srs.db version is missing the +towgs84= part. (epsg:31370 is srs:1817 there, or just go to the layers proj properties in QGIS and look)

So, I'm not sure whether this is an issue with the epsg
codes, or rather a problem with the transformation code.

incorrectly generated QGIS srs.db file I think.

But I think it's time to take this over to the relevant
QGIS and gvSIG (which has the same problem) mailing lists...

qgis bugs:
http://trac.osgeo.org/qgis/ticket/1035
http://trac.osgeo.org/qgis/ticket/1037
http://trac.osgeo.org/qgis/ticket/1079

Could you could create/find the gvSIG report?

whole lotta qualifiers,
Hamish

16

Sort of off-topic:

I like that GRASS is strict with the projection of data. You are 100%
sure that all data are there where you see them. I don't think it's
really necessary to add an on-the-fly function in GRASS. You can always
re-project and know what happens.

Nikos

Me too. I like being forced to get all my "ducks in a row", so to speak.
There's never any ambiguity about what projection your data is in, or if you
are witnessing on-the-fly projection errors or not.

~ Eric.

17

On Sun, 25 May 2008, Hamish wrote:

But the way I read the above nabble thread (and especially
Frank's quote) the issue is that some projections are missing
in the proj and thus the QGIS database because of the multiple
transform possibilities

strictly speaking the "projection" (lcc) is there and fine, but (AFAIU) the datum transform params for a given datum are left out if there are more than one option in the upstream EPSG database. But the +datum= name should still be given in that case.

There are actually only a few +datum= codes recognised by PROJ.4 - most of them are specific to GRASS. So it isn't meaningful to include a +datum=xxxx if the PROJ string is destined to be used by something other than GRASS. The OGR functions that GRASS uses to interpret EPSG codes somehow manage to return the full EPSG name for the datum, which is then compared to GRASS's datum.table (the second field in this table holds the full EPSG name). Combined with datumtransform.table that gives GRASS enough information to go on to present the user with a choice of datum transform paramteters. However the full EPSG name isn't included in the PROJ.4 string, and if that is all QGIS is using then maybe that's where its deficiency stems from. It is all very messy and complicated.

Paul

18

Moritz Lennert ha scritto:

I agree 100%. I'm not pleading for on-the-fly projection at all. However,
since I'm currently teaching a training course which also includes QGIS
and gvSIG, I also use that feature, but apparently I'd better not !

agreed: same experience for me, here in Evora.
pc

19

I have finally come around to looking at this issue again...

On Mon, May 26, 2008 05:01, Hamish wrote:

Moritz Lennert:
strictly speaking the "projection" (lcc) is there and fine, but (AFAIU)
the datum transform params for a given datum are left out if there are
more than one option in the upstream EPSG database. But the +datum= name
should still be given in that case.

but I think that is a misdirection: in this case I think QGIS's srs.db is
broken, missing *all* datum info not just the ambiguos transform params.
The issue Frank raises still applies, but I do not think it is the main
problem here.

[...]

However both of the projections _are_ in the database and
QGIS recognizes them.

but they are incomplete! compare the entry for epsg:31370 in
/usr/share/proj/epsg and in /usr/share/qgis/resources/srs.sb (using
sqlitebrowser or such). you will see the srs.db version is missing the
+towgs84= part. (epsg:31370 is srs:1817 there, or just go to the layers
proj properties in QGIS and look)

I have been able to solve the problem in QGIS by creating a custom
projection including the transform info (copy-paste of the g.proj -jf
output, replacing +a and +rf with the relevant +ellps). Now the result
looks fine.

Could you could create/find the gvSIG report?

Actually, after further investigation, I found that gvSIG allows for
transformations. When you load a file, you can give its projection info.
See [1]. Notice the option to choose a transformation at the bottom (in
this case I chose the EPSG list). You then have a series of choices of
transformations [2]. Or you can chose to go for manual and enter the
transformation parameters yourself [3].

So it's finally only QGIS which has the problem and you can sail around it
by creating a custom projection containing the transformation info.

Moritz

[1] http://geog-pc40.ulb.ac.be/grass/misc/transformation_gvSIG_1.png
[2] http://geog-pc40.ulb.ac.be/grass/misc/transformation_gvSIG_2.png
[3] http://geog-pc40.ulb.ac.be/grass/misc/transformation_gvSIG_3.png

20

Hamish pisze:

but they are incomplete! compare the entry for epsg:31370 in /usr/share/proj/epsg and in /usr/share/qgis/resources/srs.sb (using sqlitebrowser or such). you will see the srs.db version is missing the +towgs84= part.

Definitely, it is a defect in QGIS srs.db. I have posted a script to
regenerate a correct, updated srs.db to QGIS dev ML [1]. Among the other
fixes, the database the script generates has the EPSG:31370 +towgs84
definition complete too. Once the remaining issues outlined in [1] are
solved, there are chances that the fixed srs.db will be included in next
QGIS release.

[1]http://www.nabble.com/forum/ViewPost.jtp?post=17586603&framed=y

--
Maciej Sieczka
www.sieczka.org