[OSRS-PROJ] Converting to Finnish National Grid Coordinate System (KKJ)

[Seppo Enarvi]

Hello,

I got my hands on a data cable for my Garmin GPS receiver a couple of
weeks ago and wrote a program that reads the position from the GPS device
and displays it on a GeoTIFF map. Everything I know about coordinate
transformations is what I have read from the Internet during these weeks.
Currently, my program uses

- GARNIX / LIBGEO [1] to retrieve the position and convert it to the
Finnish grid coordinate system (KKJ)

- libgeotiff / libtiff [2] to display the map and my position on the map

While trying to find another coordinate conversion library that could read
all the parameters required for the conversion from an EPSG database [3],
I run into the PROJ.4 library.


The KKJ grid coordinate system and the parameters I have used are briefly
described in this document:

http://www.nls.fi/kartta/julkaisu/kkj.html

As you can see, the coordinate system is based on ED50 coordinates, but
after doing the Transverse Mercator projection, the projected coordinates
should be transformed with 2-dimensional Helmert transformation. I think
that many programs ignore the final 2-dimensional transformation,
resulting in poor accuracy. I'm not sure if there are any parameters in
PROJ.4 that could be used to specify the 2-dimensional Helmert
transformation.


Here's what I have tried:

1) Change the EPSG database to include the ED50 datum shift:

# KKJ / Finland Uniform Coordinate System
<2393> +proj=tmerc +lat_0=0.000000000 +lon_0=27.000000000 +k=1.000000
+x_0=3500000.000 +y_0=0.000 +ellps=intl
+towgs84=-93.477,-103.453,-123.431,-0.246,0.109,0.068,-2.062 +units=m
no_defs <>


2) Run cs2cs:

>cs2cs -v +proj=latlong +datum=WGS84 +to +init=epsg:2393 coordinates
# ---- From Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=latlong +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# ---- To Coordinate System ----
#Transverse Mercator
#       Cyl, Sph&Ell
# +init=epsg:2393 +proj=tmerc +lat_0=0.000000000 +lon_0=27.000000000
# +k=1.000000 +x_0=3500000.000 +y_0=0.000 +ellps=intl
# +towgs84=-93.477,-103.453,-123.431,-0.246,0.109,0. +units=m +no_defs
#--- following specified but NOT used
# +068,-2.062
3531710.65      7074338.84 -8.973


3) Do the 2D Helmert manually:

		x = 7074338.863;
		y = 3531709.713;

		y -= 3000000;

		double a = 1.0000007515;
		double b = 0.0000043933;
		double dx = -61.5805;
		double dy = 95.6691;

		x2 = (a * x) - (b * y) + dx;
		y2 = (b * x) + (a * y) + dy;

		y2 += 3000000;


This gives me exactly the correct result (I checked from the example at
http://www.nls.fi/kartta/julkaisu/kkj.html ), except that

1) For some reason PROJ.4 truncates the towgs84 parameter to 41
characters(!), and more importantly,

2) I don't want to do the 2D Helmert transformation manually.


What would you suggest?


Btw, I think it's quite amazing that cs2cs gives wrong results if I enter
latitude and longitude in wrong order even when I specify which one is
which (with N and E). It could at least display some warning message.


Thanks,

Seppo Enarvi



[1]
http://homepage.ntlworld.com/anton.helm/garnix.html

[2]
http://remotesensing.org/geotiff/geotiff.html
http://www.libtiff.org/

[3]
The projection zone that a GeoTIFF file represents is stored in the file
as an ID to the EPSG geodesy parameter database available at
http://www.epsg.org/

----------------------------------------
PROJ.4 Discussion List
See http://www.remotesensing.org/proj for subscription, unsubscription
and other information.