Skip to content

Commit fbd2e9b

Browse files
committed
Update 15 arcsec grid definition and formatting
1 parent ac00534 commit fbd2e9b

File tree

1 file changed

+42
-47
lines changed

1 file changed

+42
-47
lines changed

climada/entity/exposures/litpop/litpop.py

Lines changed: 42 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,6 @@
2727
import rasterio
2828
import shapely
2929
from shapefile import Shape
30-
import copy
3130
from affine import Affine
3231

3332
import climada.util.coordinates as u_coord
@@ -178,18 +177,14 @@ def from_countries(
178177
total_values = [None] * len(countries)
179178
elif len(total_values) != len(countries):
180179
raise ValueError(
181-
"'countries' and 'total_values' must be lists of same length"
182-
)
183-
180+
"'countries' and 'total_values' must be lists of same length")
184181
if target_grid is None:
185-
LOGGER.info(f"Creating target grid at {res_arcsec} arcsec resolution...")
186-
target_grid = cls._define_target_grid(
187-
reference_year=reference_year,
188-
gpw_version=gpw_version,
189-
data_dir=data_dir,
190-
res_arcsec=res_arcsec
191-
)
192-
182+
target_grid = cls._define_target_grid(
183+
reference_year=reference_year,
184+
gpw_version=gpw_version,
185+
data_dir=data_dir,
186+
res_arcsec=res_arcsec
187+
)
193188
# litpop_list is initiated, a list containing one Exposure instance per
194189
# country and None for countries that could not be identified:
195190
if admin1_calc: # each admin 1 region is initiated separately,
@@ -342,6 +337,13 @@ def from_nightlight_intensity(
342337
raise ValueError(
343338
"Not allowed to set both `countries` and `shape`. Aborting."
344339
)
340+
if target_grid is None:
341+
target_grid = cls._define_target_grid(
342+
reference_year=reference_year,
343+
gpw_version=GPW_VERSION,
344+
data_dir=data_dir,
345+
res_arcsec=res_arcsec
346+
)
345347
if countries is not None:
346348
exp = cls.from_countries(
347349
countries,
@@ -431,6 +433,13 @@ def from_population(
431433
raise ValueError(
432434
"Not allowed to set both `countries` and `shape`. Aborting."
433435
)
436+
if target_grid is None:
437+
target_grid = cls._define_target_grid(
438+
reference_year=reference_year,
439+
gpw_version=gpw_version,
440+
data_dir=data_dir,
441+
res_arcsec=res_arcsec
442+
)
434443
if countries is not None:
435444
exp = cls.from_countries(
436445
countries,
@@ -549,7 +558,6 @@ def from_shape_and_countries(
549558
shape_gdf = _get_litpop_single_polygon(
550559
shape,
551560
reference_year,
552-
res_arcsec,
553561
data_dir,
554562
gpw_version,
555563
exponents,
@@ -694,7 +702,6 @@ def from_shape(
694702
litpop_gdf, _ = _get_litpop_single_polygon(
695703
shape,
696704
reference_year,
697-
res_arcsec,
698705
data_dir,
699706
gpw_version,
700707
exponents,
@@ -780,43 +787,34 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec):
780787
target_grid : dict
781788
A dictionary containing metadata for the global target grid.
782789
"""
783-
res_deg = res_arcsec / 3600 # Convert arcseconds to degrees
784-
790+
res_deg = res_arcsec / 3600
785791
if res_arcsec == 15:
786-
# **Use Nightlight grid (NASA Black Marble)**
787-
file_path = nl_util.get_nasa_nl_file_path(reference_year, data_dir=data_dir, verbose=False)
788-
789-
# Load metadata
790-
with rasterio.open(file_path, "r") as src:
791-
global_transform = src.transform
792-
global_crs = src.crs
793-
global_width = src.width
794-
global_height = src.height
795-
796-
elif res_arcsec:
797-
# **Use GPW grid (Population)**
798-
LOGGER.info(f"Loading global GPW metadata for {reference_year} at {res_arcsec} arcsec...")
799-
file_path = pop_util.get_gpw_file_path(gpw_version, reference_year, data_dir=data_dir, verbose=False)
800-
801-
# Load metadata
792+
# Black Marble Nightlights Grid
793+
global_crs = rasterio.crs.CRS.from_epsg(4326) # WGS84 projection
794+
global_transform = Affine(
795+
res_deg, 0, -180,
796+
0, -res_deg, 90
797+
)
798+
global_width = round(360 / res_deg)
799+
global_height = round(180 / res_deg)
800+
else:
801+
# GPW Population Grid
802+
file_path = pop_util.get_gpw_file_path(gpw_version,
803+
reference_year,
804+
data_dir=data_dir,
805+
verbose=False)
802806
with rasterio.open(file_path, "r") as src:
803-
global_crs = src.crs # Keep CRS from GPW dataset
804-
gpw_transform = src.transform # Original GPW grid transform
805-
806-
# **Align the resolution to GPW grid**
807+
global_crs = src.crs
808+
gpw_transform = src.transform
809+
# Align grid resolution with GPW dataset
807810
aligned_lon_min = -180 + (round((gpw_transform[2] - (-180)) / res_deg) * res_deg)
808811
aligned_lat_max = 90 - (round((90 - gpw_transform[5]) / res_deg) * res_deg)
809-
810812
global_transform = Affine(
811-
res_deg, 0, aligned_lon_min, # Adjust longitude alignment
812-
0, -res_deg, aligned_lat_max # Adjust latitude alignment
813-
)
814-
815-
# Compute width & height for the new grid
813+
res_deg, 0, aligned_lon_min,
814+
0, -res_deg, aligned_lat_max)
816815
global_width = round(360 / res_deg)
817816
global_height = round(180 / res_deg)
818-
819-
# 🌍 **Create the global target grid**
817+
# Define the target grid using the computed values
820818
target_grid = {
821819
"driver": "GTiff",
822820
"dtype": "float32",
@@ -826,7 +824,6 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec):
826824
"height": global_height,
827825
"transform": global_transform,
828826
}
829-
830827
return target_grid
831828

832829
@staticmethod
@@ -893,7 +890,6 @@ def _from_country(
893890
) = _get_litpop_single_polygon(
894891
polygon,
895892
reference_year,
896-
res_arcsec,
897893
data_dir,
898894
gpw_version,
899895
exponents,
@@ -965,7 +961,6 @@ def _crop_target_grid(target_grid, polygon):
965961
def _get_litpop_single_polygon(
966962
polygon,
967963
reference_year,
968-
res_arcsec,
969964
data_dir,
970965
gpw_version,
971966
exponents,

0 commit comments

Comments
 (0)