Here :)
Am 2020-03-07 14:27, schrieb Rodrigo Severo via Therion:
‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
On Friday, March 6, 2020 9:34 AM, Benedikt Hallinger
<b...@hallinger.org> wrote:
... if it is needed, i can translate it to English. Just say so and i
try to squeeze out some time...
It would certainly help a lot.
Regards,
Rodrigo
_______________________________________________
Therion mailing list
Therion@speleo.sk
https://mailman.speleo.sk/listinfo/therion
#!/bin/bash
#
# Commands to generate therion DEM data with GDAL tools out of austrian ALS data
#
# 0. Make source data available in temporary folder
mkdir tmp
unzip -o ../../srcdata/DOM/ALS_DGM/Oberösterreich-DGM_10M.zip -d tmp/
# 1. Generate heightmap (DEM) from ALS-data
# =======================================================
# The source data is in EPSG:31255 with grid width 10m; our target is EPSG:31258 (BMN M31)
# Target size is 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-coord-system ^llx ^lly ^urx ^ury ^grid-width ^nodata-value ^interpolation ^overwrite tgt data if existing
# Verify results (shows generated boundingbox):
# The results are not always exactly the requested box coordinates which is a result of the underlaying grid.
gdalinfo tmp/dem.tif |grep -A4 "Corner Coordinates" |grep "Lower Left"
gdalinfo tmp/dem.tif |grep -A4 "Corner Coordinates" |grep "Upper Right"
# Results in this case (which fits exactly):
# Lower Left ( 469500.000, 263000.000)
# Upper Right ( 476700.000, 268400.000)
# => Result is DEM file in GEO-TIFF: tmp/dem.tif in EPSG:31258 with 10m-grid
# Lets build therion data with that (therion can't use the TIFF file directly,
# we need to convert into ASCII format):
gdal_translate -co DECIMAL_PRECISION=5 -of AAIGrid tmp/dem.tif tmp/dem.txt
# Now add therion header data.
# We will produce the grid file which can be used by therion.
# (The needed values for the "grid" command comes from inside the ASCII-grid file)
# 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
# We need to comment the informative gdal ascii header data.
# (could also be removed, but that information is nice to preserve)
sed -i '/^\(ncols\|nrows\|xllcorner\|yllcorner\|cellsize\|NODATA_value\)/s/^/#/' hirlatz-dom-10mALS.grid
# 2. Generate georeferenced picture (DOM) using QGIS
# =======================================================
# Using the exact boundingbox from the DEM (see gdalinfo above) we can now produce some
# surface image as overlay for the DEM model.
# For that we need some georeferenced image (e.g. ortophoto) covering the area.
# We can get this using an existing QGis project.
# This example assumes some project using the ArcGis WMS web service.
# You need to setup QGis projects coordinate system and also zoom level to get the picture you desire.
# Directly set QGis to use the desired target coordinate system so you won't need to convert it later on.
# If this is not possible, you can use gdalwarp to transform it.
#
# The resolution of the resulting piture is depending on:
# - the resoliution of the underlaying service (wms zoom level)
# - the target picture size in pixel
#
# In the example we want to achieve a 1m/pixel resolution which can be provided by the web service.
# So we need to calculate the needed picture size from the covered DEM area.
# EPSG:31258 is Gauss-Krueger (measurement in meters from grid origin) so this is easy:
# we just get the x and y dimensions in meters and use that as pixel target.
# 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
# ^stored QGis project ^generate picture ^pixel size ^selected bounding box in EPSG:31258
# The resiulting geo-png is not usable from therion, so lets convet to geo-jpeg:
# (this produces two files: a standard-jpeg picture and an additional XML file covering the GEO-information)
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:
# Use gdalifo to print information on the resulting picture (this also may not fit exactly the requested coordinates
# and we need the exact coordinates for the therion "bitmap"-command):
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"
# Result is (does not fit exactly, but 99,5%):
# 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. Generate a therion file, binding the DEM and DOM together
# (this is optional, yo can use the files directly. I howver found it more convinient
# to have a given file binding the results together. Also i can update the grid file independently)
#==========================================================
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
# done, clean up
rm -rf tmp/
cd -
_______________________________________________
Therion mailing list
Therion@speleo.sk
https://mailman.speleo.sk/listinfo/therion
