Bezugsrahmenwechsel: ST_Fineltra in Action

04 October 2015

Sandro Santilli hat eine erste Version der ST_Fineltra Funktion für PostGIS zum Testen freigegeben. Um die Funktion verwenden zu können, sind zwei Schritte notwendig:

Installieren der neuen Funktion
Importieren der Dreiecksvermaschung in die Datenbank

Das Installieren der Funktion funktioniert unter Ubuntu 14.04 problemlos. Hat man sowieso bereits die üblichen «GIS-Tools» (PostgreSQL/Postgis, gdal/ogr etc.) installiert, fehlt anscheinend nur das dev-Paket von liblwgeom:

sudo apt-get install liblwgeom-dev

Anschliessend kann gemäss der Installationsanleitung im README vorgegangen werden:

git clone https://github.com/strk/fineltra.git
cd fineltra
./autogen.sh
./configure
make
sudo make install

Die Funktion ist eine PostgreSQL-Extension und kann für die gewünschte Datenbank aktiviert werden:

sudo -u postgres psql -d xanadu2 -c "CREATE EXTENSION fineltra;"

Der zweite Schritt ist der Import der Dreiecksvermaschung. Eine Spatialite-Datei mit einer Tabelle mit sämtlichen Dreiecken in beiden Bezugsrahmen gibts hier. Vor Jahren habe ich diese aus den Shapedateien von Swisstopo erstellt. Der Import in die Datenbank ist ein einfacher ogr2ogr-Befehl. Leider scheint das mit einer Version kleiner 2.0 nicht wirklich gut zu funktionieren, weil die Tabelle zwei Geometriespalten aufweist. Bei mir werden nur die Sachattribute importiert. Mit ogr2ogr >= 2.0 reicht aber folgender Befehl:

ogr2ogr -f "PostgreSQL" PG:"dbname='xanadu2' host='localhost' port='5432' user='stefan' password='ziegler12'" chenyx06.sqlite chenyx06 -lco SCHEMA=av_chenyx06 -nln chenyx06_triangles

Die Tabelle mit den Dreiecken muss beide Bezugsrahmen beinhalten. Es ist nicht möglich die Funktion ST_Fineltra mit Dreiecksgeometrien aus verschiedenen Tabelle zu bedienen.

Das Schema av_chenyx06 muss bereits existieren. Ansonsten muss es manuell angelegt werden:

sudo -u postgres psql -d xanadu2 -c "CREATE SCHEMA av_chenyx06 AUTHORIZATION stefan;"

Es ist unbedingt darauf zu achten, dass die Geometrietabellen einen Geometrieindex aufweisen. Bei mir werden diese beim Import mit ogr2ogr automatisch erzeugt. Falls dies nicht der Fall ist:

sudo -u postgres psql -d xanadu2 -c "CREATE INDEX ON av_chenyx06.chenyx06_triangles USING GiST (the_geom_lv03);"

Hat das geklappt, kann es ans Testen gehen. Die Syntax ist sowohl im README wie auch im vorangegangen Beitrag erläutert.

Ok, das ist jetzt noch keine Wissenschaft. Aber es scheint zu funktionieren. Die erste Frage muss aber sein: «Stimmt die Transformation überhaupt?» Um dies zu beantworten, vergleichen wir die Resultate aus zwei unabhängigen Transformationen. Einmal transformieren wir einen Testdatensatz mit der neuen ST_Fineltra-Funktion und einmal mit einem offiziellen Transformationsdienst der Swisstopo.

Diesen Vergleich können wir auch innerhalb der Datenbank machen. Dazu verwenden wir die PostgreSQL-Extension http. Mit dieser Extension wird PostgreSQL zu einem HTTP-Client. Bevor man das jetzt produktiv wirklich nutzen will, sollte man unbedingt das Kapitel «Why This is a Bad Idea» im README lesen.

Die Installation ist simpel:

git clone https://github.com/pramsey/pgsql-http.git
cd pgsql-http
make
sudo make install

Und anschliessend:

sudo -u postgres psql -d xanadu2 -c "CREATE EXTENSION http;"

Ein Testaufruf zeigt folgendes:

Mit den JSON-Funktionen von PostgreSQL kann man einfach auf die einzelnen Rückgabewerte zugreifen. Als Testdatensatz wählen wir die Fixpunkte aus der amtlichen Vermessung. Weil jetzt für jeden Fixpunkt ein HTTP-GET-Request gemacht wird, dauert das ein Weilchen:

