download_files_from_LGL_BW.py

"""
LoD2 CityGML tiles are available for whole Baden-Württemberg, from LGL.

https://opengeodata.lgl-bw.de/#/(sidenav:product/12)

This script downloads the requires tiles for given regions
(as WKT strings, in *polygons* variable), and extracts the region.

Required:
* Python
* pyproj project (https://pypi.org/project/pyproj/)
* SimStadt installed on the Desktop (for RegionChooser)

Eric Duminil, 2024
"""
from pathlib import Path
from math import floor
import subprocess
import re
# pip install pyproj
from pyproj import CRS
from pyproj import Transformer
import urllib.request
import time
import zipfile

###### User input ##########
POLYGONS = {
    "StuttgartCenter": "POLYGON((9.175287 48.780916, 9.185501 48.777522, 9.181467 48.773704, 9.174429 48.768472, 9.168807 48.773902, 9.175287 48.780916))",
    # "AnotherRegion": "Another WKT Polygon...",
}
# Should RegionChooser extract the regions from multiple CityGMLs?
EXTRACT_REGIONS = True
############################


CITYGML_SERVER = "https://opengeodata.lgl-bw.de/data/lod2"
RASTER = 2  # [km]
BUNDESLAND = 'bw'
WGS84 = CRS.from_epsg(4326)

LOCAL_CRS = CRS.from_epsg(32632)  # UTM32N, used in BW. https://epsg.io/32632
UTM = 32

SCRIPT_DIR = Path(__file__).parent
WAIT_BETWEEN_DOWNLOADS = 5  # [s] Be nice to LGL Server.

if EXTRACT_REGIONS:
    SIMSTADT_FOLDER = next(Path.home().glob('Desktop/SimStadt*_0.*/'))
    GML_GLOB = "LoD2_*/LoD2_*.gml"


def coordinates_to_grid(longitude: float, latitude: float, local_utm: CRS) -> tuple[int, int]:
    """Returns (x, y) of the tile on CITYGML_SERVER containing a given point."""
    transformer = Transformer.from_crs(WGS84, local_utm, always_xy=True)
    x, y = transformer.transform(longitude, latitude)
    x = floor(x / 1000) - 1  # Odd x
    y = floor(y / 1000)  # Even y
    x -= x % RASTER
    y -= y % RASTER
    return (x + 1, y)


def wkt_polygon_to_grid_coords(location_name: str, wkt: str) -> tuple[int, int, int, int]:
    """Returns (x, y) of lower-left and bottom-right tiles, containing a given region."""
    if 'POLYGON' not in wkt:
        raise ValueError(f"wkt for {location_name} should be a WKT POLYGON or MULTIPOLYGON")

    coordinates = re.findall(r'\-?\d+\.\d+', wkt)

    lons = [float(lon) for lon in coordinates[::2]]
    lats = [float(lat) for lat in coordinates[1::2]]

    min_lon, max_lon = min(lons), max(lons)
    min_lat, max_lat = min(lats), max(lats)

    print("%s (%.3f°N %.3f°E -> %.3f°N %.3f°E)" % (location_name, max_lat, min_lon, min_lat, max_lon))

    x1, y1 = coordinates_to_grid(min_lon, min_lat, LOCAL_CRS)
    x2, y2 = coordinates_to_grid(max_lon, max_lat, LOCAL_CRS)

    return (x1, x2, y1, y2)


def download_rectangle(output_dir: Path, x1: int, x2: int, y1: int, y2: int) -> None:
    """Downloads every zip of a given region, to output_dir, and extracts CityGML files."""
    for x in range(x1, x2 + 1, RASTER):
        for y in range(y1, y2 + 1, RASTER):
            citygml_zip = f"LoD2_{UTM}_{x}_{y}_{RASTER}_{BUNDESLAND}.zip"
            citygml_url = f"{CITYGML_SERVER}/{citygml_zip}"
            local_zip = output_dir / citygml_zip  # .replace('.zip', '.gml')
            if local_zip.exists():
                print(f"  {local_zip.name} already in {output_dir.name}/")
            else:
                print(f"  Download {citygml_zip} to {output_dir.name}/ ", end='')
                urllib.request.urlretrieve(citygml_url, local_zip)
                time.sleep(WAIT_BETWEEN_DOWNLOADS)
                print("✓")

                print(f"    Extract {citygml_zip} to {output_dir.name}/ ", end='')
                print("✓")
                print("")
                with zipfile.ZipFile(local_zip, "r") as zip_ref:
                    zip_ref.extractall(output_dir)


def extract_region(output_dir: Path, location_name: str, wkt: str) -> None:
    """Uses RegionChooser to extract a given region from all the CityGML files found in subfolder."""
    output_file = output_dir / (location_name + '.gml')
    if output_file.exists():
        print(f"  {output_file} already exists. Not extracting.")
        return
    print(f"  Extracting {output_file}.")
    region_chooser_libs = Path(SIMSTADT_FOLDER).expanduser() / 'lib/*'
    gml_inputs = output_dir.glob(GML_GLOB)
    result = subprocess.run(' '.join(['java', '-classpath', f'"{region_chooser_libs}"',
                                      'eu.simstadt.regionchooser.RegionChooserCLI',
                                      '--input', ','.join(str(gml) for gml in gml_inputs),
                                      '--output', str(output_file),
                                      '--wkt', '-',
                                      ]),
                            input=wkt,
                            text=True,
                            capture_output=True
                            )
    if (result.stderr):
        print(result.stderr)
    print("  DONE!")


def main(polygons: dict[str, str]) -> None:
    """Downloads ZIP files, extracts CityGML files, and selects desired region."""
    for location_name, wkt in polygons.items():
        output_dir = SCRIPT_DIR / (location_name + '.proj')
        output_dir.mkdir(parents=True, exist_ok=True)
        x1, x2, y1, y2 = wkt_polygon_to_grid_coords(location_name, wkt)
        download_rectangle(output_dir, x1, x2, y1, y2)
        if EXTRACT_REGIONS:
            extract_region(output_dir, location_name, wkt)
        print()


if __name__ == '__main__':
    main(POLYGONS)