diff --git a/DSM/BE/BE_DOM_tiles.csv.gz b/DSM/BE/BE_DOM_tiles.csv.gz deleted file mode 100644 index c5ba0f4..0000000 Binary files a/DSM/BE/BE_DOM_tiles.csv.gz and /dev/null differ diff --git a/DSM/BE/BE_DSM_tindex.py b/DSM/BE/BE_DSM_tindex.py new file mode 100644 index 0000000..70477a7 --- /dev/null +++ b/DSM/BE/BE_DSM_tindex.py @@ -0,0 +1,176 @@ +############################################################################ +# +# NAME: BE_DSM_tindex.py +# +# AUTHOR(S): Anika Weinmann +# mundialis GmbH & Co. KG, Bonn +# https://www.mundialis.de +# +# PURPOSE: Create tile index of Berlin DOM/DSM txt files +# +# Data source: https://www.opengeodata.nrw.de/produkte/geobasis/dop/dop/ +# +# COPYRIGHT: (C) 2023 by Markus Neteler, Anika Weinmann, mundialis +# +# REQUIREMENTS: lynx, gdal, gzip, sed +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +############################################################################ +# Usage: +# Set Parameter `GDALTINDEX` +# `GDALTINDEX=True`: using `gdaltindex` (more correct limits of the tiles) +# `GDALTINDEX=False`: creating the tiles out of the file names (only small +# improvement in runtime, because you still have to check +# with `gdalinfo` whether the tiles are valid files at all) +# Then call script like this: +# python3 DSM/BE/BE_DSM_tindex.py +# Output: +# DSM/BE/be_dom_tindex_proj.gpkg.gz + + +import os +import json + +from osgeo import gdal +from remotezip import RemoteZip + + +# Parameter for Berlin DOM txt files +URL = ( + "https://fbinter.stadt-berlin.de/fb/berlin/service_intern.jsp?" + "id=a_dom1@senstadt&type=FEED" +) +GREP_STR = "https://fbinter.stadt-berlin.de/fb/atom/DOM/DOM1" +EPSG_CODE = 25833 +FILE_EXTENSION = ".txt" +TILE_SIZE = 2000 +OUTPUT_FILE = "be_dom_tindex_proj.gpkg.gz" +GDALTINDEX = True +VSI_PART = "/vsizip/vsicurl/" +os.chdir("DSM/BE/") + + +def create_tindex_by_filename(data_list): + # create GeoJson + geojson_dict = { + "type": "FeatureCollection", + "name": "tindex", + "crs": { + "type": "name", + "properties": {"name": f"urn:ogc:def:crs:EPSG::{EPSG_CODE}"}, + }, + "features": [] + } + + for num, data in enumerate(data_list): + ds = gdal.Open(data) + if ds is None: + print( + f"<{data}> is not a valid file and will be left out for " + "tindex creation." + ) + continue + splitted_data_name = os.path.basename( + data + ).replace(FILE_EXTENSION, "").split("_") + x1 = int(splitted_data_name[1]) * 1000 + y1 = int(splitted_data_name[2]) * 1000 + x2 = x1 + TILE_SIZE + y2 = y1 + TILE_SIZE + feat = { + "type": "Feature", "properties": { + "fid": num + 1, + "location": data, + }, + "geometry": { + "type": "Polygon", "coordinates": [[ + [x1, y1], + [x2, y1], + [x2, y2], + [x1, y2], + [x1, y1] + ]], + } + } + geojson_dict["features"].append(feat) + + with open("tindex.geojson", "w") as f: + json.dump(geojson_dict, f, indent=4) + + # create GPKG from GeoJson + tindex_gpkg = OUTPUT_FILE.rsplit(".", 1)[0] + stream = os.popen(f"ogr2ogr {tindex_gpkg} tindex.geojson") + stream.read() + return tindex_gpkg + + +def create_tindex_by_gdaltindex(data_list): + tiles_csv = "tiles.csv" + tindex_gpkg = OUTPUT_FILE.rsplit(".", 1)[0] + with open(tiles_csv, "w") as f: + for data in data_list: + f.write(f"{data}\n") + # create tile index + tindex_cmd = ( + f"gdaltindex -f GPKG -t_srs EPSG:{EPSG_CODE} {tindex_gpkg} " + f"--optfile {tiles_csv}" + ) + stream = os.popen(tindex_cmd) + stream.read() + os.remove(tiles_csv) + return tindex_gpkg + + +# check if we have lynx tool +stream = os.popen("which lynx") +output = stream.read() +if output is None or output == "": + raise Exception("lynx required, please install lynx first") + +# full tile index with 35860 NRW DOPs +get_data_cmd = ( + f"lynx -dump -nonumbers -listonly '{URL}' | grep {GREP_STR} | grep 'zip$' " + f"| sed 's+^+{VSI_PART}+g'" +) +stream = os.popen(get_data_cmd) +data_str = stream.read() +data_list = [] +for data in data_str.split(): + with RemoteZip(data.replace(VSI_PART, "")) as zip: + for zip_info in zip.infolist(): + file_name = zip_info.filename + data_list.append(f"{data}/{file_name}") + +# create tindex +if GDALTINDEX: + tindex_gpkg = create_tindex_by_gdaltindex(data_list) +else: + tindex_gpkg = create_tindex_by_filename(data_list) + +# verify +print("Verifying vector tile index:") +stream = os.popen(f"ogrinfo -so -al {tindex_gpkg}") +tindex_verification = stream.read() +print(tindex_verification) + +# package +if os.path.isfile(OUTPUT_FILE): + os.remove(OUTPUT_FILE) +stream = os.popen(f"gzip {tindex_gpkg}") +stream.read() +print(f"<{OUTPUT_FILE}> created") + +# cleanup +if os.path.isfile("tindex.geojson"): + os.remove("tindex.geojson") +if os.path.isfile(tindex_gpkg): + os.remove(tindex_gpkg) diff --git a/DSM/BE/be_dom_tindex_proj.gpkg.gz b/DSM/BE/be_dom_tindex_proj.gpkg.gz index 292a412..a0d28d4 100644 Binary files a/DSM/BE/be_dom_tindex_proj.gpkg.gz and b/DSM/BE/be_dom_tindex_proj.gpkg.gz differ diff --git a/DSM/BE/fetch_BE_DOM_list.sh b/DSM/BE/fetch_BE_DOM_list.sh deleted file mode 100644 index 1d1e7ea..0000000 --- a/DSM/BE/fetch_BE_DOM_list.sh +++ /dev/null @@ -1,57 +0,0 @@ -#!/bin/sh - -############################################################################ -# -# NAME: fetch_BE_DOM_list.sh -# -# AUTHOR(S): Markus Neteler -# Veronica Koess -# mundialis GmbH & Co. KG, Bonn -# https://www.mundialis.de -# -# PURPOSE: Fetch list of Berlin, Germany DOM (digital surface model) xyz .zip files -# -# Data source: https://fbinter.stadt-berlin.de/ -# -# COPYRIGHT: (C) 2023 by Markus Neteler, mundialis -# -# REQUIREMENTS: lynx, gdal, gTILE, sed -# -# This program is free software; you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation; either version 2 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -############################################################################ - -# Usage: -# sh fetch_BE_DOM_list.sh -# Output: -# BE_DOM_tiles.csv -######################################## - -# fail on error -set -e - -lynx -dump -nonumbers -listonly "https://fbinter.stadt-berlin.de/fb/berlin/service_intern.jsp?id=a_dom1@senstadt&type=FEED" | grep https://fbinter.stadt-berlin.de/fb/atom/DOM/DOM1 | sed 's+^+/vsizip/vsicurl/+g' > BE_DOM_tiles.csv - -# generate download script -# cat BE_DOM_tiles.csv | sed 's+^+wget -c +g' > fetch_be_dom_urls.sh - -# compress ndom50 URLs list -gzip BE_DOM_tiles.csv -echo "Generated " - -echo "" -echo "Single tile import: Import into GRASS GIS with, e.g.: -r.import input=/vsicurl/https://fbinter.stadt-berlin.de/fb/atom/DOM/DOM1_368_5808.zip output=be_dom1_368_5808" -echo "" -echo "For mosaics, better generate a VRT mosaic first (using ), then import the VRT file." -echo "" -echo "For a tile index, run -gdaltindex -f GPKG BE_DOM_tileindex.gpkg --optfile BE_DOM_tiles.csv" diff --git a/DSM/BE/requirements.txt b/DSM/BE/requirements.txt new file mode 100644 index 0000000..510d15b --- /dev/null +++ b/DSM/BE/requirements.txt @@ -0,0 +1 @@ +remotezip \ No newline at end of file diff --git a/DTM/BE/BE_DTM_tindex.py b/DTM/BE/BE_DTM_tindex.py index 1074dac..6bc0864 100644 --- a/DTM/BE/BE_DTM_tindex.py +++ b/DTM/BE/BE_DTM_tindex.py @@ -28,7 +28,9 @@ # Usage: # Set Parameter `GDALTINDEX` # `GDALTINDEX=True`: using `gdaltindex` (more correct limits of the tiles) -# `GDALTINDEX=False`: creating the tiles out of the file names (only small improvement in runtime, because you still have to check with `gdalinfo` whether the tiles are valid files at all ) +# `GDALTINDEX=False`: creating the tiles out of the file names (only small +# improvement in runtime, because you still have to check +# with `gdalinfo` whether the tiles are valid files at all) # Then call script like this: # python3 DTM/BE/BE_DTM_tindex.py # Output: @@ -42,12 +44,17 @@ # Parameter for Berlin DGM XYZ files -URL = "https://fbinter.stadt-berlin.de/fb/berlin/service_intern.jsp?id=a_dgm@senstadt&type=FEED" +URL = ( + "https://fbinter.stadt-berlin.de/fb/berlin/service_intern.jsp?" + "id=a_dgm@senstadt&type=FEED" +) GREP_STR = "https://fbinter.stadt-berlin.de/fb/atom/DGM1/DGM1_" EPSG_CODE = 25833 +FILE_EXTENSION = ".xyz" TILE_SIZE = 2000 -OUTPUT_FILE = "be_dgm_tindex_proj_gdal.gpkg.gz" +OUTPUT_FILE = "be_dgm_tindex_proj.gpkg.gz" GDALTINDEX = True +VSI_PART = "/vsizip/vsicurl/" os.chdir("DTM/BE/") @@ -73,7 +80,7 @@ def create_tindex_by_filename(data_list): continue splitted_data_name = os.path.basename( data - ).replace(".xyz", "").split("_") + ).replace(FILE_EXTENSION, "").split("_") x1 = int(splitted_data_name[1]) * 1000 y1 = int(splitted_data_name[2]) * 1000 x2 = x1 + TILE_SIZE @@ -129,11 +136,14 @@ def create_tindex_by_gdaltindex(data_list): raise Exception("lynx required, please install lynx first") # full tile index with 35860 NRW DOPs -get_data_cmd = f"lynx -dump -nonumbers -listonly '{URL}' | grep {GREP_STR} | grep 'zip$' | sed 's+^+/vsizip/vsicurl/+g'" +get_data_cmd = ( + f"lynx -dump -nonumbers -listonly '{URL}' | grep {GREP_STR} | grep 'zip$' " + f"| sed 's+^+{VSI_PART}+g'" +) stream = os.popen(get_data_cmd) data_str = stream.read() data_list = [ - f"{x}/{os.path.basename(x).replace('.zip', '.xyz')}" + f"{x}/{os.path.basename(x).replace('.zip', FILE_EXTENSION)}" for x in data_str.split() ]