Geo Coordinate Converter

About

gk-slo is a converter between geographic cartesian coordinates (Gauss-Krueger/D48, Transverse Mercator/D96) and geodetic coordinates (latitude/longitude on ETRS89/WGS84) for Slovenia. It can be used as a replacement for the official conversion program SiTra (with Helmert parameters for the whole Slovenia, no regional parameters) or more precise when used with the built-in affine/triangle-based transformation using the reference virtual tie points v3.0 (for detailed description see AFT), which is similar as the official conversion program 3Tra.

GCC Cover Image

Program can read files in SiTraNet format (ASCII XYZ), LIDAR (ASCII XYZ with semicolon, .asc) or ESRI shapefile (ArcGIS .shp format, use gk-shp).

The following transformations are available (in both directions):

1.	xy (D96/TM)	⇔	φλ (ETRS89)
2.	xy (D48/GK)	⇔	xy (D96/TM)	Helmert transformation
3.	xy (D48/GK)	⇔	φλ (ETRS89)	Helmert transformation
4.	xy (D48/GK)	⇔	xy (D96/TM)	Affine transformation
5.	xy (D48/GK)	⇔	φλ (ETRS89)	Affine transformation

For calculating heights with the help of geoid model two absolute geoid models for Slovenia are available: Slo2000 and EGM2008.

It's written in C language and can be compiled and used on all major operating systems. Coordinate conversion routines can be easily adapted to locations other than Slovenia (via definition of ellipsoid, projection and Helmert parameters).

Detailed description of coordinate conversion routines and their API (in module "geo.c") is in file geo_api.md.

Comments

If you have any comments or suggestions, please leave them on a Geo Coordinate Converter's Facebook page

Usage

$ gk-slo [<options>] [<inpname> ...]
  -d                enable debug output
  -ht               calculate output height with 7-params Helmert trans.
  -hc               copy input height unchanged to output
  -hg               calculate output height from geoid model (default)
  -g slo|egm        select geoid model (Slo2000 or EGM2008)
                    default: Slo2000
  -dms              display fila in DMS format after height
  -t <n>            select transformation:
                     1: xy   (d96tm)  --> fila (etrs89), hg?, default
                     2: fila (etrs89) --> xy   (d96tm),  hg
                     3: xy   (d48gk)  --> fila (etrs89), ht
                     4: fila (etrs89) --> xy   (d48gk),  hg
                     5: xy   (d48gk)  --> xy   (d96tm),  hg(hc)
                     6: xy   (d96tm)  --> xy   (d48gk),  ht(hc)
                     7: xy   (d48gk)  --> xy   (d96tm),  hc, affine trans.
                     8: xy   (d96tm)  --> xy   (d48gk),  hc, affine trans.
                     9: xy   (d48gk)  --> fila (etrs89), hg, affine trans.
                    10: fila (etrs89) --> xy   (d48gk),  hg, affine trans.
  -r                reverse parsing order of xy/fila
                    (warning is displayed if y < 200000 or la > 17.0)
  <inpname>         parse and convert input data from <inpname>
                    <inpname> "-" means stdin, use "--" before
  -o -|=|<outname>  write output data to:
                    -: stdout (default)
                    =: append ".out" to each <inpname> and
                       write output to these separate files
                    <outname>: write all output to 1 file <outname>

Typical input data format (SiTra .xyz or LIDAR .asc):
[<label> ]<fi|x> <la|y> <h|H>
[<label>;]<fi|x>;<la|y>;<h|H>

A new program gk-shp, able to read ESRI shapefiles (ArcGIS .shp format), has similar syntax:

