[gelöst] Laserhöhen Österreich mit Versatz

Wie ich aus der Wochennotiz 244 gelesen habe, stehen für Österreich laservermessene Höhendaten (10m-Auflösung) unter CC-BY-Lizenz zur Verfügung.
Ich habe mal versucht diese Daten als Höhenrelief in meine eigenen Karte einzubauen. Ziel war es die hochauflösenden Laserdaten (www.data.gv.at) mit bestehenden SRTM1v3-Daten zu einem gemeinsamen Höhenrelief zusammenzubauen.

Mein Problem
Ich habe einen Versatz von ca. 100m, den ich mir nicht erklären kann. Hier die Lasterdaten und STRM-Daten im direkten Vergleich.

Normale SRTM-Daten

Laserdaten

Hat jemand eine Idee oder so was schon mal hinbekommen?

Grüße Andreas

Was ich gemacht habe…

Downloads…

  1. Download der Höhendaten (3 GB große ZIP-Datei mit 3 GB-GeoTIF-Inhalt dhm_lamb_10m.tif).
  2. Download anderen SRTM1v3-Höhendaten für den außerösterreichischen Bereich als *.hgt

GeoTif aus normalen STRM-Daten erstellen…
3) gdal_merge.py -v -o srtm_raw.tif -ul_lr 10 48 15 45 *.hgt
4) gdal_translate -of GTiff -co “TILED=YES” -a_srs “+proj=latlong” srtm_raw.tif srtm_adapted.tif
5) gdalwarp -of GTiff -co “TILED=YES” -srcnodata 32767 -t_srs “+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m” -rcs -order 3 -tr 30 30 -multi srtm_adapted.tif srtm_warped.tif

Laserdaten-GeoTif projektionsmäßig umwandeln…
6) gdal_translate -of GTiff -c “TILED=YES” -a_srs “+proj=latlong” dhm_lamb_10m.tif at_adapted.tif
7) gdalwarp -of GTiff -co “TILES=YES” -srcnodata -32768 -t_srs “+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m” -s_srs " +proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.33333333333333 +x_0=400000 +y_0=400000 +datum=hermannskogel +units=m +no_defs +ellps=bessel +towgs84=653.0,-212.0,449.0" -rcs -order 3 -r bilinear -multi at_adopted.tif at_warped.tif

Beide Höhendatenquellen zu einer Datei zusammenfassen…
8) gdalwarp -of GTiff -co “TILES=YES” -srcnodata -32768 -t_srs “+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m” -rcs -order 3 -multi -tr 10 10 -r bilinear srtm_warped.tif at_warped.tif ele.tif

Höhenrelief erstellen…
9) gdaldem hillshade ele.tif hillshade.tif -z 2

Kürzlich gelesen: http://www.mail-archive.com/talk-at@openstreetmap.org/msg07228.html

Andreas denk auch daran, dass die Positionen der Gipfel und Wege nicht unbedingt mit GPS-Geräten bestimmt wurden und von einem Luftbild abgemalt wurden.
Alleine die Positionsunterschiede aus den Lufbildern BING/Geoimage sind teilweise beträchtlich.

Elmauer Halt und Hinterer Goinger Halt sind mMn. durch Tracks passabel festgelegt.

@maxbe
Hi Max, danke für den Tip. Habe den Transformbefehl mal wie im Talk-Beitrag geändert. Leider ist der Versatz auch hier zu sehen

Durchgeführte Schritte (habe mal der Einfachheit halber den Merge mit den SRTM-Daten nicht durcheführt und nur die Laserdaten verwendet):

  1. gdal_translate -a_srs EPSG:31287 -co “TILED=YES” -co “BLOCKXSIZE=128” -co “BLOCKYSIZE=128” -co “COMPRESS=LZW” dhm_lamb_10m.tif at_adapted.tif
  2. gdalwarp -of GTiff -co “TILES=YES” -srcnodata -32768 -t_srs “+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m” -s_srs " +proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.33333333333333 +x_0=400000 +y_0=400000 +datum=hermannskogel +units=m +no_defs +ellps=bessel +towgs84=653.0,-212.0,449.0" -rcs -order 3 -r bilinear -multi at_adapted.tif at_warped.tif
  3. gdaldem hillshade at_warped.tif hillshade.tif -z 2

@noframe
Hi Horst, den Versatz glaube ich dadurch zu sehen, dass die STRM-Daten mit den GSS-Tracks&OSM-Daten zusammenpassen, während die Laserdaten um ca. 100m versetzt GPS-Tracks&OSM erscheinen. Für mich sind die meisten Projektionsparameter spanische Dörfer und denke, dass ich da noch irgendwas falsch umwandle :slight_smile:

Grüße
Andreas

Ich habs mal ohne den ganzen s_srs-Teil in deinem Punkt 2 versucht: Sieht genauso verschoben aus… Für “EPSG 31287” hab ich die Werte von spatialreference.org genommen.

Ich würde auf talk-at fragen, irgendwer dort hat sicher schon die Schummerung in der Schublade liegen und weiss den Weg dorthin…

Nachtrag: Schön sind sie ja schon, die 10-Meter-Daten… die Berge wirken irgendwie kantiger im Vergleich zu den sanft geglätteten 90-Mater-Daten

Ist denn für das dhm_lamb_10m.tif die korrekte Projektion eingestellt? Guck das mal mit GDALInfo in die Datei rein, was da als Koordinatensystem/Projektion definiert ist.

Vielleicht haben die auch irgendwo nen komischen Ellipsoiden verwendet. Im Zweifelsfall mal @Nakaner fragen, der kennt sich mit sowas qua Ausbildung aus.

Sieht ziemlich richtig aus:

Driver: GTiff/GeoTIFF
Files: dhm_lamb_10m.tif
Size is 58808, 30323
Coordinate System is:
PROJCS["PCS Name = MGI_Austria_Lambert",
    GEOGCS["GCS Name = GCS_MGI|Ellipsoid = Bessel_1841|Primem = Greenwich|",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel 1841",6377397.155,299.1528128000033,
                AUTHORITY["EPSG","7004"]],
            AUTHORITY["EPSG","6312"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",46],
    PARAMETER["standard_parallel_2",49],
    PARAMETER["latitude_of_origin",47.5],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",400000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (106549.267203768889885,576922.512073625810444)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  106549.267,  576922.512) (  9d19'7.63"E, 49d 1'24.32"N)
Lower Left  (  106549.267,  273692.512) (  9d31'19.50"E, 46d17'55.02"N)
Upper Right (  694629.267,  576922.512) ( 17d21'50.31"E, 49d 1'22.34"N)
Lower Right (  694629.267,  273692.512) ( 17d 9'35.51"E, 46d17'53.15"N)
Center      (  400589.267,  425307.512) ( 13d20'28.29"E, 47d43'39.80"N)

Hi Max,
danke fürs testen, dann weiß ich wenigstens, dass es nicht an meinem Rechner oder nur an mir liegt. Ich habe noch eine aktuellere GDAL-Version (1.11.1) versucht, aber das Ergebnis ist das selbe. Auf talk-at nachfragen muss ich mir noch überlegen. Ich möchte nicht noch einen Kommunikationskanal…
Wenn Du einen talk-Zugang hast: Meinst Du, Du könntest für mich einen Eintrag reinsetzen mit einem Link auf dieses Forum? Ich weiß, ich bin unverschämt und faul :slight_smile:

Genau aus diesem Grund wollte ich sie ja verwenden :slight_smile:

Grüße
Andreas

Bringt es irgendwas, da statt “hermannskogel” “bessel” (oder “rauenberg” / “potsdam”) einzutragen? Ist aber mehr “stochern im Nebel”.

Nachdem ich so ziemlich alle Kombinationen durchprobiert habe (ausser preussischen Datümern in österreichischen Projektionen :wink: )… Nicht schön, aber passender: http://geo.dianacht.de/bilder/10m-leuchsturm.png

Lass doch mal in “2) gdalwarp …” das ganze Zeug hinter t_srs und s_srs weg und trag dir stattdessen die Projektionen in die /usr/share/proj/epsg ein (oder wo gdal sonst bei dir die EPSG-Codes hernimmt):

<31287> +proj=lcc +lat_1=49 +lat_2=46 +lat_0=47.5 +lon_0=13.33333333333333 +x_0=400000 +y_0=400000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs
<900913> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs <>

und schreib dann “gdalwarp -of GTiff … -t_srs EPSG:900913 -s_srs EPSG:31287 -rcs -order 3 -r bilinear -multi …”

Grüße, Max

Zu den Projektionen habe ich auch noch einen kleinen Vorschlag:

Ich befürchte, dass der Dateiname “dhm_lamb_10m.tif” einen Hinweis auf eine Lambert-Projektion beinhaltet - eventuell muss man das noch irgendwo berücksichtigen (ohne jetzt gdal näher zu kennen…)

Nachtrag hierzu auch noch:
http://www.esri-austria.at/downloads/coords_at.html

Das wird mit “proj=lcc” ausgedrückt: “Lamberts conformal conic”

Die Daten will ich haben:


(oben die sanften Hügel aus den Aster-Daten, unten ein Gebirge aus den Daten von geoland.at)

Ein Projekt für verregnete Sommer…

@maxbe
Hi Max,

bist ein Schatz :slight_smile: Sieht gut aus. Ich mach mir jetzt erst mal ein Bier auf. Keine Ahnung, wie oft ich heute unterschiedlichste Projektionen ausprobiert habe und gewarped habe…
Wenn ich eigenen Screenshot und reproduzierbar Schritte habe, aktualisiere ich den Beitrag.

@gormo

Hat leider am Versatz nichts geändert. Ich habe “bessel” versucht.

Grüße Andreas

Da mich das Thema Höhendaten schon einige Zeit beschäftigt, habe ich diesen Faden zum Anlass genommen, mich ein bischen einzulesen.

Meine Erkenntnisse bisher:
gdlwarp dient dazu, georeferenzierte Rasterdaten von einem in ein anderes räumlichen Bezugssystem zu transformieren.
Das Ausgangsbezugsystem wird mit dem Kommandozeilenparameter -s_srs, das Zielsystem mit** -t_srs** anggegeben. Die Angaben sind u.a. im PROJ.4-Format möglich.

Das Bezugssystem für das österreichische Höhenmodell, also das Ausgangssystem, ist mit EPSG:31287 angegeben, die PROJ.4-Definition für -s_srs müsste demnach “+proj=lcc +lat_1=49 +lat_2=46 +lat_0=47.5 +lon_0=13.33333333333333 +x_0=400000 +y_0=400000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs” sein.
Das Bezugssystem für die SRTM-Daten ist nach meinen Recherchen EPSG 4326, die PROJ.4-Definition für -t_srs also “+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs”

Grüße
Franz

Zum Testen ist QGIS ganz nützlich. Hierzu zuerst

gdalinfo ogd-10m-at/dhm_lamb_10m.tif

um die bereits genannte Ausgabe zu erhalten, nämlich:

Am Anfang findet sich “PROJCS[“MGI_Austria_Lambert””. Nun kann man das Bild dhm_lamb_10m.tif mit QGIS öffnen, Rechtsklick auf den Layer und Properties wählen. Im öffnenden Popup kann unter “Coordinate reference system” das Bezugssystem gewählt werden. Wenn man im Filter “austria” eingibt, findet sich in der Liste schnell ein “MGI / Austria Lambert”. Das passt von der Benennung her zur Ausgabe von gdalinfo und QGIS zeigt nun auch den Parameter für GDAL an: +proj=lcc +lat_1=49 +lat_2=46 +lat_0=47.5 +lon_0=13.33333333333333 +x_0=400000 +y_0=400000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs

Um zu verifizieren, dass das so passt habe ich nun noch in QGIS einen WMS-Layer ergänzt und mit dhm_lamb_10m.tif verglichen. Die Grate und Gipfel überlagern sich genau.
Als WMS hatte ich http://gis.lebensministerium.at/wmsgw/?key=4d80de696cd562a63ce463a58a61488d&FORMAT=image/jpeg&VERSION=1.1.1&SERVICE=WMS&REQUEST=GetMap&Layers=Luftbild_MR,Luftbild_1m,Luftbild_8m,Satellitenbild_30m&SRS={proj}&WIDTH={width}&HEIGHT={height}&BBOX={bbox} verwendet, den ich aus der Luftbildliste von JOSM habe.

Der Betrag enthält zwar keine neue Information, aber vielleicht helfen ja die Schritte/dieser Weg jemandem bei der Bestimmung der Projektionsparameter für andere Daten.

Hat etwas gedauert, aber jetzt klappt es auch bei mir :slight_smile:

Mein letztes Problem war, dass ich zwar die zwei Zeilen von Max in die Projektionsdatei epsg eintragen konnte, aber dies nichts gebracht hat. Nur wenn ich “neue” Projektionsnummern vergeben habe, klappte es. Nach langer Fehlersuche behaupte ich mal, dass gdal die Einstellungen aus den Projektionsdateien nur dann verwendet, wenn in den Dateien pcs.csv und cubewerx_extra.wkt (die sich bei meine Windows-GDAL außerhalb des proj-Verzeichnisses im Verzeichnis gdal-data befinden) die Projektionen nicht schon definiert sind. Und in den beiden Dateien standen bei mir noch die alten Werte für 31287 und 900913.

Ich habe folgende Schritte durchgeführt, zum gewünschten Ergebnis geführt haben:

Downloads…

  1. Download der Höhendaten (3 GB große ZIP-Datei mit 3 GB-GeoTIF-Inhalt dhm_lamb_10m.tif).
  2. Download anderen SRTM1v3-Höhendaten für den außerösterreichischen Bereich als *.hgt

GeoTif aus normalen STRM-Daten erstellen…
3) gdal_merge.py -v -o srtm_raw.tif -ul_lr 10 48 15 45 *.hgt
4) gdal_translate -of GTiff -co “TILED=YES” -a_srs “+proj=latlong” srtm_raw.tif srtm_adapted.tif
5) gdalwarp -of GTiff -co “TILED=YES” -srcnodata 32767 -t_srs “+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m” -rcs -order 3 -tr 30 30 -multi srtm_adapted.tif srtm_warped.tif

Projektion der Laserdaten-GeoTif umwandeln…
6) Einfügung folgender zwei Zeilen in die Datei epgs im proj\Share-Verzeichnis von gdal:
<31287000> +proj=lcc +lat_1=49 +lat_2=46 +lat_0=47.5 +lon_0=13.33333333333333 +x_0=400000 +y_0=400000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs <>
<900913000> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs <>
7) gdalwarp -of GTiff -co “TILES=YES” -srcnodata -32768 -t_srs EPSG:900913000 -s_srs EPSG:31287000 -rcs -order 3 -r bilinear -multi dhm_lamb_10m.tif at_warped.tif

Beide Höhendatenquellen zu einer Datei zusammenfassen…
8) gdalwarp -of GTiff -co “TILES=YES” -srcnodata -32768 -t_srs “+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m” -rcs -order 3 -multi -tr 10 10 -r bilinear srtm_warped.tif at_warped.tif dem.tif

Höhenrelief erstellen…
9) gdaldem hillshade dem.tif hillshade.tif -z 2

Danke an alle Beteiligten und speziell an Max :slight_smile:

Grüße
Andreas
Attributierung geändert

Nachtrag:
Hat mir keine Ruhe gelassen, warum die Manipulation der Proj-Dateien notwendig war und habe eine Lösung gefunden, mit der es auch ohne geht (im Vergleich zu meinem allerersten Versuch scheinen die Parameter +nadgrids und +towgs84 den Ausschlag zum Erfolg zu geben).


gdalwarp -of GTiff -t_srs "+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m +nadgrids=@null" -s_srs "+proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.33333333333333 +x_0=400000 +y_0=400000 +datum=hermannskogel +units=m +no_defs +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232" ... -multi dhm_lamb_10m.tif at_warped.tif

und als Kurzfassung

gdalwarp -of GTiff -t_srs EPSG:3857 -s_srs EPSG:31287 ... dhm_lamb_10m.tif at_warped.tif

Gesamter Ablauf, wenns jemand “nachkochen” möchte:

Downloads…

  1. Download der Höhendaten (3 GB große ZIP-Datei mit 3 GB-GeoTIF-Inhalt dhm_lamb_10m.tif).
  2. Download anderen SRTM1v3-Höhendaten für den außerösterreichischen Bereich als *.hgt
  3. Download der Österreichgrenze als shp-Datei. Die shp-Datei wird in Schritt 7 verwendet und sorgt dafür, dass die zum Teil sichtbaren Übergänge zwischen SRTM- und Laserdaten genau auf der Landesgrenze liegen und somit nicht mehr so auffallen.