WITH fineltra AS (
  SELECT nummer, nbident,
    ST_X(wkb_geometry) as x_lv03, ST_Y(wkb_geometry) as y_lv03,
    ST_X((ST_Fineltra(wkb_geometry, 'av_chenyx06.chenyx06_triangles', 'the_geom_lv03', 'the_geom_lv95'))) as x_lv95_pg,
    ST_Y((ST_Fineltra(wkb_geometry, 'av_chenyx06.chenyx06_triangles', 'the_geom_lv03', 'the_geom_lv95'))) as y_lv95_pg
  FROM av_avdpool_test.fixpunktekategorie__lfp
  --WHERE fid = 3515717
  LIMIT 100
),
swisstopo AS (
  SELECT nummer, nbident, x_lv95_pg, y_lv95_pg,
   cast(content::json->>'easting' as double precision) as x_lv95_lt,
   cast(content::json->>'northing' as double precision) as y_lv95_lt
  FROM fineltra as f, http_get('http://geodesy.geo.admin.ch/reframe/lv03tolv95?easting=' || f.x_lv03 || '&northing=' || f.y_lv03 || '&format=json')
)
SELECT nummer, nbident, (x_lv95_pg - x_lv95_lt)*1000 as x_diff_mm, (y_lv95_pg - y_lv95_lt)*1000 as y_diff_mm
FROM swisstopo;

Das Resultat sieht vielversprechend aus:

Die Differenzen bewegen sich also in einem völlig irrelevanten Bereich. Um sicher zu sein, müssen auch noch weitere Geometrietypen (Linien und Polygone) überprüft werden.

Zu guter Letzt soll die Performanz getestet werden. Dazu werden alle Tabellen der amtlichen Vermessung des Kantons Solothurn in einem MOpublic-ähnlichen Datenmodell in die Datenbank importiert. Eine GeoPackage-Datei des gesamten Kantons kann hier heruntergeladen und mit folgendem ogr2ogr-Befehl in die Datenbank importiert werden:

ogr2ogr -f "PostgreSQL" PG:"dbname='xanadu2' host='localhost' port='5432' user='stefan' password='ziegler12'" -lco SCHEMA=av_avdpool_test kanton.gpkg

In das Schema av_avdpool_test wurden 33 Tabellen importiert. Alle im Bezugrahmen LV03:

In einem Groovy-Skript reicht jetzt eigentlich eine einzige For-Schleife. Aus der View geometry_columns werden alle zu transformierenden Tabellen identifiziert und diese transformiert:

@Grapes([
   @Grab('org.postgresql:postgresql:9.4-1201-jdbc41'),
   @GrabConfig(systemClassLoader = true)
])

import groovy.sql.*

def dbhost = "localhost"
def dbport = "5432"
def dbdatabase = "xanadu2"
def dbusr = "stefan"
def dbpwd = "ziegler12"
def dbschema = "av_avdpool_test"

def dburl = "jdbc:postgresql://${dbhost}:${dbport}/${dbdatabase}?user=${dbusr}&password=${dbpwd}"

def query = "SELECT f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, type " +
  " FROM geometry_columns " +
  " WHERE f_table_schema = '$dbschema'" +
  " AND srid = 21781;"

def sql = Sql.newInstance(dburl)

def startTime = Calendar.instance.time
def endTime
println "Start: ${startTime}."

sql.withTransaction {
  sql.eachRow(query) {row ->
    def tableName = row.f_table_name
    def geomColumn = row.f_geometry_column

    def geomType = row.type
    if (row.coord_dimension == 3) {
      geomType += 'Z'
    }

    def alterQuery = "ALTER TABLE ${Sql.expand(dbschema)}.${Sql.expand(tableName)}" +
      " ALTER COLUMN ${Sql.expand(geomColumn)} TYPE geometry(${Sql.expand(geomType)},2056)" +
      " USING ST_Fineltra(${Sql.expand(geomColumn)}, 'av_chenyx06.chenyx06_triangles', 'the_geom_lv03', 'the_geom_lv95');"

    println "--- $dbschema.$tableName ---"
    sql.execute(alterQuery)

    endTime = Calendar.instance.time
    println "Elapsed time: ${(endTime.time - startTime.time)} ms"

  }
}

endTime = Calendar.instance.time
println "End: ${endTime}."
println "Total elapsed time: ${(endTime.time - startTime.time)} ms"

sql.connection.close()
sql.close()

Die Transformation einer Datenbanktabelle ist die Query in den Zeilen 38 - 40. Dies ist der einfachst mögliche Fall. Gibt es aber Rules und/oder Triggers etc. in der Tabelle, müssen diese wahrscheinlich vorher ausgeschaltet und nach der Transformation eingeschaltet werden. Das genaue Vorgehen muss geprüft werden.

Die gleichen 33 Tabellen liegen nun im Bezugsrahmen LV95 vor:

Die Geschwindigkeit ist verblüffend: Für sämtliche Tabellen braucht das Skript bloss circa 140 Sekunden. Die Transformation der Bodenbedeckung mit circa 280'000 Polygonen dauert nur 25 Sekunden (Ubuntu 14.04 in einer VirtualBox auf einem 2jährigen iMac mit SSD).

Posted by Stefan Ziegler. | PostGIS , Bezugsrahmenwechsel , LV95 , Fineltra