You are not logged in.

#1 2015-04-01 15:10:02

Andreas Binder
Member
From: Bavaria
Registered: 2010-06-26
Posts: 540

[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
21455374qy.jpg

Laserdaten
21455432mp.jpg

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

Last edited by Andreas Binder (2015-04-07 21:30:35)

Offline

#2 2015-04-01 15:13:54

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

Offline

#3 2015-04-01 15:30:08

Noframe
Member
From: Seelze
Registered: 2009-11-08
Posts: 340

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

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.


Gruß
Noframe

Offline

#4 2015-04-01 18:29:11

Andreas Binder
Member
From: Bavaria
Registered: 2010-06-26
Posts: 540

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

@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

21457835di.jpg

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 :-)

Grüße
Andreas

Offline

#5 2015-04-01 22:36:53

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

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

Last edited by maxbe (2015-04-01 22:57:13)

Offline

#6 2015-04-02 11:05:55

gormo
Member
Registered: 2013-08-01
Posts: 2,096
Website

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

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.


OSM hat nicht das Ziel bis Ende des Monats einen vollständigen Datensatz der Welt zu enthalten.
(nach S.W.) - Aber weil die Welt vielfältig ist, weil sie auch im Detail interessant ist, mag ich genaue Karten (nach C.)

Offline

#7 2015-04-02 11:47:35

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

gormo wrote:

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.

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)

Offline

#8 2015-04-02 12:57:56

Andreas Binder
Member
From: Bavaria
Registered: 2010-06-26
Posts: 540

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

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 :-)

"Schön sind sie ja schon, die 10-Meter-Daten..."

Genau aus diesem Grund wollte ich sie ja verwenden :-)

Grüße
Andreas

Offline

#9 2015-04-02 13:04:38

gormo
Member
Registered: 2013-08-01
Posts: 2,096
Website

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

Andreas Binder wrote:

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

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

Last edited by gormo (2015-04-02 13:05:30)


OSM hat nicht das Ziel bis Ende des Monats einen vollständigen Datensatz der Welt zu enthalten.
(nach S.W.) - Aber weil die Welt vielfältig ist, weil sie auch im Detail interessant ist, mag ich genaue Karten (nach C.)

Offline

#10 2015-04-02 13:53:14

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

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

Offline

#11 2015-04-02 14:04:04

tunnelbauer
Member
Registered: 2012-01-13
Posts: 1,011
Website

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

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

Last edited by tunnelbauer (2015-04-02 14:05:28)


Grüße
Thomas

Offline

#12 2015-04-02 14:05:39

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

tunnelbauer wrote:

Hinweis auf eine Lambert-Projektion beinhaltet - eventuell muss man das noch irgendwo berücksichtigen

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

Offline

#13 2015-04-02 17:53:57

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

Die Daten will ich haben:

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

Ein Projekt für verregnete Sommer...

Offline

#14 2015-04-02 18:21:36

Andreas Binder
Member
From: Bavaria
Registered: 2010-06-26
Posts: 540

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

@maxbe
Hi Max,

... Nicht schön, aber passender: http://geo.dianacht.de/bilder/10m-leuchsturm.png...

bist ein Schatz :-) 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

Bringt es irgendwas, da statt "hermannskogel" "bessel" (oder "rauenberg" / "potsdam") einzutragen?

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

Grüße Andreas

Offline

#15 2015-04-03 09:35:29

kartler175
Member
Registered: 2012-09-10
Posts: 272

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

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

Offline

#16 2015-04-03 15:01:06

Augustus Kling
Member
Registered: 2009-05-11
Posts: 119

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

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:

Driver: GTiff/GeoTIFF
Files: ogd-10m-at/dhm_lamb_10m.tif                                                                                                                                                                                                                                             
       ogd-10m-at/dhm_lamb_10m.tfw                                                                                                                                                                                                                                             