GeoTif aus normalen STRM-Daten erstellen…
4) gdal_merge.py -v -o srtm_raw.tif -ul_lr 10 48 15 45 *.hgt
5) gdal_translate -of GTiff -a_srs “+proj=latlong” srtm_raw.tif srtm_adapted.tif
6) gdalwarp -of GTiff -srcnodata 32767 -t_srs EPSG:3857 -rcs -order 3 -multi -r bilinear srtm_adapted.tif srtm_warped.tif

Projektion der Laserdaten-GeoTif umwandeln…
7) gdalwarp -of GTiff -dstnodata -32768 -t_srs EPSG:3857 -s_srs EPSG:31287 -rcs -order 3 -r bilinear -ot int16 -multi -crop_to_cutline -cutline Boundaries_20150407_081106.shp dhm_lamb_10m.tif at_warped.tif

Beide Höhendatenquellen zu einer Datei zusammenfassen…
8) gdalwarp -of GTiff -dstnodata 32767 -rcs -order 3 -multi srtm_warped.tif at_warped.tif dem.tif

Höhenrelief erstellen…
9) gdaldem hillshade dem.tif hillshade.tif -z 2

Getestet habe ich das Ganze mit gdal 1.11.1 (win64), gdal 1.9 (linux)

Grüße
Andreas

Danke für die ausführliche Beschreibung! Das ist klasse!

Wenn man deine Beschreibung gelesen hat, findet man danach sogar den Hinweis in den FAQ von proj.4: “Note the strategic addition of +nadgrids=@null to the spherical projection definition”.

Bei mir ist nur irgendwann mal die nächste FAQ hängengeblieben: Mann muss auch scheinbar selbstverständliches Zeug angeben, wenn man Ellipsoide transformieren will. Andernfalls macht proj.4 keine Transformation, sagt einem aber auch nicht, dass da was fehlt. Deshalb verwende ich gerne die Daten aus der epsg-Datei, wenn was nicht klappt. Da vergesse ich nichts.