[GRASSLIST:3818] Projection lat/lon (gps) to local system

Morten_Sickel1 · July 6, 2004, 9:20pm

Hi,

I was using grass a bit a couple of years ago, but I have been away for a
while. Now I'm back and have started looking at 5.7, quite impressive, kudos
to the developers.

But..

I need to project some gps data to use them together with some raster maps in
a local (Norwegian) projection :

GRASS 5.7.cvs:~/NinaSather/Kart >g.proj -p
-PROJ_INFO-------------------------------------------------
name : Transverse Mercator
datum : ngo48
towgs84 : 278.3,93.0,474.5,7.889,0.050,-6.610,6.21
proj : tmerc
ellps : modif_bessel
a : 6377492.0175999999
es : 0.0066743722
f : 299.1528128500
lat_0 : 58.0000000000
lon_0 : 8.3895830000
k_0 : 1.0000000000
x_0 : 0.0000000000
y_0 : 0.0000000000

As fas as I can see, those values are correct, but anyhow, the gps data ends
up about 300 meters off from where they should be (see
http://sickel.net/gps.html for an illustration) It looks like some kind of
datum / ellipse error, but I have no idea how to fix that (Except for
identifying a few pass points and write a transformation script) I have a
table of dx,dy,.dz values for my local projection and I have a feeling that
those are somehow connected to the towgs84 parameter, but I don't see how to
set those.

Btw, if anybody among you are reading 'Landscape and Urban planning' the maps
in
Hanne Sickel, Margareta Ihse, Ann Norderhaug and Morten A. K. Sickel: How to
monitor semi-natural key habitats in relation to grazing preferences of
cattle in mountain summer farming areas: An aerial photo and GPS method study
Landscape and Urban Planning, Volume 67, Issues 1-4, 15 March 2004, Pages
67-77
are made using grass. Each map made using two outputs using ps.map (5.0.0) and
combined in the gimp, making the upper layer semi-transparent. Feel free to
contact me if anybody have any questions regarding the production line.

--
Morten Sickel
http://sickel.net/
Drøbak, Norway

Paul_Kelly · July 6, 2004, 11:30pm

Hello Morten

Morten Sickel wrote:

[...]

I need to project some gps data to use them together with some raster maps in
a local (Norwegian) projection :

GRASS 5.7.cvs:~/NinaSather/Kart >g.proj -p
-PROJ_INFO-------------------------------------------------
name : Transverse Mercator
datum : ngo48
towgs84 : 278.3,93.0,474.5,7.889,0.050,-6.610,6.21
proj : tmerc
ellps : modif_bessel
a : 6377492.0175999999
es : 0.0066743722
f : 299.1528128500
lat_0 : 58.0000000000
lon_0 : 8.3895830000
k_0 : 1.0000000000
x_0 : 0.0000000000
y_0 : 0.0000000000

As fas as I can see, those values are correct, but anyhow, the gps data ends
up about 300 meters off from where they should be (see

Please describe how you re-projected your GPS points into the ngo48
datum so we can spot if there is anything wrong there---this seems to be
where the problem is.

http://sickel.net/gps.html for an illustration) It looks like some kind of
datum / ellipse error, but I have no idea how to fix that (Except for
identifying a few pass points and write a transformation script) I have a
table of dx,dy,.dz values for my local projection and I have a feeling that
those are somehow connected to the towgs84 parameter, but I don't see how to
set those.

Edit the PROJ_INFO file, and either change it to have a line
towgs84: dx,dy,dz
or
dx:
dy:
dz:
with the 3 numbers.

But don't do this until you get the datum transformation working because
there is obviously some problem there. Once you have this done you can
compare your values with the general ones for Norway g.setproj has put
in and see which one makes your data line up better.

Paul

Morten_Sickel1 · July 7, 2004, 8:12am

Thanks Paul,

On Wednesday 07 July 2004 01:30, Paul Kelly wrote:

Please describe how you re-projected your GPS points into the ngo48
datum so we can spot if there is anything wrong there---this seems to be
where the problem is.

Are there more ways to do it? (I know it used to be some kind of 'sites'
projection utililty, but I didn't find that in 5.7) I just made a lat/lon
location into which I imported the gps-points, and then I used v.proj to
reproject them to my 'local' location.

Morten

--
Morten Sickel
http://sickel.net/
Drøbak, Norway

Morten_Sickel1 · July 7, 2004, 1:04pm

Thanks again, Paul, I did believe that I had given all relevant info, but
here you are:

PROJ_INFO in the local data set:

GRASS 5.7.cvs:~/NinaSather/Kart > \
cat /home/grassdata/ninahs/PERMANENT/PROJ_INFO
name: Transverse Mercator
datum: ngo48
towgs84: 278.3,93.0,474.5,7.889,0.050,-6.610,6.21
proj: tmerc
ellps: modif_bessel
a: 6377492.0175999999
es: 0.0066743722
f: 299.1528128500
lat_0: 58.0000000000
lon_0: 8.3895830000
k_0: 1.0000000000
x_0: 0.0000000000
y_0: 0.0000000000

and in the from-dataset:
GRASS 5.7.cvs:~/NinaSather/Kart > \
cat /home/grassdata/ninahsgps/PERMANENT/PROJ_INFO
name: Latitude-Longitude
proj: ll
ellps: wgs84

(Should there be a datum here as well?)

the command-line:

GRASS 5.7.cvs:~/NinaSather/Kart > v.proj input=skogstad070702 \
location=ninahsgps

The gps-data were originating from a Magellan GPS315 with files like:
GRASS 5.7.cvs:~/NinaSather/Kart > head ../gps/071002Skogstad.txt
Datum,WGS84,WGS84
TP,DM, 61.03175, 008.51871,07/09/2002,14:57:37 ,1,0000000M,A,TRAK01
TP,DM, 61.03178, 008.51823,07/09/2002,14:57:37 ,0,0000000M,A,TRAK01
TP,DM, 61.03180, 008.51777,07/09/2002,14:57:37 ,0,0000000M,A,TRAK01
TP,DM, 61.03193, 008.51710,07/09/2002,14:57:37 ,0,0000000M,A,TRAK01
TP,DM, 61.03210, 008.51659,07/09/2002,14:57:37 ,0,0000000M,A,TRAK01
TP,DM, 61.03214, 008.51709,07/09/2002,14:57:37 ,0,0000000M,A,TRAK01
TP,DM, 61.03174, 008.51708,07/09/2002,14:57:37 ,0,0000000M,A,TRAK01
TP,DM, 61.03154, 008.51724,07/09/2002,14:57:37 ,0,0000000M,A,TRAK01
TP,DM, 61.03145, 008.51771,07/09/2002,14:57:37 ,0,0000000M,A,TRAK01

And have been piped through:
GRASS 5.7.cvs:~/NinaSather/Kart > cat ~/perl/gps2sites.pl
#!/usr/bin/perl -w
sub deg2dec{
  my $in = shift;
  $deg=int($in);
  $in -= $deg;
  $in /=0.6;
  $in +=$deg;
}

my $i;
while(<>){
  my @fields=split /,\s*/;
  next if $#fields<3;
  next if $field[1] ne 'DM';
  local $,=' ';
  print deg2dec($fields[3]),deg2dec($fields[2]),++$i,"\n";
}

(the DM in the GPS-file shows that the data are on the form degrees and
decimal minutes)

I cannot recall any other data treatment that should be of significance., but
wold love to be corrected.

In case that should matter, my versjon of proj is
GRASS 5.7.cvs:~/NinaSather/Kart > proj
Rel. 4.4.7, 31 March 2003
usage: proj [ -beEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

Well, one possible problem that I suddeny recognice, (I do not believe it has
happended, but you newer know) is that the data-collectors might have managed
to change the datum of the GPS, as you can see from the dates in the file,
the data have been "maturing" for a while, so it is not possible to see if
that is the case, but I do believe that the GPS in all cases outputs data
reffered to WGS84.

Morten

--
Morten Sickel
http://sickel.net/
Drøbak, Norway

Paul_Kelly · July 7, 2004, 2:38pm

Hello Morten

Morten Sickel wrote:

[...]

and in the from-dataset:
GRASS 5.7.cvs:~/NinaSather/Kart > \
cat /home/grassdata/ninahsgps/PERMANENT/PROJ_INFO
name: Latitude-Longitude
proj: ll
ellps: wgs84

(Should there be a datum here as well?)

Yes that is probably your problem. Both locations need to have datum
information in them or the datum transformation will not take place. If
you add a line
datum: wgs84
to the end of this file /home/grassdata/ninahsgps/PERMANENT/PROJ_INFO
then it should work.

Please let us know how you get on.

Paul

Morten_Sickel1 · July 7, 2004, 2:56pm

On Wednesday 07 July 2004 16:38, Paul Kelly wrote:

Morten Sickel wrote:
[...]
> GRASS 5.7.cvs:~/NinaSather/Kart > \
> cat /home/grassdata/ninahsgps/PERMANENT/PROJ_INFO
> name: Latitude-Longitude
> proj: ll
> ellps: wgs84
>
> (Should there be a datum here as well?)

Yes that is probably your problem. Both locations need to have datum
information in them or the datum transformation will not take place. If
you add a line
datum: wgs84
to the end of this file /home/grassdata/ninahsgps/PERMANENT/PROJ_INFO
then it should work.

Please let us know how you get on.

That did it. Works as close to perfectly as the gpssystem goes! Thanks a lot!

A proposal to the developers: Is it possible to issue a warning if one of the
locations (involved in a projection) has a datum, whereas the other has not?
e.g in my case a warning like:

No datum found in ninahsgps, the datum information in ninahs will not be used.

or maybe set wgs84 as a default datum in a ll-location. Even if I answer wgs84
when asked for a datum when setting up a ll-location, only the ellipsoide is
set.

(using 5.7)

Morten

--
Morten Sickel
http://sickel.net/
Drøbak, Norway

H_B · July 9, 2004, 4:37am

> Please describe how you re-projected your GPS points into the ngo48
> datum so we can spot if there is anything wrong there---this seems
> to be where the problem is.
>
Are there more ways to do it? (I know it used to be some kind of
'sites' projection utililty, but I didn't find that in 5.7) I just
made a lat/lon location into which I imported the gps-points, and then
I used v.proj to reproject them to my 'local' location.

GRASS 5.3: m.proj2 -i

GRASS 5.7: cs2cs + `g.proj -jf`
(cs2cs is from PROJ4)

see: http://article.gmane.org/gmane.comp.gis.grass.devel/4013

also in 5.3, s.in.garmin for Garmin GPSs will automatically project +
transform sites. I haven't finished updating v.in.garmin yet for
5.3 or 5.7 but they're mostly there..

Hamish