Size is 58808, 30323                                                                                                                                                                                                                                                           
Coordinate System is:                                                                                                                                                                                                                                                         
PROJCS["MGI_Austria_Lambert",                                                                                                                                                                                                                                                 
    GEOGCS["GCS_MGI",                                                                                                                                                                                                                                                         
        DATUM["Militar_Geographische_Institute",                                                                                                                                                                                                                               
            SPHEROID["Bessel_1841",6377397.155,299.1528128,                                                                                                                                                                                                                   
                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)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  NoData Value=-3.40282306073709653e+38

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/? … BOX={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.

Offline

#17 2015-04-03 18:59:04

Andreas Binder
Member
From: Bavaria
Registered: 2010-06-26
Posts: 540

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

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

21551306mn.jpg

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 :-)

Grüße
Andreas
Attributierung geändert

Last edited by Andreas Binder (2015-04-10 18:56:14)

Offline

#18 2015-04-07 21:27:55

Andreas Binder
Member
From: Bavaria
Registered: 2010-06-26
Posts: 540

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

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

Last edited by Andreas Binder (2015-04-07 21:31:34)

Offline

#19 2015-04-08 08:16:15

gormo
Member
Registered: 2013-08-01
Posts: 2,096
Website

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

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


OSM hat nicht das Ziel bis Ende des Monats einen vollständigen Datensatz der Welt zu enthalten.
(nach S.W.) - Aber weil die Welt vielfältig ist, weil sie auch im Detail interessant ist, mag ich genaue Karten (nach C.)

Offline

#20 2015-04-08 08:43:42

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

Andreas Binder wrote:

im Vergleich zu meinem allerersten Versuch scheinen die Parameter +nadgrids und +towgs84 den Ausschlag zum Erfolg zu geben

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.

Offline

#21 2015-04-10 19:05:41

Andreas Binder
Member
From: Bavaria
Registered: 2010-06-26
Posts: 540

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

Ich habe bezüglich der korrekten Lizenzangabe/Attributierung beim Datenanbieter (österreichisches Bundeskanzleramt) nachgefragt, da ich auf der Webseite nicht wirklich schlau wurde.

Aus der Antwort:

...
Eine mögliche Lizenznennung wäre: "Datenquelle: geoland.at - data.gv.at CC by sa AT 3.0"
Sollten Sie oder Ihre Kollegen Applikationen oder Visualisierungen basierend auf Open Data einreichen, können Sie diese auch gern zurück melden: https://www.data.gv.at/anwendung-einreichen/
...

Offline

#22 2015-07-27 14:11:34

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

Ich hab die Daten jetzt mal bei mir eingebaut. Die Verbesserung ist schon deutlich: Die Berge wirken irgendwie lebhafter. Warum der östliche Ausläufer der Ellmauer Halt "Kopftörlgrat" heisst, wirkt auch klarer und die Scharte zwischen der Hinteren Karlspitze und dem Totenkirchl liegt nicht mehr so unsinnig auf dem flachen Hang rum:

kaiserastergeoland.png

Besonders schön ist, dass die ganzen Fehler durch Fehlsichtigkeit und Blendung des Satelliten weg sind. Der Roßkopf ist leider doch kein Vulkan. Man sieht da auch, wie gut das Zeug von viewfinderpanoramas.org in der OpenTopoMap ist. Auch da fehlt der Krater zwischen den Gipfeln. Bei den SRTM-Daten von CGIAR-CSI, wie sie die Wanderreitkarte verwendet und den Höhendaten der OpenCycleMap (welchen?) ist der Berg interessanter, aber meines Wissens leider falsch:

4rosskopfvulkane.png

Grüße und Danke, Max (der sich den Ausflug zum Schneeloch sparen wird, das ist gar kein 3000m tiefer Schlot, eigentlich nicht mal ein Loch)

Offline

#23 2015-07-27 14:22:20

TobWen
Member
From: Ruhrgebiet
Registered: 2009-03-31
Posts: 1,108

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

Kurz zwei Sachen zur Schummerung:
1. Mittlerweile arbeiten die Kartographen häufig nicht mehr mit einer mit der einseitigen Belichtung aus Nordwesten, sondern verrechnen Informationen aus verschiedenen Richtungen.
2. ASTER v2 hat inzwischen eine 30-meter Auflösung, allerdings mit Artefakten.


Was macht der RVR mit OpenStreetMap? https://forum.openstreetmap.org/viewtopic.php?id=63052
Aktuelle Luftbilder des RVRs (Ruhrgebiet) http://forum.openstreetmap.org/viewtopic.php?id=28511

Offline

#24 2015-07-27 14:47:01

maxbe
Member
Registered: 2010-01-19
Posts: 3,150
Website

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

TobWen wrote:

1. Mittlerweile arbeiten die Kartographen häufig nicht mehr mit einer mit der einseitigen Belichtung aus Nordwesten, sondern verrechnen Informationen aus verschiedenen Richtungen.

Ich hätte ja so gerne Schraffen. Aber ich schau mir auch anderes gern mal an wink

TobWen wrote:

2. ASTER v2 hat inzwischen eine 30-meter Auflösung, allerdings mit Artefakten.

Mir ist die "Richtigkeit" eigentlich wichtiger als die Auflösung, vor allem die Vermeidung von Phantasiebergen und -Tälern. Tatsächlich ist das Bild da oben auch nur mit ca. 25m Auflösung gerechnet, weil ich habe jetzt diese 10m-Werte aus Österreich, 50m aus Bayern und den Rest in 90m ASTER und musste irgendeinen Kompromiss finden. Sobald ich die 2.5m-Daten aus Südtirol in einem handhabbaren Format habe, denke ich noch mal drüber nach.

Grüße, Max

Offline

Board footer

Powered by FluxBB