From 887e9e09dfb3115606b1e1b049b16acd1549140f Mon Sep 17 00:00:00 2001
From: Laurens Oostwegel <laurens@gfz-potsdam.de>
Date: Tue, 4 Apr 2023 11:31:21 +0200
Subject: [PATCH] Add arcsec changes to dev branch
---
borderline/arcsec.py | 11 ++++++--
borderline/exposure2boundary.py | 50 ++++++++++++++++++---------------
2 files changed, 37 insertions(+), 24 deletions(-)
diff --git a/borderline/arcsec.py b/borderline/arcsec.py
index 6726c7b..a43d270 100644
--- a/borderline/arcsec.py
+++ b/borderline/arcsec.py
@@ -61,6 +61,7 @@ class ArcsecBoundaryParser:
lon=lon_origin - 0.5 * steps, lat=lat_origin - 0.5 * steps
)
self.tiles = set()
+ self._tiles = {}
def create_grid(self, lons: List[float], lats: List[float]):
"""
@@ -159,7 +160,7 @@ class ArcsecBoundaryParser:
lat = self.origin_coordinate.lat + (tile_y * (self.arcsecs / 3600))
return lon, lat
- def get_grid_as_geodataframe(self):
+ def get_grid_as_geodataframe(self, id_generator):
"""
Retrieve all tiles and convert these to a GeoDataFrame with the tile indices and
the bounding box geometry as attributes.
@@ -168,13 +169,19 @@ class ArcsecBoundaryParser:
A GeoDataFrame of all tiles that are in the grid.
"""
- boundaries = {"geometry": [], "TILE_X": [], "TILE_Y": []}
+ boundaries = {"geometry": [], "TILE_X": [], "TILE_Y": [], "BOUNDARY_ID": []}
if not self.tiles:
logger.warning("Grid not filled yet. First use fill_grid.")
for tile in self.tiles:
+ if tile not in self._tiles:
+ new_id = id_generator()
+ self._tiles[tile] = new_id
+
boundaries["TILE_X"].append(tile[0])
boundaries["TILE_Y"].append(tile[1])
boundaries["geometry"].append(self.get_tile_geometry(*tile))
+ boundaries["BOUNDARY_ID"].append(self._tiles[tile])
+
boundaries_gdf = geopandas.GeoDataFrame(boundaries, geometry="geometry", crs=4326)
return boundaries_gdf
diff --git a/borderline/exposure2boundary.py b/borderline/exposure2boundary.py
index 87e085e..91bec84 100644
--- a/borderline/exposure2boundary.py
+++ b/borderline/exposure2boundary.py
@@ -103,8 +103,9 @@ class Exposure2Boundary:
# Get country as geopandas
country_gpd = self.country2boundary()
complete_boundary_df = country_gpd
+ arcsec_boundary = {}
- for exposure_type in ["residential", "commercial", "industrial"]:
+ for exposure_type in ["industrial", "commercial", "residential"]:
# Skip type if it is not there
if exposure_type not in self.config["exposure"]:
continue
@@ -128,22 +129,7 @@ class Exposure2Boundary:
arcsec_grid=arcsec_grid,
)
- clipped_boundaries = boundaries.clip(country_gpd)
-
- # Get all ArcSec cells that are outside of the clipped area and set the geometry
- # to `POINT EMPTY`
- boundaries_outside_country = boundaries.iloc[
- boundaries.index.difference(clipped_boundaries.index)
- ]
- boundaries_outside_country["geometry"] = wkt.loads("POINT EMPTY")
- boundaries_outside_country.set_geometry("geometry", crs=boundaries.crs)
-
- # merge all boundaries
- complete_boundary_df = gpd.GeoDataFrame(
- pd.concat(
- [complete_boundary_df, clipped_boundaries, boundaries_outside_country]
- )
- )
+ arcsec_boundary[arcsec_grid] = boundaries
exposure.drop("geometry", axis=1, inplace=True)
elif exposure_settings["method"] == "district-boundary":
@@ -221,6 +207,24 @@ class Exposure2Boundary:
exposure.to_csv(exposure_settings["output-file"], index=False)
+ # del self.arcsec_boundary_parsers
+ for arcsec_grid, boundaries in arcsec_boundary.items():
+ # Get all ArcSec cells that are outside of the clipped area and set the geometry
+ # to `POINT EMPTY`
+ clipped_boundaries = boundaries.clip(country_gpd)
+ boundaries_outside_country = boundaries.iloc[
+ boundaries.index.difference(clipped_boundaries.index)
+ ]
+ boundaries_outside_country["geometry"] = wkt.loads("POINT EMPTY")
+ boundaries_outside_country.set_geometry("geometry", crs=boundaries.crs)
+
+ # merge all boundaries
+ complete_boundary_df = gpd.GeoDataFrame(
+ pd.concat(
+ [complete_boundary_df, clipped_boundaries, boundaries_outside_country]
+ )
+ )
+ logger.info(f"{complete_boundary_df.shape[0]} boundaries found")
return complete_boundary_df
def country2boundary(self):
@@ -280,13 +284,15 @@ class Exposure2Boundary:
exposure[["TILE_X", "TILE_Y"]] = tiles
# Get boundaries of grid as GeoDataFrame
- boundaries = self.arcsec_boundary_parsers[arcsec_grid].get_grid_as_geodataframe()
-
- # Create boundary IDs
- boundaries["BOUNDARY_ID"] = boundaries.apply(
- lambda feat: self.get_new_boundary_id(), axis=1
+ boundaries = self.arcsec_boundary_parsers[arcsec_grid].get_grid_as_geodataframe(
+ self.get_new_boundary_id
)
+ # # Create boundary IDs
+ # boundaries["BOUNDARY_ID"] = boundaries.apply(
+ # lambda feat: self.get_new_boundary_id(), axis=1
+ # )
+
# Merge exposure with boundary so exposure also has the boundary IDs
exposure = pd.merge(
exposure, boundaries, how="inner", on=["TILE_X", "TILE_Y"], suffixes=(None, "_2")
--
GitLab