$ gk-shp [<options>] <inpname> <outname>
  -d                enable debug output
  -ht               calculate output height with 7-params Helmert trans.
  -hc               copy input height unchanged to output
  -hg               calculate output height from geoid model (default)
  -g slo|egm        select geoid model (Slo2000 or EGM2008)
                    default: Slo2000
  -t <n>            select transformation:
                     1: xy   (d96tm)  --> fila (etrs89), hg?, default
                     2: fila (etrs89) --> xy   (d96tm),  hg
                     3: xy   (d48gk)  --> fila (etrs89), ht
                     4: fila (etrs89) --> xy   (d48gk),  hg
                     5: xy   (d48gk)  --> xy   (d96tm),  hg(hc)
                     6: xy   (d96tm)  --> xy   (d48gk),  ht(hc)
                     7: xy   (d48gk)  --> xy   (d96tm),  hc, affine trans.
                     8: xy   (d96tm)  --> xy   (d48gk),  hc, affine trans.
                     9: xy   (d48gk)  --> fila (etrs89), hg, affine trans.
                    10: fila (etrs89) --> xy   (d48gk),  hg, affine trans.
  -r                reverse parsing order of xy/fila
                    (warning is displayed if y < 200000 or la > 17.0)
  <inpname>         parse and convert input data from <inpname>
  <outname>         write output data to <outname>

Input data format:
ESRI Shapefile (ArcGIS)

Example 1 (D48)

Input file VTG2225.XYZ (DMV, in SiTraNet format, Gauss-Krueger/D48):

0000001 509000.000 76000.000 343.30
0000002 509005.000 76000.000 342.80
0000003 509010.000 76000.000 342.30

Convert to new coordinate system (Transverse Mercator/D96); heights should be copied, not calculated:

$ gk-slo -t 5 -hc VTG2225.XYZ
VTG2225.XYZ: possibly reversed x/y
0000001 509487.490 575640.546 343.300
0000002 509492.490 575640.546 342.800
0000003 509497.490 575640.546 342.300

If you see the warning "possibly reversed x/y", use the option "-r" to get the correct conversion:

$ gk-slo -t 5 -hc -r VTG2225.XYZ
0000001 76484.893 508628.990 343.300
0000002 76484.893 508633.991 342.800
0000003 76484.893 508638.991 342.300

Convert the same file to ETRS89/WGS84 coordinates. This time the height will be recalculated using Helmert transformation (default):

$ gk-slo -t 3 -r VTG2225.XYZ
0000001 45.8281655853 15.1110624092 389.063
0000002 45.8281655218 15.1111267639 388.563
0000003 45.8281654582 15.1111911187 388.063

For better readout you can include the option "-dms":

$ gk-slo -t 3 -r -dms VTG2225.XYZ
0000001 45.8281655853 15.1110624092 389.063 45 49 41.39611 15  6 39.82467
0000002 45.8281655218 15.1111267639 388.563 45 49 41.39588 15  6 40.05635
0000003 45.8281654582 15.1111911187 388.063 45 49 41.39565 15  6 40.28803

Store result in a file:

$ gk-slo -t 3 -r VTG2225.XYZ -o VTG2225.flh
(creates file VTG2225.flh)

Example 2 (D96)

Input file VTC0512.XYZ (DMV, in SiTraNet format, Transverse Mercator/D96):

0000001 412250 97000 606.2
0000002 412250 96995 606.9
0000003 412250 96990 607.9

Convert to ETRS89/WGS84 coordinates, height will be recalculated using geoid model Slo2000 (default):

$ gk-slo -t 1 -r VTC0512.XYZ
0000001 46.0071929697 13.8669428837 652.772
0000002 46.0071479903 13.8669438021 653.472
0000003 46.0071030110 13.8669447206 654.472

If you want to use the EGM2008 geoid model, use the "-g egm" option:

$ gk-slo -t 1 -r -g egm VTC0512.XYZ
0000001 46.0071929697 13.8669428837 652.660
0000002 46.0071479903 13.8669438021 653.359
0000003 46.0071030110 13.8669447206 654.359

Example 3 (ETRS89/WGS84)

Convert ETRS89/WGS84 coordinates to Transverse Mercator/D96 via keyboard (ignoring height, use "--" to stop parsing options):

$ gk-slo -t 2 -- -
46.0071929697 13.8669428837 0 <Enter>
97000.000 412250.000 -46.572

Convert file VTG2225.flh (with ETRS89/WGS84 coordinates, see Example 1) to Gauss-Krueger/D48:

$ gk-slo -t 4 VTG2225.flh
0000001 76000.000 509000.000 342.896
0000002 76000.000 509005.000 342.396
0000003 76000.000 509010.000 341.896

If you compare the output with the original VTG2225.XYZ, you'll notice a small difference in heights. This is because the heights in VTG2225.flh were calculated using Helmert transformation but the default height calculation for "-t 4" is using geoid model. Which height calculation you use for conversion depends also on input data. There are some recommended defaults for each type of conversion (see Usage).

Example 4 (processing many files)

If gk-slo was compiled with MinGW on Windows or is being used on Unix, you can process many files with one command. Input files ∗.XYZ (DMV, Transverse Mercator/D96):

VTH0720.XYZ
VTH0721.XYZ
VTH0722.XYZ
...

Convert all files to ETRS89/WGS84 coordinates (with some debug info to see what's going on):

$ gk-slo -t 1 -r -d ∗.XYZ -o =
Processing VTH0720.xyz
Processing time: 4.854913
Processing VTH0721.xyz
Processing time: 4.846438
Processing VTH0722.xyz
Processing time: 4.846361
...

Results of conversion for each file are written to a new file with an extension ".out":

VTH0720.XYZ.out
VTH0721.XYZ.out
VTH0722.XYZ.out
...

Example 5 (ESRI shapefiles)

Input file RABA_20151031.shp (GERK, in ESRI shapefile format, Gauss-Krueger/D48), convert to ETRS89/WGS84 coordinates using affine transformation (with debug info):

$ gk-shp -t 9 -dd RABA_20151031.shp raba_conv.shp
Processing RABA_20151031.shp
Shapefile type: Polygon, number of shapes: 1601832
Shape: 678 (0.04%) ...

Result of conversion is a set of files according to ESRI shapefile format:

raba_conv.cpg
raba_conv.dbf
raba_conv.prj
raba_conv.shp
raba_conv.shx

In file raba_conv.prj an output projection (WGS84) is stored, so converted files can be easily opened by GIS programs.

Affine/Triangle-based Transformation

Theory

Affine transformation only works with cartesian coordinates (GK x,y ⇔ TM x,y). You need minimum 3 control points (for example a triangle), each with 4 coordinates (Xs, Ys, Xt, Yt), where Xs,Ys are coordinates in the source system and Xt,Yt coordinates in the target system.

For each new point Xi,Yi transformation is calculated via the following equation:

Xo = a∗Xi + b∗Yi + c
Yo = d∗Xi + e∗Yi + f

where parameters a..f have the following meaning:

a: Scale X
e: Scale Y
d: Rotation X
b: Rotation Y
c: Translation X
f: Translation Y

To calculate parameters a..f for a given triangle you have to solve the system of linear equations with 6 unknowns, represented with the following matrices:

| Xs1 Ys1 1 0 0 0 | | a |   | Xt1 |
| Xs2 Ys2 1 0 0 0 | | b |   | Xt2 |
| Xs3 Ys3 1 0 0 0 | | c | = | Xt3 |
| 0 0 0 Xs1 Ys1 1 | | d |   | Yt1 |
| 0 0 0 Xs2 Ys2 1 | | e |   | Yt2 |
| 0 0 0 Xs3 Ys3 1 | | f |   | Yt3 |

This system can be solved using elimination method (row operations ⇒ identity matrix) or with other methods (Jacobi, Gauss-Seidel, ...).

Pre-calculation

For Slovenia we have 899 reference virtual tie points (v3.0) in both coordinate systems (GK and TM), which form 1776 triangles:

We can pre-calculate parameters a..f for every triangle and store their values in a table to be directly used in a program.

To calculate triangles from the initial list of reference points we use the existing program triangle from Jonathan Shewchuk (UCB), which generates exact Delaunay triangulations.

Exact description of process

get the reference virtual tie points (already supplied):
⇒ file virtualne_vezne_tocke_v3.0.txt
extract GK and TM coordinates from the reference file
use triangle program to create triangles
use the supplied program ctt to create tables to be included in a program:
⇒ files aft_gktm.h and aft_tmgk.h
ctt solves the above mentioned system of linear equations for each triangle.

All this can be automated using the supplied scripts cvvt.sh (for Unix) or cvvt.bat (for Windows).

Usage in a program

For each cartesian coordinate you have to loop through the array of 1776 pre-calculated triangles and check if the coordinate is contained in the triangle. This check should be very accurate, otherwise you'll end up with a point, which lies somewhere "between" the triangles (see geo_api.md).

When/if such triangle is found, a simple transformation using parameters a..f is applied (see Theory above).