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.
$ 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)
0000001 509000.000 76000.000 343.30 0000002 509005.000 76000.000 342.80 0000003 509010.000 76000.000 342.30Convert 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.300If 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.300Convert 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.063For 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.28803Store result in a file:
$ gk-slo -t 3 -r VTG2225.XYZ -o VTG2225.flh (creates file VTG2225.flh)
0000001 412250 97000 606.2 0000002 412250 96995 606.9 0000003 412250 96990 607.9Convert 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.472If 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
$ gk-slo -t 2 -- - 46.0071929697 13.8669428837 0 <Enter> 97000.000 412250.000 -46.572Convert 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.896If 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).
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 ...
$ 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.shxIn file raba_conv.prj an output projection (WGS84) is stored, so converted files can be easily opened by GIS programs.
(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 + fwhere parameters
a..f
have the following meaning:
a: Scale X e: Scale Y d: Rotation X b: Rotation Y c: Translation X f: Translation YTo 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, ...).
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.
virtualne_vezne_tocke_v3.0.txt
ctt
to create tables to be included in a program: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).
When/if such triangle is found, a simple transformation using parameters
a..f
is applied (see Theory above).
Last modified: January 10, 2020 Matjaz Rihtar <matjaz@eunet.si> | Copyright © 2014-2020 Matjaz Rihtar |