Commit 8ebe58c8 authored by Eric Duminil's avatar Eric Duminil
Browse files

Define Transformer only once

parent aad32b55
Showing with 8 additions and 8 deletions
+8 -8
...@@ -41,9 +41,11 @@ EXTRACT_REGIONS = True ...@@ -41,9 +41,11 @@ EXTRACT_REGIONS = True
CITYGML_SERVER = "https://opengeodata.lgl-bw.de/data/lod2" CITYGML_SERVER = "https://opengeodata.lgl-bw.de/data/lod2"
RASTER = 2 # [km] RASTER = 2 # [km]
BUNDESLAND = 'bw' BUNDESLAND = 'bw'
WGS84 = CRS.from_epsg(4326)
WGS84 = CRS.from_epsg(4326)
LOCAL_CRS = CRS.from_epsg(32632) # UTM32N, used in BW. https://epsg.io/32632 LOCAL_CRS = CRS.from_epsg(32632) # UTM32N, used in BW. https://epsg.io/32632
TO_LOCAL_CRS = Transformer.from_crs(WGS84, LOCAL_CRS, always_xy=True)
UTM = 32 UTM = 32
SCRIPT_DIR = Path(__file__).parent SCRIPT_DIR = Path(__file__).parent
...@@ -62,10 +64,9 @@ if EXTRACT_REGIONS: ...@@ -62,10 +64,9 @@ if EXTRACT_REGIONS:
GML_GLOB = "LoD2_*/LoD2_*.gml" GML_GLOB = "LoD2_*/LoD2_*.gml"
def coordinates_to_grid(longitude: float, latitude: float, local_utm: CRS) -> tuple[int, int]: def coordinates_to_grid(longitude: float, latitude: float) -> tuple[int, int]:
"""Returns (x, y) of the tile on CITYGML_SERVER containing a given point.""" """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 = TO_LOCAL_CRS.transform(longitude, latitude)
x, y = transformer.transform(longitude, latitude)
x = floor(x / 1000) - 1 # Odd x x = floor(x / 1000) - 1 # Odd x
y = floor(y / 1000) # Even y y = floor(y / 1000) # Even y
x -= x % RASTER x -= x % RASTER
...@@ -88,8 +89,8 @@ def wkt_polygon_to_grid_coords(location_name: str, wkt: str) -> tuple[int, int, ...@@ -88,8 +89,8 @@ def wkt_polygon_to_grid_coords(location_name: str, wkt: str) -> tuple[int, int,
print("%s (%.3f°N %.3f°E -> %.3f°N %.3f°E)" % (location_name, max_lat, min_lon, min_lat, max_lon)) 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) x1, y1 = coordinates_to_grid(min_lon, min_lat)
x2, y2 = coordinates_to_grid(max_lon, max_lat, LOCAL_CRS) x2, y2 = coordinates_to_grid(max_lon, max_lat)
return (x1, x2, y1, y2) return (x1, x2, y1, y2)
...@@ -181,8 +182,7 @@ def main(polygons: dict[str, str]) -> None: ...@@ -181,8 +182,7 @@ def main(polygons: dict[str, str]) -> None:
def show_coordinates(match): def show_coordinates(match):
longitude, latitude = match.groups() longitude, latitude = match.groups()
transformer = Transformer.from_crs(WGS84, LOCAL_CRS, always_xy=True) x, y = TO_LOCAL_CRS.transform(longitude, latitude)
x, y = transformer.transform(longitude, latitude)
return f"{x} {y}" return f"{x} {y}"
def convert_wkt_to_local(wkt): def convert_wkt_to_local(wkt):
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment