#!/bin/bash
#
# Kommandos zum generieren der Daten, als Beispielreferenz für GDAL tools mit ALS Daten
#

# 0. Quelldaten bereitstellen
mkdir tmp
unzip -o ../../srcdata/DOM/ALS_DGM/Oberösterreich-DGM_10M.zip -d tmp/


# 1. Generierung der Höhendaten (DEM) aus den ALS-Daten
# =======================================================
# Diese sind in EPSG:31255 mit Gitterweite 10m;  Zielsystem ist EPSG:31258 (BMN M31)
# Zielgröße ist die Boundingbox in EPSG:31258:
#  #xllcorner    469000
#  #yllcorner    263500
#  #xurcorner    476700
#  #yurcorner    268400
gdalwarp -t_srs EPSG:31258   -te 469500 263000 476700 268400  -tr 10 10     -dstnodata 9999 -r cubic       -overwrite  tmp/DGM_10M.tif tmp/dem.tif
#         ^tgt-Koordinatensystem ^llx   ^lly   ^urx   ^ury    ^Gitterweite  ^nodata-wert    ^Interpolation ^Zieldaten Überschreiben 

# Ergebnis prüfen (zeigt generierte Boundingbox):
# Uns interessieren die Eck-Koordinaten. Das sind nicht immer die angeforderten Werte, sondern manchmal auch leicht abweichende.
# Das liegt am Raster.
gdalinfo tmp/dem.tif |grep -A4 "Corner Coordinates" |grep "Lower Left"
gdalinfo tmp/dem.tif |grep -A4 "Corner Coordinates" |grep "Upper Right"
# Ergebnis ist (in dem Fall passts genau):
#   Lower Left  (  469500.000,  263000.000)
#   Upper Right (  476700.000,  268400.000)
# => Ergebnis ist also: tmp/dem.tif in EPSG:31258 mit 10m-Raster ist genau passend


# Jetzt die Theriondaten da draus bauen (Therion kann das noch nicht direkt verwenden)
# als erstes müssen wir es ins ASCII-Format bringen:
gdal_translate -co DECIMAL_PRECISION=5 -of AAIGrid tmp/dem.tif tmp/dem.txt

# Nun ergänzen wir die Therion-Koorindatensystem-headerdaten.
# und erzeugen dabei gleich das grid file draus, das im therionfile dann eingebunden wird.
# Die Daten für grid bekommen wir direkt aus der Quelldatei:
#   ncols        720
#   nrows        540
#   xllcorner    469500.000000000000
#   yllcorner    263000.000000000000
#   cellsize     10.000000000000
cat - tmp/dem.txt > hirlatz-dom-10mALS.grid <<THDATA
cs EPSG:31258
grid 469500 263000 10 10 720 540

THDATA

# Die informativen ersten header-Zeilen kommentiert man einfach aus.
sed -i '/^\(ncols\|nrows\|xllcorner\|yllcorner\|cellsize\|NODATA_value\)/s/^/#/' hirlatz-dom-10mALS.grid



# 2. Generierung der Bilddaten (DOM) aus den Ortophotos
# =======================================================
# Mit den genauen Eckkoordinaten aus dem Geländemodell (gdalinfo oben) können wir nun das
# Oberflächenbild aus einem Ortophoto (Bilddatei+Geoinformationen, z.b. XML) extrahieren.
# Das Ortophoto lässt sich z.b. mit QGis erzeugen; hier benutzen wir ein fertiges QGis-Projekt, das
# auf ArcGis mit Dachstein eingestellt wurde.
# Die auflösung wird bestimmt:
#  - einerseits durch die grundlegenden Daten des Service/Bildes (höher geht nicht; QGis-Zoomstufe beachten!)
#  - andererseits durch die Zielgröße des Bildes.
# Wir wollen hier eine 1m Auflösung haben, d.h. wir brauchen die Anzahl Pixel der Meterausdehnung:
#   Lower Left  (  469500.000,  263000.000)
#   Upper Right (  476700.000,  268400.000)
#  x=469500-476700=7200;  y=268400-263000=5400
qgis --project ../../tool/qgis/Dachstein_ArcGIS-2.qgs --snapshot tmp/surface.png --width 7200 --height 5400 --extent 469500,263000,476700,268400

# PNG geht in therion nicht, also Konvertieren in Geo-jpeg:
gdal_translate tmp/surface.png -of jpeg hirlatz-dop-1m_geoland.jpeg


# Ausgabe der Bilddimensionen, damit wir für das bitmap-Kommando die Daten bekommen:
gdalinfo hirlatz-dop-1m_geoland.jpeg |grep -A2 "Size is"
gdalinfo hirlatz-dop-1m_geoland.jpeg |grep -A4 "Corner Coordinates" |grep "Lower Left"
gdalinfo hirlatz-dop-1m_geoland.jpeg |grep -A4 "Corner Coordinates" |grep "Upper Right"
# Ergebnis:
# Size is 7200, 5400
# Origin = (469499.666543164115865,268400.000000000000000)
# Pixel Size = (1.000370507595406,-1.000370507595406)
# Lower Left  (  469499.667,  262997.999) 
# Upper Right (  476702.334,  268400.000)


# 3. Theriondatei bauen. Diese bindet alles zusammen
#==========================================================
cat > hirlatz-dop-1m_geoland-dom-10mALS.th <<THDATA
surface
  input hirlatz-dom-10mALS.grid

  cs EPSG:31258
  bitmap hirlatz-dop-1m_geoland.jpeg [0 0 469499.667 262997.999 7200 5400 476702.334 268400.000
endsurface
THDATA
    
    
    
# Ende, aufräumen
rm -rf tmp/
cd -