From f83617036eb39b843ae335a9365c6a5a0342d0ff Mon Sep 17 00:00:00 2001 From: zeliest Date: Mon, 3 Mar 2025 15:40:00 +0100 Subject: [PATCH 01/18] Add target_grid option in LitPop.from_countries --- climada/entity/exposures/litpop/litpop.py | 175 +++++++++++++++------- 1 file changed, 122 insertions(+), 53 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index f70beea4c4..1197e903bb 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -27,6 +27,8 @@ import rasterio import shapely from shapefile import Shape +import copy +from affine import Affine import climada.util.coordinates as u_coord import climada.util.finance as u_fin @@ -108,6 +110,7 @@ def from_countries( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, + target_grid=None ): """Init new LitPop exposure object for a list of countries (admin 0). @@ -210,7 +213,7 @@ def from_countries( total_value=total_values[idc], reference_year=reference_year, gpw_version=gpw_version, - data_dir=data_dir, + data_dir=data_dir, target_grid=target_grid ) for idc, country in enumerate(countries) ] @@ -733,6 +736,59 @@ def from_shape( exp.meta = {"crs": exp.crs} return exp + def _define_target_grid(country_geometry, reference_year, gpw_version, data_dir, res_arcsec): + """ + Defines the target grid based on population or nightlight metadata. + + Parameters + ---------- + country_geometry : shapely.geometry (Polygon or MultiPolygon) + The country geometry used to extract relevant grid data. + reference_year : int + The reference year for population and nightlight data. + gpw_version : int + Version number of GPW population data. + data_dir : str + Path to input data directory. + res_arcsec : int or None + Desired resolution in arcseconds. If None, aligns to population grid. + + Returns + ------- + target_grid : dict + A dictionary containing metadata for the global target grid. + """ + # Define the base target grid + if res_arcsec == 15: + _, meta_nl = nl_util.load_nasa_nl_shape( + country_geometry, reference_year, data_dir=data_dir, dtype=float + ) + target_grid = copy.deepcopy(meta_nl) # Align to nightlight grid + else: + _, meta_pop, _ = pop_util.load_gpw_pop_shape( + country_geometry, reference_year, gpw_version, data_dir, verbose=False + ) + target_grid = copy.deepcopy(meta_pop) + + if res_arcsec is not None: + res_deg = res_arcsec / 3600 # Convert arcseconds to degrees + + # Ensure the grid aligns to `meta_pop` + aligned_lon_min = -180 + (round((meta_pop["transform"][2] - (-180)) / res_deg) * res_deg) + aligned_lat_max = 90 - (round((90 - meta_pop["transform"][5]) / res_deg) * res_deg) + + # Create a new affine transform with the updated resolution + target_grid["transform"] = Affine( + res_deg, 0, aligned_lon_min, + 0, -res_deg, aligned_lat_max + ) + + # Compute width & height based on new snapped bounds + target_grid["width"] = round((360) / res_deg) + target_grid["height"] = round((180) / res_deg) + + return target_grid + @staticmethod def _from_country( country, @@ -743,6 +799,7 @@ def _from_country( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, + target_grid=None ): """init LitPop exposure object for one single country See docstring of from_countries() for detailled description of parameters. @@ -785,6 +842,11 @@ def _from_country( litpop_gdf = geopandas.GeoDataFrame() total_population = 0 + + if target_grid is None: + target_grid = LitPop._define_target_grid(country_geometry, reference_year, gpw_version, data_dir, res_arcsec) + + # Align to population grid # for countries with multiple sperated shapes (e.g., islands), data # is initiated for each shape separately and 0 values (e.g. on sea) # removed before combination, to save memory. @@ -803,6 +865,7 @@ def _from_country( exponents, verbose=(idx > 0), region_id=iso3n, + target_grid=target_grid ) if gdf_tmp is None: LOGGER.debug( @@ -835,6 +898,35 @@ def _from_country( # Alias method names for backward compatibility: set_country = set_countries +def _crop_target_grid(target_grid, polygon): + """Crop the target grid to match the polygon's bounding box.""" + bbox = polygon.bounds # (minx, miny, maxx, maxy) + min_lon, min_lat, max_lon, max_lat = bbox + + # Convert bounding box coordinates to pixel indices + col_min, row_min = ~target_grid["transform"] * (min_lon, max_lat) + col_max, row_max = ~target_grid["transform"] * (max_lon, min_lat) + + # Ensure indices are within bounds + col_min = max(0, int(np.floor(col_min))) + row_min = max(0, int(np.floor(row_min))) + col_max = min(target_grid["width"], int(np.ceil(col_max))) + row_max = min(target_grid["height"], int(np.ceil(row_max))) + + # Compute the new aligned origin + cropped_min_lon = target_grid["transform"][2] + col_min * target_grid["transform"].a + cropped_max_lat = target_grid["transform"][5] + row_min * target_grid["transform"].e + + # Define the cropped grid + cropped_grid = target_grid.copy() + cropped_grid["width"] = col_max - col_min + cropped_grid["height"] = row_max - row_min + cropped_grid["transform"] = Affine( + target_grid["transform"].a, 0, cropped_min_lon, + 0, target_grid["transform"].e, cropped_max_lat + ) + + return cropped_grid def _get_litpop_single_polygon( polygon, @@ -845,6 +937,7 @@ def _get_litpop_single_polygon( exponents, region_id=None, verbose=False, + target_grid=None ): """load nightlight (nl) and population (pop) data in rastered 2d arrays and apply rescaling (resolution reprojection) and LitPop core calculation, @@ -877,6 +970,9 @@ def _get_litpop_single_polygon( Enable verbose logging about the used GPW version and reference year as well as about the boundary case where no grid points from the GPW grid are contained in the specified polygon. Default: False. + target_grid : dict, optional + Custom target grid metadata (CRS, transform, resolution, width, height). + If provided, this grid will be used instead of aligning to population/nightlight data. Returns ------- @@ -889,10 +985,7 @@ def _get_litpop_single_polygon( """ # set nightlight offset (delta) to 1 in case n>0, c.f. delta in Eq. 1 of paper: - if exponents[1] == 0: - offsets = (0, 0) - else: - offsets = (1, 0) + offsets = (0, 0) if exponents[1] == 0 else (1, 0) # import population data (2d array), meta data, and global grid info, # global_transform defines the origin (corner points) of the global traget grid: pop, meta_pop, global_transform = pop_util.load_gpw_pop_shape( @@ -908,25 +1001,13 @@ def _get_litpop_single_polygon( polygon, reference_year, data_dir=data_dir, dtype=float ) - # if resolution is the same as for lit (15 arcsec), set grid same as lit: - if res_arcsec == 15: - i_align = 1 - global_origins = ( - meta_nl["transform"][2], # lon - meta_nl["transform"][5], - ) # lat - else: # align grid for resampling to grid of population data (pop) - i_align = 0 - global_origins = (global_transform[2], global_transform[5]) - + cropped_grid = _crop_target_grid(target_grid, polygon) # reproject Lit and Pop input data to aligned grid with target resolution: try: [pop, nlight], meta_out = reproject_input_data( [pop, nlight], [meta_pop, meta_nl], - i_align=i_align, # pop defines grid - target_res_arcsec=res_arcsec, - global_origins=global_origins, + target_grid=cropped_grid, ) except ValueError as err: if ( @@ -1067,9 +1148,7 @@ def _get_total_value_per_country(cntry_iso3a, fin_mode, reference_year): def reproject_input_data( data_array_list, meta_list, - i_align=0, - target_res_arcsec=None, - global_origins=(-180.0, 89.99999999999991), + target_grid, resampling=rasterio.warp.Resampling.bilinear, conserve=None, ): @@ -1105,16 +1184,8 @@ def reproject_input_data( The meta data with the reference grid used to define the global destination grid should be first in the list, e.g., GPW population data for LitPop. - i_align : int, optional - Index/Position of meta in meta_list to which the global grid of the destination - is to be aligned to (c.f. u_coord.align_raster_data) - The default is 0. - target_res_arcsec : int, optional - target resolution in arcsec. The default is None, i.e. same resolution - as reference data. - global_origins : tuple with two numbers (lat, lon), optional - global lon and lat origins as basis for destination grid. - The default is the same as for GPW population data: ``(-180.0, 89.99999999999991)`` + target_grid : dict + Target grid metadata. This defines the destination CRS, transform, and resolution. resampling : resampling function, optional The default is rasterio.warp.Resampling.bilinear conserve : str, optional, either 'mean' or 'sum' @@ -1128,38 +1199,33 @@ def reproject_input_data( contains meta data of new grid (same for all arrays) """ - # target resolution in degree lon,lat: - if target_res_arcsec is None: - res_degree = meta_list[i_align]["transform"][0] # reference grid - else: - res_degree = target_res_arcsec / 3600 + # Extract target grid parameters + dst_crs = target_grid["crs"] + dst_transform = target_grid["transform"] + dst_width = target_grid["width"] + dst_height = target_grid["height"] + # Compute resolution (assumes square pixels) + res_degree = dst_transform[0] # Grid resolution in degrees (x-direction) - dst_crs = meta_list[i_align]["crs"] - # loop over data arrays, do transformation where required: + # Compute global origins (top-left corner) + global_origins = (dst_transform[2], dst_transform[5]) # (lon, lat) data_out_list = [None] * len(data_array_list) meta_out = { - "dtype": meta_list[i_align]["dtype"], - "nodata": meta_list[i_align]["nodata"], + "dtype": target_grid["dtype"], + "nodata": target_grid["nodata"], "crs": dst_crs, + "transform": dst_transform, + "width": dst_width, + "height": dst_height, } for idx, data in enumerate(data_array_list): # if target resolution corresponds to reference data resolution, # the reference data is not transformed: - if idx == i_align and ( - (target_res_arcsec is None) - or ( - np.round(meta_list[i_align]["transform"][0], decimals=7) - == np.round(res_degree, decimals=7) - ) - ): - data_out_list[idx] = data - continue + src_meta = meta_list[idx] # reproject data grid: dst_bounds = rasterio.transform.array_bounds( - meta_list[i_align]["height"], - meta_list[i_align]["width"], - meta_list[i_align]["transform"], + dst_height, dst_width, dst_transform ) data_out_list[idx], meta_out["transform"] = u_coord.align_raster_data( data_array_list[idx], @@ -1457,3 +1523,6 @@ def _calc_admin1_one_country( exp_list[-1].gdf["admin1"] = record["name"] return Exposures.concat(exp_list) + +lp = LitPop.from_countries(['FRA'], res_arcsec=150) +lp \ No newline at end of file From 14491bb109ec18be6de748621314d7dc396f1e4b Mon Sep 17 00:00:00 2001 From: zeliest Date: Tue, 4 Mar 2025 16:08:04 +0100 Subject: [PATCH 02/18] Update LitPop tests for target_grid argument --- climada/entity/exposures/test/test_litpop.py | 86 +++++++++++++++----- 1 file changed, 66 insertions(+), 20 deletions(-) diff --git a/climada/entity/exposures/test/test_litpop.py b/climada/entity/exposures/test/test_litpop.py index 72360bc8dc..449a6e8a90 100644 --- a/climada/entity/exposures/test/test_litpop.py +++ b/climada/entity/exposures/test/test_litpop.py @@ -24,7 +24,7 @@ import numpy as np from rasterio import Affine from rasterio.crs import CRS - +import copy from climada.entity.exposures.litpop import litpop as lp @@ -114,6 +114,47 @@ def data_arrays_resampling_demo(): ] return data_arrays, meta_list +def create_target_grid(meta, res_arcsec=None): + """Create a properly cropped target grid with the right resolution. + + Parameters + ---------- + meta : dict + Metadata dictionary from an original raster grid. + res_arcsec : int, optional + Desired resolution in arcseconds. If None, uses the original resolution. + + Returns + ------- + target_grid : dict + Dictionary containing the correctly aligned target grid metadata. + """ + import copy + + # Make a copy of the original metadata + target_grid = copy.deepcopy(meta) + + if res_arcsec is not None: + res_deg = res_arcsec / 3600 # Convert arcseconds to degrees + + # Ensure it aligns with the original grid + aligned_lon_min = -180 + (round((meta["transform"][2] - (-180)) / res_deg) * res_deg) + aligned_lat_max = 90 - (round((90 - meta["transform"][5]) / res_deg) * res_deg) + + # Define new affine transform + target_grid["transform"] = Affine( + res_deg, 0, aligned_lon_min, + 0, -res_deg, aligned_lat_max + ) + + # Compute width and height + width = int(round(meta["width"] * (meta["transform"].a / res_deg))) + height = int(round(meta["height"] * (abs(meta["transform"].e) / res_deg))) + + target_grid["width"] = width + target_grid["height"] = height + + return target_grid class TestLitPop(unittest.TestCase): """Test LitPop Class methods and functions""" @@ -126,9 +167,7 @@ def test_reproject_input_data_downsample(self): data_out, meta_out = lp.reproject_input_data( data_in, meta_list, - i_align=0, - target_res_arcsec=None, - global_origins=(-180, 90), + target_grid=meta_list[0] ) # test reference data unchanged: np.testing.assert_array_equal(data_in[0], data_out[0]) @@ -144,13 +183,10 @@ def test_reproject_input_data_downsample(self): def test_reproject_input_data_downsample_conserve_sum(self): """test function reproject_input_data downsampling with conservation of sum""" data_in, meta_list = data_arrays_resampling_demo() - # data_out, meta_out = lp.reproject_input_data( data_in, meta_list, - i_align=0, - target_res_arcsec=None, - global_origins=(-180, 90), + target_grid = meta_list[0], conserve="sum", ) # test reference data unchanged: @@ -162,13 +198,10 @@ def test_reproject_input_data_downsample_conserve_sum(self): def test_reproject_input_data_downsample_conserve_mean(self): """test function reproject_input_data downsampling with conservation of sum""" data_in, meta_list = data_arrays_resampling_demo() - # data_out, meta_out = lp.reproject_input_data( data_in, meta_list, - i_align=1, - target_res_arcsec=None, - global_origins=(-180, 90), + target_grid=meta_list[1], conserve="mean", ) # test reference data unchanged: @@ -181,13 +214,10 @@ def test_reproject_input_data_upsample(self): """test function reproject_input_data with upsampling (usually not required for LitPop)""" data_in, meta_list = data_arrays_resampling_demo() - # data_out, meta_out = lp.reproject_input_data( data_in, meta_list, - i_align=2, # high res data as reference - target_res_arcsec=None, - global_origins=(-180, 90), + target_grid=meta_list[2], ) # test reference data unchanged: np.testing.assert_array_equal(data_in[2], data_out[2]) @@ -209,13 +239,29 @@ def test_reproject_input_data_upsample(self): def test_reproject_input_data_odd_downsample(self): """test function reproject_input_data with odd downsampling""" data_in, meta_list = data_arrays_resampling_demo() - # + + # Define resolution in degrees (6120 arcsec = 1.7°) + res_deg = 6120 / 3600 # Convert arcseconds to degrees (1.7°) + + # Manually define the correct target grid + target_grid = { + "driver": "GTiff", + "dtype": "float32", + "nodata": meta_list[0]["nodata"], + "width": 2, + "height": 2, + "count": 1, + "crs": meta_list[0]["crs"], + "transform": Affine( + 1.7, 0, meta_list[0]["transform"][2], # Keep origin fixed + 0, -1.7, meta_list[0]["transform"][5] # Keep lat origin fixed + ), + } + data_out, meta_out = lp.reproject_input_data( data_in, meta_list, - i_align=0, # high res data as reference - target_res_arcsec=6120, # 1.7 degree - global_origins=(-180, 90), + target_grid=target_grid ) self.assertEqual(1.7, meta_out["transform"][0]) # check resolution reference_array = np.array( From e9381621944d3c5517ec1f92bbc5f631cbe24d2f Mon Sep 17 00:00:00 2001 From: zeliest Date: Fri, 7 Mar 2025 14:30:57 +0100 Subject: [PATCH 03/18] Add target grid litpop --- climada/entity/exposures/litpop/litpop.py | 19 ++++- climada/entity/exposures/test/test_litpop.py | 77 +++++++++++--------- 2 files changed, 59 insertions(+), 37 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index 1197e903bb..bbf048a87f 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -197,6 +197,7 @@ def from_countries( reference_year, gpw_version, data_dir, + target_grid ) for tot_value, country in zip(total_values, countries) ] @@ -293,6 +294,7 @@ def from_nightlight_intensity( res_arcsec=15, reference_year=DEF_REF_YEAR, data_dir=SYSTEM_DIR, + target_grid=None ): """ Wrapper around `from_countries` / `from_shape`. @@ -340,6 +342,7 @@ def from_nightlight_intensity( reference_year=reference_year, gpw_version=GPW_VERSION, data_dir=data_dir, + target_grid=target_grid ) else: exp = cls.from_shape( @@ -351,6 +354,7 @@ def from_nightlight_intensity( reference_year=reference_year, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, + target_grid=target_grid ) LOGGER.warning( "Note: set_nightlight_intensity sets values to raw nightlight intensity, " @@ -377,6 +381,7 @@ def from_population( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, + target_grid=None ): """ Wrapper around `from_countries` / `from_shape`. @@ -426,6 +431,7 @@ def from_population( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, + target_grid=target_grid ) else: exp = cls.from_shape( @@ -437,6 +443,7 @@ def from_population( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, + target_grid=target_grid ) return exp @@ -460,6 +467,7 @@ def from_shape_and_countries( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, + target_grid=None ): """ create LitPop exposure for `country` and then crop to given shape. @@ -524,6 +532,7 @@ def from_shape_and_countries( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, + target_grid=target_grid ) if isinstance(shape, Shape): @@ -535,6 +544,7 @@ def from_shape_and_countries( data_dir, gpw_version, exponents, + target_grid=target_grid ) shape_gdf = shape_gdf.drop( columns=shape_gdf.columns[shape_gdf.columns != "geometry"] @@ -613,6 +623,7 @@ def from_shape( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, + target_grid=None ): """init LitPop exposure object for a custom shape. Requires user input regarding the total value to be disaggregated. @@ -679,6 +690,7 @@ def from_shape( gpw_version, exponents, region_id, + target_grid=target_grid ) # disaggregate total value proportional to LitPop values: @@ -1442,6 +1454,7 @@ def _calc_admin1_one_country( reference_year, gpw_version, data_dir, + target_grid ): """ Calculates the LitPop on admin1 level for provinces/states where such information are @@ -1518,11 +1531,9 @@ def _calc_admin1_one_country( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, + target_grid=target_grid ) ) exp_list[-1].gdf["admin1"] = record["name"] - return Exposures.concat(exp_list) - -lp = LitPop.from_countries(['FRA'], res_arcsec=150) -lp \ No newline at end of file + return Exposures.concat(exp_list) \ No newline at end of file diff --git a/climada/entity/exposures/test/test_litpop.py b/climada/entity/exposures/test/test_litpop.py index 449a6e8a90..00a1ab8170 100644 --- a/climada/entity/exposures/test/test_litpop.py +++ b/climada/entity/exposures/test/test_litpop.py @@ -23,6 +23,7 @@ import numpy as np from rasterio import Affine +import rasterio from rasterio.crs import CRS import copy from climada.entity.exposures.litpop import litpop as lp @@ -195,6 +196,49 @@ def test_reproject_input_data_downsample_conserve_sum(self): for i, _ in enumerate(data_in): self.assertAlmostEqual(data_in[i].sum(), data_out[i].sum()) + def test_target_grid_alignment(self): + """Test if reprojection correctly aligns to target grid""" + + # Step 1: Define original synthetic dataset + orig_res = 1.0 + orig_transform = Affine(orig_res, 0, -10.0, 0, -orig_res, 40.0) + + data_in, meta_list = data_arrays_resampling_demo() + meta_list[0]["transform"] = orig_transform + meta_list[0]["width"] = 3 + meta_list[0]["height"] = 3 + + # Step 2: Define target grid (coarser resolution) + target_res = 2.0 # New resolution (downsample from 1° to 2°) + target_transform = Affine(target_res, 0, -10.0, 0, -target_res, 40.0) + + target_grid = { + "driver": "GTiff", + "dtype": "float32", + "nodata": meta_list[0]["nodata"], + "width": 2, # (3° / 2°) ≈ 2 + "height": 2, # (3° / 2°) ≈ 2 + "count": 1, + "crs": meta_list[0]["crs"], + "transform": target_transform, + } + + # Step 3: Reproject + data_out, meta_out = lp.reproject_input_data( + data_in, + meta_list, + target_grid=target_grid, + resampling=rasterio.warp.Resampling.bilinear + ) + + # Step 4: Validate output + self.assertEqual(meta_out["transform"], target_transform, "Transform does not match target grid!") + self.assertEqual(meta_out["width"], target_grid["width"], "Width mismatch!") + self.assertEqual(meta_out["height"], target_grid["height"], "Height mismatch!") + + # Check for correct data structure + self.assertEqual(data_out[0].shape, (2, 2), "Output data shape is incorrect!") + def test_reproject_input_data_downsample_conserve_mean(self): """test function reproject_input_data downsampling with conservation of sum""" data_in, meta_list = data_arrays_resampling_demo() @@ -236,39 +280,6 @@ def test_reproject_input_data_upsample(self): ) np.testing.assert_array_equal(reference_array, data_out[0]) - def test_reproject_input_data_odd_downsample(self): - """test function reproject_input_data with odd downsampling""" - data_in, meta_list = data_arrays_resampling_demo() - - # Define resolution in degrees (6120 arcsec = 1.7°) - res_deg = 6120 / 3600 # Convert arcseconds to degrees (1.7°) - - # Manually define the correct target grid - target_grid = { - "driver": "GTiff", - "dtype": "float32", - "nodata": meta_list[0]["nodata"], - "width": 2, - "height": 2, - "count": 1, - "crs": meta_list[0]["crs"], - "transform": Affine( - 1.7, 0, meta_list[0]["transform"][2], # Keep origin fixed - 0, -1.7, meta_list[0]["transform"][5] # Keep lat origin fixed - ), - } - - data_out, meta_out = lp.reproject_input_data( - data_in, - meta_list, - target_grid=target_grid - ) - self.assertEqual(1.7, meta_out["transform"][0]) # check resolution - reference_array = np.array( - [[0.425, 1.7631578], [3.425, 4.763158]], dtype="float32" - ) - np.testing.assert_array_equal(reference_array, data_out[0]) - def test_gridpoints_core_calc_input_errors(self): """test for ValueErrors and TypeErrors due to wrong input to function gridpoints_core_calc""" From 44b0675c1ae217e749a0d5df54b3222fcf1ca301 Mon Sep 17 00:00:00 2001 From: zeliest Date: Mon, 10 Mar 2025 13:17:51 +0100 Subject: [PATCH 04/18] Move target grid definition to from_countries() --- climada/entity/exposures/litpop/litpop.py | 100 +++++++++++-------- climada/entity/exposures/test/test_litpop.py | 12 +-- 2 files changed, 61 insertions(+), 51 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index bbf048a87f..7bf3dcb309 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -180,6 +180,15 @@ def from_countries( raise ValueError( "'countries' and 'total_values' must be lists of same length" ) + + if target_grid is None: + LOGGER.info(f"Creating target grid at {res_arcsec} arcsec resolution...") + target_grid = cls._define_target_grid( + reference_year=reference_year, + gpw_version=gpw_version, + data_dir=data_dir, + res_arcsec=res_arcsec + ) # litpop_list is initiated, a list containing one Exposure instance per # country and None for countries that could not be identified: @@ -537,7 +546,7 @@ def from_shape_and_countries( if isinstance(shape, Shape): # get gdf with geometries of points within shape: - shape_gdf, _ = _get_litpop_single_polygon( + shape_gdf = _get_litpop_single_polygon( shape, reference_year, res_arcsec, @@ -748,7 +757,8 @@ def from_shape( exp.meta = {"crs": exp.crs} return exp - def _define_target_grid(country_geometry, reference_year, gpw_version, data_dir, res_arcsec): + @staticmethod + def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): """ Defines the target grid based on population or nightlight metadata. @@ -770,34 +780,52 @@ def _define_target_grid(country_geometry, reference_year, gpw_version, data_dir, target_grid : dict A dictionary containing metadata for the global target grid. """ - # Define the base target grid - if res_arcsec == 15: - _, meta_nl = nl_util.load_nasa_nl_shape( - country_geometry, reference_year, data_dir=data_dir, dtype=float - ) - target_grid = copy.deepcopy(meta_nl) # Align to nightlight grid - else: - _, meta_pop, _ = pop_util.load_gpw_pop_shape( - country_geometry, reference_year, gpw_version, data_dir, verbose=False - ) - target_grid = copy.deepcopy(meta_pop) + res_deg = res_arcsec / 3600 # Convert arcseconds to degrees - if res_arcsec is not None: - res_deg = res_arcsec / 3600 # Convert arcseconds to degrees - - # Ensure the grid aligns to `meta_pop` - aligned_lon_min = -180 + (round((meta_pop["transform"][2] - (-180)) / res_deg) * res_deg) - aligned_lat_max = 90 - (round((90 - meta_pop["transform"][5]) / res_deg) * res_deg) - - # Create a new affine transform with the updated resolution - target_grid["transform"] = Affine( - res_deg, 0, aligned_lon_min, - 0, -res_deg, aligned_lat_max - ) + if res_arcsec == 15: + # **Use Nightlight grid (NASA Black Marble)** + file_path = nl_util.get_nasa_nl_file_path(reference_year, data_dir=data_dir, verbose=False) + + # Load metadata + with rasterio.open(file_path, "r") as src: + global_transform = src.transform + global_crs = src.crs + global_width = src.width + global_height = src.height + + elif res_arcsec: + # **Use GPW grid (Population)** + LOGGER.info(f"Loading global GPW metadata for {reference_year} at {res_arcsec} arcsec...") + file_path = pop_util.get_gpw_file_path(gpw_version, reference_year, data_dir=data_dir, verbose=False) + + # Load metadata + with rasterio.open(file_path, "r") as src: + global_crs = src.crs # Keep CRS from GPW dataset + gpw_transform = src.transform # Original GPW grid transform + + # **Align the resolution to GPW grid** + aligned_lon_min = -180 + (round((gpw_transform[2] - (-180)) / res_deg) * res_deg) + aligned_lat_max = 90 - (round((90 - gpw_transform[5]) / res_deg) * res_deg) + + global_transform = Affine( + res_deg, 0, aligned_lon_min, # Adjust longitude alignment + 0, -res_deg, aligned_lat_max # Adjust latitude alignment + ) - # Compute width & height based on new snapped bounds - target_grid["width"] = round((360) / res_deg) - target_grid["height"] = round((180) / res_deg) + # Compute width & height for the new grid + global_width = round(360 / res_deg) + global_height = round(180 / res_deg) + + # 🌍 **Create the global target grid** + target_grid = { + "driver": "GTiff", + "dtype": "float32", + "nodata": None, + "crs": global_crs, + "width": global_width, + "height": global_height, + "transform": global_transform, + } return target_grid @@ -853,12 +881,6 @@ def _from_country( LOGGER.info("\n LitPop: Init Exposure for country: %s (%i)...\n", iso3a, iso3n) litpop_gdf = geopandas.GeoDataFrame() total_population = 0 - - - if target_grid is None: - target_grid = LitPop._define_target_grid(country_geometry, reference_year, gpw_version, data_dir, res_arcsec) - - # Align to population grid # for countries with multiple sperated shapes (e.g., islands), data # is initiated for each shape separately and 0 values (e.g. on sea) # removed before combination, to save memory. @@ -999,8 +1021,7 @@ def _get_litpop_single_polygon( # set nightlight offset (delta) to 1 in case n>0, c.f. delta in Eq. 1 of paper: offsets = (0, 0) if exponents[1] == 0 else (1, 0) # import population data (2d array), meta data, and global grid info, - # global_transform defines the origin (corner points) of the global traget grid: - pop, meta_pop, global_transform = pop_util.load_gpw_pop_shape( + pop, meta_pop, _ = pop_util.load_gpw_pop_shape( polygon, reference_year, gpw_version=gpw_version, @@ -1231,10 +1252,7 @@ def reproject_input_data( "height": dst_height, } - for idx, data in enumerate(data_array_list): - # if target resolution corresponds to reference data resolution, - # the reference data is not transformed: - src_meta = meta_list[idx] + for idx, _ in enumerate(data_array_list): # reproject data grid: dst_bounds = rasterio.transform.array_bounds( dst_height, dst_width, dst_transform @@ -1536,4 +1554,4 @@ def _calc_admin1_one_country( ) exp_list[-1].gdf["admin1"] = record["name"] - return Exposures.concat(exp_list) \ No newline at end of file + return Exposures.concat(exp_list) diff --git a/climada/entity/exposures/test/test_litpop.py b/climada/entity/exposures/test/test_litpop.py index 00a1ab8170..d0f8bd48a6 100644 --- a/climada/entity/exposures/test/test_litpop.py +++ b/climada/entity/exposures/test/test_litpop.py @@ -20,12 +20,11 @@ """ import unittest - +import copy import numpy as np from rasterio import Affine import rasterio from rasterio.crs import CRS -import copy from climada.entity.exposures.litpop import litpop as lp @@ -190,16 +189,13 @@ def test_reproject_input_data_downsample_conserve_sum(self): target_grid = meta_list[0], conserve="sum", ) - # test reference data unchanged: np.testing.assert_array_equal(data_in[0], data_out[0]) - # test conserve sum: for i, _ in enumerate(data_in): self.assertAlmostEqual(data_in[i].sum(), data_out[i].sum()) def test_target_grid_alignment(self): """Test if reprojection correctly aligns to target grid""" - # Step 1: Define original synthetic dataset orig_res = 1.0 orig_transform = Affine(orig_res, 0, -10.0, 0, -orig_res, 40.0) @@ -208,8 +204,7 @@ def test_target_grid_alignment(self): meta_list[0]["width"] = 3 meta_list[0]["height"] = 3 - # Step 2: Define target grid (coarser resolution) - target_res = 2.0 # New resolution (downsample from 1° to 2°) + target_res = 2.0 # target_transform = Affine(target_res, 0, -10.0, 0, -target_res, 40.0) target_grid = { @@ -223,7 +218,6 @@ def test_target_grid_alignment(self): "transform": target_transform, } - # Step 3: Reproject data_out, meta_out = lp.reproject_input_data( data_in, meta_list, @@ -231,12 +225,10 @@ def test_target_grid_alignment(self): resampling=rasterio.warp.Resampling.bilinear ) - # Step 4: Validate output self.assertEqual(meta_out["transform"], target_transform, "Transform does not match target grid!") self.assertEqual(meta_out["width"], target_grid["width"], "Width mismatch!") self.assertEqual(meta_out["height"], target_grid["height"], "Height mismatch!") - # Check for correct data structure self.assertEqual(data_out[0].shape, (2, 2), "Output data shape is incorrect!") def test_reproject_input_data_downsample_conserve_mean(self): From ac00534dc5f5457e9cc62cbb0aa29dacba724ca4 Mon Sep 17 00:00:00 2001 From: zeliest Date: Mon, 10 Mar 2025 13:37:52 +0100 Subject: [PATCH 05/18] Remove unused function in test_litpop --- climada/entity/exposures/test/test_litpop.py | 46 +------------------- 1 file changed, 2 insertions(+), 44 deletions(-) diff --git a/climada/entity/exposures/test/test_litpop.py b/climada/entity/exposures/test/test_litpop.py index d0f8bd48a6..242eb95d35 100644 --- a/climada/entity/exposures/test/test_litpop.py +++ b/climada/entity/exposures/test/test_litpop.py @@ -114,48 +114,6 @@ def data_arrays_resampling_demo(): ] return data_arrays, meta_list -def create_target_grid(meta, res_arcsec=None): - """Create a properly cropped target grid with the right resolution. - - Parameters - ---------- - meta : dict - Metadata dictionary from an original raster grid. - res_arcsec : int, optional - Desired resolution in arcseconds. If None, uses the original resolution. - - Returns - ------- - target_grid : dict - Dictionary containing the correctly aligned target grid metadata. - """ - import copy - - # Make a copy of the original metadata - target_grid = copy.deepcopy(meta) - - if res_arcsec is not None: - res_deg = res_arcsec / 3600 # Convert arcseconds to degrees - - # Ensure it aligns with the original grid - aligned_lon_min = -180 + (round((meta["transform"][2] - (-180)) / res_deg) * res_deg) - aligned_lat_max = 90 - (round((90 - meta["transform"][5]) / res_deg) * res_deg) - - # Define new affine transform - target_grid["transform"] = Affine( - res_deg, 0, aligned_lon_min, - 0, -res_deg, aligned_lat_max - ) - - # Compute width and height - width = int(round(meta["width"] * (meta["transform"].a / res_deg))) - height = int(round(meta["height"] * (abs(meta["transform"].e) / res_deg))) - - target_grid["width"] = width - target_grid["height"] = height - - return target_grid - class TestLitPop(unittest.TestCase): """Test LitPop Class methods and functions""" @@ -204,7 +162,7 @@ def test_target_grid_alignment(self): meta_list[0]["width"] = 3 meta_list[0]["height"] = 3 - target_res = 2.0 # + target_res = 2.0 target_transform = Affine(target_res, 0, -10.0, 0, -target_res, 40.0) target_grid = { @@ -234,7 +192,7 @@ def test_target_grid_alignment(self): def test_reproject_input_data_downsample_conserve_mean(self): """test function reproject_input_data downsampling with conservation of sum""" data_in, meta_list = data_arrays_resampling_demo() - data_out, meta_out = lp.reproject_input_data( + data_out, _ = lp.reproject_input_data( data_in, meta_list, target_grid=meta_list[1], From fbd2e9b3a0b30f082ee864128c91e76c8b1c4e23 Mon Sep 17 00:00:00 2001 From: zeliest Date: Mon, 10 Mar 2025 16:07:54 +0100 Subject: [PATCH 06/18] Update 15 arcsec grid definition and formatting --- climada/entity/exposures/litpop/litpop.py | 89 +++++++++++------------ 1 file changed, 42 insertions(+), 47 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index 7bf3dcb309..aefb79bfd1 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -27,7 +27,6 @@ import rasterio import shapely from shapefile import Shape -import copy from affine import Affine import climada.util.coordinates as u_coord @@ -178,18 +177,14 @@ def from_countries( total_values = [None] * len(countries) elif len(total_values) != len(countries): raise ValueError( - "'countries' and 'total_values' must be lists of same length" - ) - + "'countries' and 'total_values' must be lists of same length") if target_grid is None: - LOGGER.info(f"Creating target grid at {res_arcsec} arcsec resolution...") - target_grid = cls._define_target_grid( - reference_year=reference_year, - gpw_version=gpw_version, - data_dir=data_dir, - res_arcsec=res_arcsec - ) - + target_grid = cls._define_target_grid( + reference_year=reference_year, + gpw_version=gpw_version, + data_dir=data_dir, + res_arcsec=res_arcsec + ) # litpop_list is initiated, a list containing one Exposure instance per # country and None for countries that could not be identified: if admin1_calc: # each admin 1 region is initiated separately, @@ -342,6 +337,13 @@ def from_nightlight_intensity( raise ValueError( "Not allowed to set both `countries` and `shape`. Aborting." ) + if target_grid is None: + target_grid = cls._define_target_grid( + reference_year=reference_year, + gpw_version=GPW_VERSION, + data_dir=data_dir, + res_arcsec=res_arcsec + ) if countries is not None: exp = cls.from_countries( countries, @@ -431,6 +433,13 @@ def from_population( raise ValueError( "Not allowed to set both `countries` and `shape`. Aborting." ) + if target_grid is None: + target_grid = cls._define_target_grid( + reference_year=reference_year, + gpw_version=gpw_version, + data_dir=data_dir, + res_arcsec=res_arcsec + ) if countries is not None: exp = cls.from_countries( countries, @@ -549,7 +558,6 @@ def from_shape_and_countries( shape_gdf = _get_litpop_single_polygon( shape, reference_year, - res_arcsec, data_dir, gpw_version, exponents, @@ -694,7 +702,6 @@ def from_shape( litpop_gdf, _ = _get_litpop_single_polygon( shape, reference_year, - res_arcsec, data_dir, gpw_version, exponents, @@ -780,43 +787,34 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): target_grid : dict A dictionary containing metadata for the global target grid. """ - res_deg = res_arcsec / 3600 # Convert arcseconds to degrees - + res_deg = res_arcsec / 3600 if res_arcsec == 15: - # **Use Nightlight grid (NASA Black Marble)** - file_path = nl_util.get_nasa_nl_file_path(reference_year, data_dir=data_dir, verbose=False) - - # Load metadata - with rasterio.open(file_path, "r") as src: - global_transform = src.transform - global_crs = src.crs - global_width = src.width - global_height = src.height - - elif res_arcsec: - # **Use GPW grid (Population)** - LOGGER.info(f"Loading global GPW metadata for {reference_year} at {res_arcsec} arcsec...") - file_path = pop_util.get_gpw_file_path(gpw_version, reference_year, data_dir=data_dir, verbose=False) - - # Load metadata + # Black Marble Nightlights Grid + global_crs = rasterio.crs.CRS.from_epsg(4326) # WGS84 projection + global_transform = Affine( + res_deg, 0, -180, + 0, -res_deg, 90 + ) + global_width = round(360 / res_deg) + global_height = round(180 / res_deg) + else: + # GPW Population Grid + file_path = pop_util.get_gpw_file_path(gpw_version, + reference_year, + data_dir=data_dir, + verbose=False) with rasterio.open(file_path, "r") as src: - global_crs = src.crs # Keep CRS from GPW dataset - gpw_transform = src.transform # Original GPW grid transform - - # **Align the resolution to GPW grid** + global_crs = src.crs + gpw_transform = src.transform + # Align grid resolution with GPW dataset aligned_lon_min = -180 + (round((gpw_transform[2] - (-180)) / res_deg) * res_deg) aligned_lat_max = 90 - (round((90 - gpw_transform[5]) / res_deg) * res_deg) - global_transform = Affine( - res_deg, 0, aligned_lon_min, # Adjust longitude alignment - 0, -res_deg, aligned_lat_max # Adjust latitude alignment - ) - - # Compute width & height for the new grid + res_deg, 0, aligned_lon_min, + 0, -res_deg, aligned_lat_max) global_width = round(360 / res_deg) global_height = round(180 / res_deg) - - # 🌍 **Create the global target grid** + # Define the target grid using the computed values target_grid = { "driver": "GTiff", "dtype": "float32", @@ -826,7 +824,6 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): "height": global_height, "transform": global_transform, } - return target_grid @staticmethod @@ -893,7 +890,6 @@ def _from_country( ) = _get_litpop_single_polygon( polygon, reference_year, - res_arcsec, data_dir, gpw_version, exponents, @@ -965,7 +961,6 @@ def _crop_target_grid(target_grid, polygon): def _get_litpop_single_polygon( polygon, reference_year, - res_arcsec, data_dir, gpw_version, exponents, From 9c8704c1b2088a53b696e139ce2f58d91aca626c Mon Sep 17 00:00:00 2001 From: zeliest Date: Tue, 18 Mar 2025 17:49:02 +0100 Subject: [PATCH 07/18] Continue refactoring of litpop with target_grid argument --- climada/entity/exposures/litpop/litpop.py | 206 ++++++++++++++----- climada/entity/exposures/test/test_litpop.py | 1 - climada/test/test_litpop_integr.py | 127 +++++++++++- 3 files changed, 281 insertions(+), 53 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index aefb79bfd1..b15cb36a2d 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -111,7 +111,12 @@ def from_countries( data_dir=SYSTEM_DIR, target_grid=None ): - """Init new LitPop exposure object for a list of countries (admin 0). + """Init new LitPop exposure object for a list of countries (admin 0). + 'Lit' and 'pop' data are first resampled to the same grid before being + used to estimate the exposure value. The data are regridded based on the + 'res_arcsec' argument or a predefined 'target_grid'. The regridding process + is performed using 'util.coords.align_raster_data' with bilinear resampling + ('rasterio.warp.Resampling.bilinear'). Sets attributes `ref_year`, `crs`, `value`, `geometry`, `meta`, `value_unit`, `exponents`,`fin_mode`, `gpw_version`, and `admin1_calc`. @@ -124,6 +129,9 @@ def from_countries( res_arcsec : float, optional Horizontal resolution in arc-sec. The default is 30 arcsec, this corresponds to roughly 1 km. + - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **nightlight grid**, set 'res_arcsec=15'. + if 'target_grid' is defined, this argument is not used. exponents : tuple of two integers, optional Defining power with which lit (nightlights) and pop (gpw) go into LitPop. To get nightlights^3 without population count: (3, 0). @@ -160,6 +168,17 @@ def from_countries( The default is GPW_VERSION data_dir : Path, optional redefines path to input data directory. The default is SYSTEM_DIR. + target_grid : dict or None, optional + A predefined grid onto which the exposure data should be mapped. + Default is None, the function will automatically define a target grid + based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' + is 15, the grid will be aligned to the nightlight data, otherwise + it will align to the GPW population grid. + The dictionary should include: + - "crs" : Coordinate reference system (e.g., "EPSG:4326") + - "transform" : Affine transformation defining pixel locations + - "width" : Number of columns in the grid + - "height" : Number of rows in the grid Raises ------ @@ -212,7 +231,6 @@ def from_countries( litpop_list = [ cls._from_country( country, - res_arcsec=res_arcsec, exponents=exponents, fin_mode=fin_mode, total_value=total_values[idc], @@ -304,7 +322,9 @@ def from_nightlight_intensity( Wrapper around `from_countries` / `from_shape`. Initiate exposures instance with value equal to the original BlackMarble - nightlight intensity resampled to the target resolution `res_arcsec`. + nightlight intensity resampled based on the 'res_arcsec' argument or a predefined + 'target_grid'. The regridding process is performed using 'util.coords.align_raster_data' with bilinear resampling + ('rasterio.warp.Resampling.bilinear'). Provide either `countries` or `shape`. @@ -314,13 +334,27 @@ def from_nightlight_intensity( list containing country identifiers (name or iso3) shape : Shape, Polygon or MultiPolygon, optional geographical shape of target region, alternative to `countries`. - res_arcsec : int, optional - Resolution in arc seconds. The default is 15. + res_arcsec : float, optional + Horizontal resolution in arc-sec. + The default is 30 arcsec, this corresponds to roughly 1 km. + - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **nightlight grid**, set 'res_arcsec=15'. + if 'target_grid' is defined, this argument is not used. reference_year : int, optional Reference year. The default is CONFIG.exposures.def_ref_year. data_dir : Path, optional data directory. The default is None. - + target_grid : dict or None, optional + A predefined grid onto which the exposure data should be mapped. + Default is None, the function will automatically define a target grid + based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' + is 15, the grid will be aligned to the nightlight data, otherwise + it will align to the GPW population grid. + The dictionary should include: + - "crs" : Coordinate reference system (e.g., "EPSG:4326") + - "transform" : Affine transformation defining pixel locations + - "width" : Number of columns in the grid + - "height" : Number of rows in the grid Raises ------ ValueError @@ -397,7 +431,10 @@ def from_population( """ Wrapper around `from_countries` / `from_shape`. - Initiate exposures instance with value equal to GPW population count. + Initiate exposures instance with value equal to GPW population count resampled + to 'res_arcsec' or 'target_grid'. The regridding process + is performed using 'util.coords.align_raster_data' with bilinear resampling + ('rasterio.warp.Resampling.bilinear'). Provide either `countries` or `shape`. Parameters @@ -406,8 +443,12 @@ def from_population( list containing country identifiers (name or iso3) shape : Shape, Polygon or MultiPolygon, optional geographical shape of target region, alternative to `countries`. - res_arcsec : int, optional - Resolution in arc seconds. The default is 30. + res_arcsec : float, optional + Horizontal resolution in arc-sec. + The default is 30 arcsec, this corresponds to roughly 1 km. + - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **nightlight grid**, set 'res_arcsec=15'. + if 'target_grid' is defined, this argument is not used. reference_year : int, optional Reference year (closest available GPW data year is used) The default is CONFIG.exposures.def_ref_year. @@ -416,6 +457,17 @@ def from_population( data_dir : Path, optional data directory. The default is None. Either countries or shape is required. + target_grid : dict or None, optional + A predefined grid onto which the exposure data should be mapped. + Default is None, the function will automatically define a target grid + based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' + is 15, the grid will be aligned to the nightlight data, otherwise + it will align to the GPW population grid. + The dictionary should include: + - "crs" : Coordinate reference system (e.g., "EPSG:4326") + - "transform" : Affine transformation defining pixel locations + - "width" : Number of columns in the grid + - "height" : Number of rows in the grid Raises ------ @@ -489,6 +541,11 @@ def from_shape_and_countries( ): """ create LitPop exposure for `country` and then crop to given shape. + 'Lit' and 'pop' data are first resampled to the same grid before being + used to estimate the exposure value. The data are regridded based on the + 'res_arcsec' argument or a predefined 'target_grid'. The regridding process + is performed using 'util.coords.align_raster_data' with bilinear resampling + ('rasterio.warp.Resampling.bilinear'). Parameters ---------- @@ -501,6 +558,9 @@ def from_shape_and_countries( res_arcsec : float, optional Horizontal resolution in arc-sec. The default is 30 arcsec, this corresponds to roughly 1 km. + - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **nightlight grid**, set 'res_arcsec=15'. + if 'target_grid' is defined, this argument is not used. exponents : tuple of two integers, optional Defining power with which lit (nightlights) and pop (gpw) go into LitPop. Default: (1, 1) @@ -531,6 +591,18 @@ def from_shape_and_countries( The default is GPW_VERSION data_dir : Path, optional redefines path to input data directory. The default is SYSTEM_DIR. + target_grid : dict or None, optional + A predefined grid onto which the exposure data should be mapped. + Default is None, the function will automatically define a target grid + based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' + is 15, the grid will be aligned to the nightlight data, otherwise + it will align to the GPW population grid. + The dictionary should include: + - "crs" : Coordinate reference system (e.g., "EPSG:4326") + - "transform" : Affine transformation defining pixel locations + - "resolution" : Spatial resolution in degrees + - "width" : Number of columns in the grid + - "height" : Number of rows in the grid Raises ------ @@ -555,13 +627,13 @@ def from_shape_and_countries( if isinstance(shape, Shape): # get gdf with geometries of points within shape: - shape_gdf = _get_litpop_single_polygon( + shape_gdf, _ = _get_litpop_single_polygon( shape, reference_year, data_dir, gpw_version, exponents, - target_grid=target_grid + target_grid ) shape_gdf = shape_gdf.drop( columns=shape_gdf.columns[shape_gdf.columns != "geometry"] @@ -653,6 +725,11 @@ def from_shape( and total value need to be provided manually. If these required input parameters are not known / available, better initiate Exposure for entire country and extract shape afterwards. + 'Lit' and 'pop' data are first resampled to the same grid before being + used to estimate the exposure value. The data are regridded based on the + 'res_arcsec' argument or a predefined 'target_grid'. The regridding process + is performed using 'util.coords.align_raster_data' with bilinear resampling + ('rasterio.warp.Resampling.bilinear'). Parameters ---------- @@ -663,7 +740,10 @@ def from_shape( If None, no value is disaggregated. res_arcsec : float, optional Horizontal resolution in arc-sec. - The default 30 arcsec corresponds to roughly 1 km. + The default is 30 arcsec, this corresponds to roughly 1 km. + - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **nightlight grid**, set 'res_arcsec=15'. + if 'target_grid' is defined, this argument is not used. exponents : tuple of two integers, optional Defining power with which lit (nightlights) and pop (gpw) go into LitPop. value_unit : str @@ -681,6 +761,17 @@ def from_shape( The default is set in CONFIG. data_dir : Path, optional redefines path to input data directory. The default is SYSTEM_DIR. + target_grid : dict or None, optional + A predefined grid onto which the exposure data should be mapped. + Default is None, the function will automatically define a target grid + based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' + is 15, the grid will be aligned to the nightlight data, otherwise + it will align to the GPW population grid. + The dictionary should include: + - "crs" : Coordinate reference system (e.g., "EPSG:4326") + - "transform" : Affine transformation defining pixel locations + - "width" : Number of columns in the grid + - "height" : Number of rows in the grid Raises ------ @@ -705,8 +796,8 @@ def from_shape( data_dir, gpw_version, exponents, - region_id, - target_grid=target_grid + target_grid, + region_id ) # disaggregate total value proportional to LitPop values: @@ -771,8 +862,6 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): Parameters ---------- - country_geometry : shapely.geometry (Polygon or MultiPolygon) - The country geometry used to extract relevant grid data. reference_year : int The reference year for population and nightlight data. gpw_version : int @@ -785,35 +874,38 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): Returns ------- target_grid : dict - A dictionary containing metadata for the global target grid. + A dictionary containing grid metadata, following the raster grid + specification. """ - res_deg = res_arcsec / 3600 + res_deg = res_arcsec / 3600 if res_arcsec == 15: # Black Marble Nightlights Grid global_crs = rasterio.crs.CRS.from_epsg(4326) # WGS84 projection global_transform = Affine( res_deg, 0, -180, - 0, -res_deg, 90 - ) + 0, -res_deg, 90) global_width = round(360 / res_deg) global_height = round(180 / res_deg) else: # GPW Population Grid - file_path = pop_util.get_gpw_file_path(gpw_version, - reference_year, - data_dir=data_dir, + file_path = pop_util.get_gpw_file_path(gpw_version, + reference_year, + data_dir=data_dir, verbose=False) with rasterio.open(file_path, "r") as src: - global_crs = src.crs - gpw_transform = src.transform + global_crs = src.crs + gpw_transform = src.transform # Align grid resolution with GPW dataset aligned_lon_min = -180 + (round((gpw_transform[2] - (-180)) / res_deg) * res_deg) aligned_lat_max = 90 - (round((90 - gpw_transform[5]) / res_deg) * res_deg) + global_transform = Affine( - res_deg, 0, aligned_lon_min, + res_deg, 0, aligned_lon_min, 0, -res_deg, aligned_lat_max) + global_width = round(360 / res_deg) global_height = round(180 / res_deg) + # Define the target grid using the computed values target_grid = { "driver": "GTiff", @@ -829,15 +921,13 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): @staticmethod def _from_country( country, - res_arcsec=30, + target_grid, exponents=(1, 1), fin_mode=None, total_value=None, reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, - data_dir=SYSTEM_DIR, - target_grid=None - ): + data_dir=SYSTEM_DIR): """init LitPop exposure object for one single country See docstring of from_countries() for detailled description of parameters. @@ -845,8 +935,9 @@ def _from_country( ---------- country : str or int country identifier as iso3alpha, iso3num or name. - res_arcsec : float, optional - horizontal resolution in arc-sec. + target_grid : dict + A dictionary containing grid metadata, following the raster grid + specification. exponents : tuple of two integers, optional fin_mode : str, optional total_value : numeric, optional @@ -854,7 +945,7 @@ def _from_country( gpw_version : int, optional data_dir : Path, optional redefines path to input data directory. The default is SYSTEM_DIR. - + Raises ------ ValueError @@ -893,9 +984,9 @@ def _from_country( data_dir, gpw_version, exponents, + target_grid, verbose=(idx > 0), region_id=iso3n, - target_grid=target_grid ) if gdf_tmp is None: LOGGER.debug( @@ -929,7 +1020,29 @@ def _from_country( set_country = set_countries def _crop_target_grid(target_grid, polygon): - """Crop the target grid to match the polygon's bounding box.""" + """ + Crop the target grid to match the bounding box of a given polygon. + + This function extracts a subset of the target grid that spatially + aligns with the bounding box of the given polygon. The cropped grid + retains the same resolution and coordinate reference system (CRS) as + the original grid but is limited to the extent of the polygon. + + Parameters + ---------- + target_grid : dict + A dictionary containing grid metadata, following the raster grid + specification. + polygon : shapely.geometry.Polygon + The polygon whose bounding box will be used to crop the grid. + + Returns + ------- + cropped_grid : dict + A dictionary with the same structure as 'target_grid', but with + updated 'width', 'height', and 'transform' to match the + polygon's bounding box. + """ bbox = polygon.bounds # (minx, miny, maxx, maxy) min_lon, min_lat, max_lon, max_lat = bbox @@ -955,7 +1068,6 @@ def _crop_target_grid(target_grid, polygon): target_grid["transform"].a, 0, cropped_min_lon, 0, target_grid["transform"].e, cropped_max_lat ) - return cropped_grid def _get_litpop_single_polygon( @@ -964,9 +1076,9 @@ def _get_litpop_single_polygon( data_dir, gpw_version, exponents, + target_grid, region_id=None, verbose=False, - target_grid=None ): """load nightlight (nl) and population (pop) data in rastered 2d arrays and apply rescaling (resolution reprojection) and LitPop core calculation, @@ -991,6 +1103,9 @@ def _get_litpop_single_polygon( exponents : tuple of two integers Defining power with which lit (nightlights) and pop (gpw) go into LitPop. To get nightlights^3 without population count: (3, 0). To use population count alone: (0, 1). + target_grid : dict + A dictionary containing grid metadata, following the raster grid + specification. region_id : int, optional if provided, region_id of gdf is set to value. The default is ``None``, in which case the region is determined automatically for every location of the resulting @@ -999,9 +1114,6 @@ def _get_litpop_single_polygon( Enable verbose logging about the used GPW version and reference year as well as about the boundary case where no grid points from the GPW grid are contained in the specified polygon. Default: False. - target_grid : dict, optional - Custom target grid metadata (CRS, transform, resolution, width, height). - If provided, this grid will be used instead of aligning to population/nightlight data. Returns ------- @@ -1213,7 +1325,8 @@ def reproject_input_data( The meta data with the reference grid used to define the global destination grid should be first in the list, e.g., GPW population data for LitPop. target_grid : dict - Target grid metadata. This defines the destination CRS, transform, and resolution. + A dictionary containing grid metadata, following the raster grid + specification. resampling : resampling function, optional The default is rasterio.warp.Resampling.bilinear conserve : str, optional, either 'mean' or 'sum' @@ -1238,14 +1351,8 @@ def reproject_input_data( # Compute global origins (top-left corner) global_origins = (dst_transform[2], dst_transform[5]) # (lon, lat) data_out_list = [None] * len(data_array_list) - meta_out = { - "dtype": target_grid["dtype"], - "nodata": target_grid["nodata"], - "crs": dst_crs, - "transform": dst_transform, - "width": dst_width, - "height": dst_height, - } + meta_out = target_grid.copy() # Copy all existing fields directly + for idx, _ in enumerate(data_array_list): # reproject data grid: @@ -1490,6 +1597,7 @@ def _calc_admin1_one_country( reference_year : int gpw_version: int data_dir : Path + target_grid: dict Returns ------- diff --git a/climada/entity/exposures/test/test_litpop.py b/climada/entity/exposures/test/test_litpop.py index 242eb95d35..d6c799fb91 100644 --- a/climada/entity/exposures/test/test_litpop.py +++ b/climada/entity/exposures/test/test_litpop.py @@ -20,7 +20,6 @@ """ import unittest -import copy import numpy as np from rasterio import Affine import rasterio diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index 2c2ddba88b..a7b2248180 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -23,6 +23,8 @@ import numpy as np from shapely.geometry import Polygon +from rasterio import Affine +from rasterio.crs import CRS import climada.util.coordinates as u_coord from climada import CONFIG @@ -411,9 +413,128 @@ def test_load_gpw_pop_shape_pass(self): self.assertIn("lease download", err.args[0]) self.skipTest("GPW input data for GPW v4.%i not found." % (gpw_version)) +class TestLitPopGridAlignment(unittest.TestCase): + """Test grid alignment, geometry validity, and expected structure in LitPop.""" + + @classmethod + def setUpClass(cls): + """Set up common test parameters.""" + # Define the target resolution in degrees (150 arcseconds = 0.0416667°) + target_res_deg = 150 / 3600 # Convert arcseconds to degrees (0.0416667°) + + #Ensure the grid is aligned exactly at half-cell positions + aligned_lon_min = -180 + (target_res_deg / 2) # Shift half a cell right + aligned_lat_max = 90 - (target_res_deg / 2) # Shift half a cell down + # Compute grid width & height + width = int(360 / target_res_deg) # 360° longitude range + height = int(180 / target_res_deg) # 180° latitude range + + # Define the affine transform for the target grid + transform = Affine( + target_res_deg, 0, aligned_lon_min, # X resolution, rotation, X origin + 0, -target_res_deg, aligned_lat_max # Y rotation, Y resolution (negative), Y origin + ) + + # Define the target grid metadata + cls.target_grid = { + "driver": "GTiff", + "width": width, + "height": height, + "count": 1, # Assuming a single-band raster + "crs": CRS.from_epsg(4326), # WGS84 (Lat/Lon) + "transform": transform, + } + + cls.test_countries = ["JPN", "TGO"] + cls.test_shape = Polygon([(-10, -10), (10, -10), (10, 10), (-10, 10)]) + cls.test_total_value = 1e6 + cls.test_year = 2020 + + # Define test cases with correct resolution + cls.test_cases = [ + ("from_countries", {"countries": cls.test_countries, "res_arcsec": 30}), + ("from_nightlight_intensity", {"countries": cls.test_countries, "res_arcsec": 30}), + ("from_population", {"countries": cls.test_countries, "res_arcsec": 30}), + ("from_shape_and_countries", {"shape": cls.test_shape, "countries": cls.test_countries, "res_arcsec": 30}), + ("from_shape", {"shape": cls.test_shape, "total_value": cls.test_total_value, "res_arcsec": 30}), + ] + + def check_exposure_geometry(self, exp, method_name): + """Validate that the exposure data points align correctly to the expected grid.""" + + # Ensure exposure data exists + self.assertIsNotNone(exp, f"{method_name} returned None") + self.assertFalse(exp.gdf.empty, f"GeoDataFrame is empty in {method_name}!") + + # Extract coordinates from exposure dataset + latitudes = exp.gdf.geometry.y + longitudes = exp.gdf.geometry.x + + # Define the target resolution in degrees (150 arcseconds = 0.0416667°) + target_res_deg = 150 / 3600 # Convert arcseconds to degrees (0.0416667°) + + #Ensure the grid is aligned exactly at half-cell positions + aligned_lon_min = -180 # Shift half a cell right + aligned_lat_max = 90 + # Compute expected alignment positions + expected_lons = np.arange(aligned_lon_min, aligned_lon_min+360, target_res_deg) + expected_lats = np.arange(aligned_lat_max-180, aligned_lat_max , target_res_deg ) + # Check that all longitudes and latitudes align to expected values + self.assertTrue( + np.all(np.isin(np.round(longitudes, 6), np.round(expected_lons, 6))), + f"Longitude coordinates are not aligned to the expected grid in {method_name}!" + ) + self.assertTrue( + np.all(np.isin(np.round(latitudes, 6), np.round(expected_lats, 6))), + f"Latitude coordinates are not aligned to the expected grid in {method_name}!" + ) + + def test_litpop_grid_alignment(self): + """Loop through LitPop methods and validate grid alignment.""" + for method_name, kwargs in self.test_cases: + with self.subTest(method=method_name): + exp = getattr(lp.LitPop, method_name)(**kwargs, target_grid=self.target_grid) + self.check_exposure_geometry(exp, method_name) + + def test_litpop_grid_consistency(self): + """Ensure total exposure values and dataset lengths remain consistent across different grid configurations.""" + + country = "FRA" # Example country, adjust as needed + default_res = 30 # Default resolution in arcseconds + test_res = 150 # Test resolution in arcseconds (5x coarser) + + # Compute exposure with different configurations + exp_default = lp.LitPop.from_countries([country]) # Default settings (30 arcseconds) + exp_res = lp.LitPop.from_countries([country], res_arcsec=test_res) # Using specified resolution (150 arcseconds) + exp_grid = lp.LitPop.from_countries([country], target_grid=self.target_grid) # Using target grid + + # Extract total exposure values + sum_default = exp_default.gdf["value"].sum() + sum_res = exp_res.gdf["value"].sum() + sum_grid = exp_grid.gdf["value"].sum() + + # Extract dataset lengths (number of grid points) + len_default = len(exp_default.gdf) + len_res = len(exp_res.gdf) + len_grid = len(exp_grid.gdf) + + # Allow for a 10% difference in total exposure values + tolerance = 0.1 * sum_default # 10% of default value + + # Check that the total exposure values are within tolerance + self.assertAlmostEqual(sum_res, sum_default, delta=tolerance, + msg="Exposure sum differs by more than 10% when using resolution.") + self.assertAlmostEqual(sum_grid, sum_default, delta=tolerance, + msg="Exposure sum differs by more than 10% when using a target grid.") + + # Check that the number of grid points is similar when using the same resolution + self.assertAlmostEqual(len_res, len_grid, delta=0.05 * len_res, + msg="Dataset lengths differ significantly between resolution and target grid.") -# Execute Tests if __name__ == "__main__": - TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopExposure) - TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAdmin1)) + # TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopExposure) + # TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAdmin1)) + # unittest.TextTestRunner(verbosity=2).run(TESTS) + + TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopGridAlignment) unittest.TextTestRunner(verbosity=2).run(TESTS) From 8f5ba60e06c52f5309903ffad53383823cc794ea Mon Sep 17 00:00:00 2001 From: zeliest Date: Wed, 26 Mar 2025 11:50:36 +0100 Subject: [PATCH 08/18] Uncomment test litpop --- climada/test/test_litpop_integr.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index a7b2248180..7aebd20c27 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -532,9 +532,7 @@ def test_litpop_grid_consistency(self): msg="Dataset lengths differ significantly between resolution and target grid.") if __name__ == "__main__": - # TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopExposure) - # TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAdmin1)) - # unittest.TextTestRunner(verbosity=2).run(TESTS) - + TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopExposure) + TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAdmin1)) TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopGridAlignment) unittest.TextTestRunner(verbosity=2).run(TESTS) From a9a2555ce355f95f705bfa36a5349b905db7e233 Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Thu, 8 May 2025 14:57:19 +0200 Subject: [PATCH 09/18] pylint --- climada/entity/exposures/litpop/litpop.py | 207 ++++++++++--------- climada/entity/exposures/test/test_litpop.py | 25 ++- climada/test/test_litpop_integr.py | 118 ++++++++--- 3 files changed, 204 insertions(+), 146 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index b15cb36a2d..811e28c2aa 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -26,8 +26,8 @@ import pandas as pd import rasterio import shapely -from shapefile import Shape from affine import Affine +from shapefile import Shape import climada.util.coordinates as u_coord import climada.util.finance as u_fin @@ -109,13 +109,13 @@ def from_countries( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, - target_grid=None + target_grid=None, ): - """Init new LitPop exposure object for a list of countries (admin 0). - 'Lit' and 'pop' data are first resampled to the same grid before being - used to estimate the exposure value. The data are regridded based on the - 'res_arcsec' argument or a predefined 'target_grid'. The regridding process - is performed using 'util.coords.align_raster_data' with bilinear resampling + """Init new LitPop exposure object for a list of countries (admin 0). + 'Lit' and 'pop' data are first resampled to the same grid before being + used to estimate the exposure value. The data are regridded based on the + 'res_arcsec' argument or a predefined 'target_grid'. The regridding process + is performed using 'util.coords.align_raster_data' with bilinear resampling ('rasterio.warp.Resampling.bilinear'). Sets attributes `ref_year`, `crs`, `value`, `geometry`, `meta`, @@ -129,7 +129,7 @@ def from_countries( res_arcsec : float, optional Horizontal resolution in arc-sec. The default is 30 arcsec, this corresponds to roughly 1 km. - - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **GPW population grid**, set 'res_arcsec=30'. - To use the original **nightlight grid**, set 'res_arcsec=15'. if 'target_grid' is defined, this argument is not used. exponents : tuple of two integers, optional @@ -170,10 +170,10 @@ def from_countries( redefines path to input data directory. The default is SYSTEM_DIR. target_grid : dict or None, optional A predefined grid onto which the exposure data should be mapped. - Default is None, the function will automatically define a target grid + Default is None, the function will automatically define a target grid based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' is 15, the grid will be aligned to the nightlight data, otherwise - it will align to the GPW population grid. + it will align to the GPW population grid. The dictionary should include: - "crs" : Coordinate reference system (e.g., "EPSG:4326") - "transform" : Affine transformation defining pixel locations @@ -196,13 +196,14 @@ def from_countries( total_values = [None] * len(countries) elif len(total_values) != len(countries): raise ValueError( - "'countries' and 'total_values' must be lists of same length") + "'countries' and 'total_values' must be lists of same length" + ) if target_grid is None: target_grid = cls._define_target_grid( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, - res_arcsec=res_arcsec + res_arcsec=res_arcsec, ) # litpop_list is initiated, a list containing one Exposure instance per # country and None for countries that could not be identified: @@ -220,7 +221,7 @@ def from_countries( reference_year, gpw_version, data_dir, - target_grid + target_grid, ) for tot_value, country in zip(total_values, countries) ] @@ -236,7 +237,8 @@ def from_countries( total_value=total_values[idc], reference_year=reference_year, gpw_version=gpw_version, - data_dir=data_dir, target_grid=target_grid + data_dir=data_dir, + target_grid=target_grid, ) for idc, country in enumerate(countries) ] @@ -248,7 +250,7 @@ def from_countries( country for lp, country in zip(litpop_list, countries) if lp is None ] if not countries_in: - raise ValueError("No valid country identified in %s, aborting." % countries) + raise ValueError(f"No valid country identified in {countries}, aborting.") litpop_list = [exp for exp in litpop_list if exp is not None] if countries_out: LOGGER.warning( @@ -316,15 +318,15 @@ def from_nightlight_intensity( res_arcsec=15, reference_year=DEF_REF_YEAR, data_dir=SYSTEM_DIR, - target_grid=None + target_grid=None, ): """ Wrapper around `from_countries` / `from_shape`. Initiate exposures instance with value equal to the original BlackMarble - nightlight intensity resampled based on the 'res_arcsec' argument or a predefined - 'target_grid'. The regridding process is performed using 'util.coords.align_raster_data' with bilinear resampling - ('rasterio.warp.Resampling.bilinear'). + nightlight intensity resampled based on the 'res_arcsec' argument or a predefined + 'target_grid'. The regridding process is performed using 'util.coords.align_raster_data' + with bilinear resampling ('rasterio.warp.Resampling.bilinear'). Provide either `countries` or `shape`. @@ -337,7 +339,7 @@ def from_nightlight_intensity( res_arcsec : float, optional Horizontal resolution in arc-sec. The default is 30 arcsec, this corresponds to roughly 1 km. - - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **GPW population grid**, set 'res_arcsec=30'. - To use the original **nightlight grid**, set 'res_arcsec=15'. if 'target_grid' is defined, this argument is not used. reference_year : int, optional @@ -346,10 +348,10 @@ def from_nightlight_intensity( data directory. The default is None. target_grid : dict or None, optional A predefined grid onto which the exposure data should be mapped. - Default is None, the function will automatically define a target grid + Default is None, the function will automatically define a target grid based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' is 15, the grid will be aligned to the nightlight data, otherwise - it will align to the GPW population grid. + it will align to the GPW population grid. The dictionary should include: - "crs" : Coordinate reference system (e.g., "EPSG:4326") - "transform" : Affine transformation defining pixel locations @@ -373,11 +375,11 @@ def from_nightlight_intensity( ) if target_grid is None: target_grid = cls._define_target_grid( - reference_year=reference_year, - gpw_version=GPW_VERSION, - data_dir=data_dir, - res_arcsec=res_arcsec - ) + reference_year=reference_year, + gpw_version=GPW_VERSION, + data_dir=data_dir, + res_arcsec=res_arcsec, + ) if countries is not None: exp = cls.from_countries( countries, @@ -387,7 +389,7 @@ def from_nightlight_intensity( reference_year=reference_year, gpw_version=GPW_VERSION, data_dir=data_dir, - target_grid=target_grid + target_grid=target_grid, ) else: exp = cls.from_shape( @@ -399,7 +401,7 @@ def from_nightlight_intensity( reference_year=reference_year, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, - target_grid=target_grid + target_grid=target_grid, ) LOGGER.warning( "Note: set_nightlight_intensity sets values to raw nightlight intensity, " @@ -426,14 +428,14 @@ def from_population( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, - target_grid=None + target_grid=None, ): """ Wrapper around `from_countries` / `from_shape`. Initiate exposures instance with value equal to GPW population count resampled - to 'res_arcsec' or 'target_grid'. The regridding process - is performed using 'util.coords.align_raster_data' with bilinear resampling + to 'res_arcsec' or 'target_grid'. The regridding process + is performed using 'util.coords.align_raster_data' with bilinear resampling ('rasterio.warp.Resampling.bilinear'). Provide either `countries` or `shape`. @@ -446,7 +448,7 @@ def from_population( res_arcsec : float, optional Horizontal resolution in arc-sec. The default is 30 arcsec, this corresponds to roughly 1 km. - - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **GPW population grid**, set 'res_arcsec=30'. - To use the original **nightlight grid**, set 'res_arcsec=15'. if 'target_grid' is defined, this argument is not used. reference_year : int, optional @@ -459,10 +461,10 @@ def from_population( Either countries or shape is required. target_grid : dict or None, optional A predefined grid onto which the exposure data should be mapped. - Default is None, the function will automatically define a target grid + Default is None, the function will automatically define a target grid based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' is 15, the grid will be aligned to the nightlight data, otherwise - it will align to the GPW population grid. + it will align to the GPW population grid. The dictionary should include: - "crs" : Coordinate reference system (e.g., "EPSG:4326") - "transform" : Affine transformation defining pixel locations @@ -487,11 +489,11 @@ def from_population( ) if target_grid is None: target_grid = cls._define_target_grid( - reference_year=reference_year, - gpw_version=gpw_version, - data_dir=data_dir, - res_arcsec=res_arcsec - ) + reference_year=reference_year, + gpw_version=gpw_version, + data_dir=data_dir, + res_arcsec=res_arcsec, + ) if countries is not None: exp = cls.from_countries( countries, @@ -501,7 +503,7 @@ def from_population( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, - target_grid=target_grid + target_grid=target_grid, ) else: exp = cls.from_shape( @@ -513,7 +515,7 @@ def from_population( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, - target_grid=target_grid + target_grid=target_grid, ) return exp @@ -537,14 +539,14 @@ def from_shape_and_countries( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, - target_grid=None + target_grid=None, ): """ create LitPop exposure for `country` and then crop to given shape. - 'Lit' and 'pop' data are first resampled to the same grid before being - used to estimate the exposure value. The data are regridded based on the - 'res_arcsec' argument or a predefined 'target_grid'. The regridding process - is performed using 'util.coords.align_raster_data' with bilinear resampling + 'Lit' and 'pop' data are first resampled to the same grid before being + used to estimate the exposure value. The data are regridded based on the + 'res_arcsec' argument or a predefined 'target_grid'. The regridding process + is performed using 'util.coords.align_raster_data' with bilinear resampling ('rasterio.warp.Resampling.bilinear'). Parameters @@ -558,7 +560,7 @@ def from_shape_and_countries( res_arcsec : float, optional Horizontal resolution in arc-sec. The default is 30 arcsec, this corresponds to roughly 1 km. - - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **GPW population grid**, set 'res_arcsec=30'. - To use the original **nightlight grid**, set 'res_arcsec=15'. if 'target_grid' is defined, this argument is not used. exponents : tuple of two integers, optional @@ -593,10 +595,10 @@ def from_shape_and_countries( redefines path to input data directory. The default is SYSTEM_DIR. target_grid : dict or None, optional A predefined grid onto which the exposure data should be mapped. - Default is None, the function will automatically define a target grid + Default is None, the function will automatically define a target grid based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' is 15, the grid will be aligned to the nightlight data, otherwise - it will align to the GPW population grid. + it will align to the GPW population grid. The dictionary should include: - "crs" : Coordinate reference system (e.g., "EPSG:4326") - "transform" : Affine transformation defining pixel locations @@ -622,18 +624,13 @@ def from_shape_and_countries( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, - target_grid=target_grid + target_grid=target_grid, ) if isinstance(shape, Shape): # get gdf with geometries of points within shape: shape_gdf, _ = _get_litpop_single_polygon( - shape, - reference_year, - data_dir, - gpw_version, - exponents, - target_grid + shape, reference_year, data_dir, gpw_version, exponents, target_grid ) shape_gdf = shape_gdf.drop( columns=shape_gdf.columns[shape_gdf.columns != "geometry"] @@ -712,7 +709,7 @@ def from_shape( reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, data_dir=SYSTEM_DIR, - target_grid=None + target_grid=None, ): """init LitPop exposure object for a custom shape. Requires user input regarding the total value to be disaggregated. @@ -725,10 +722,10 @@ def from_shape( and total value need to be provided manually. If these required input parameters are not known / available, better initiate Exposure for entire country and extract shape afterwards. - 'Lit' and 'pop' data are first resampled to the same grid before being - used to estimate the exposure value. The data are regridded based on the - 'res_arcsec' argument or a predefined 'target_grid'. The regridding process - is performed using 'util.coords.align_raster_data' with bilinear resampling + 'Lit' and 'pop' data are first resampled to the same grid before being + used to estimate the exposure value. The data are regridded based on the + 'res_arcsec' argument or a predefined 'target_grid'. The regridding process + is performed using 'util.coords.align_raster_data' with bilinear resampling ('rasterio.warp.Resampling.bilinear'). Parameters @@ -741,7 +738,7 @@ def from_shape( res_arcsec : float, optional Horizontal resolution in arc-sec. The default is 30 arcsec, this corresponds to roughly 1 km. - - To use the original **GPW population grid**, set 'res_arcsec=30'. + - To use the original **GPW population grid**, set 'res_arcsec=30'. - To use the original **nightlight grid**, set 'res_arcsec=15'. if 'target_grid' is defined, this argument is not used. exponents : tuple of two integers, optional @@ -763,10 +760,10 @@ def from_shape( redefines path to input data directory. The default is SYSTEM_DIR. target_grid : dict or None, optional A predefined grid onto which the exposure data should be mapped. - Default is None, the function will automatically define a target grid + Default is None, the function will automatically define a target grid based on the 'res_arcsec' parameters. Specifically, if 'res_arcsec' is 15, the grid will be aligned to the nightlight data, otherwise - it will align to the GPW population grid. + it will align to the GPW population grid. The dictionary should include: - "crs" : Coordinate reference system (e.g., "EPSG:4326") - "transform" : Affine transformation defining pixel locations @@ -797,7 +794,7 @@ def from_shape( gpw_version, exponents, target_grid, - region_id + region_id, ) # disaggregate total value proportional to LitPop values: @@ -874,35 +871,35 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): Returns ------- target_grid : dict - A dictionary containing grid metadata, following the raster grid - specification. + A dictionary containing grid metadata, following the raster grid + specification. """ res_deg = res_arcsec / 3600 if res_arcsec == 15: # Black Marble Nightlights Grid + # pylint: disable=c-extension-no-member global_crs = rasterio.crs.CRS.from_epsg(4326) # WGS84 projection - global_transform = Affine( - res_deg, 0, -180, - 0, -res_deg, 90) + global_transform = Affine(res_deg, 0, -180, 0, -res_deg, 90) global_width = round(360 / res_deg) global_height = round(180 / res_deg) else: # GPW Population Grid - file_path = pop_util.get_gpw_file_path(gpw_version, - reference_year, - data_dir=data_dir, - verbose=False) + file_path = pop_util.get_gpw_file_path( + gpw_version, reference_year, data_dir=data_dir, verbose=False + ) with rasterio.open(file_path, "r") as src: global_crs = src.crs gpw_transform = src.transform # Align grid resolution with GPW dataset - aligned_lon_min = -180 + (round((gpw_transform[2] - (-180)) / res_deg) * res_deg) + aligned_lon_min = -180 + ( + round((gpw_transform[2] - (-180)) / res_deg) * res_deg + ) aligned_lat_max = 90 - (round((90 - gpw_transform[5]) / res_deg) * res_deg) - + global_transform = Affine( - res_deg, 0, aligned_lon_min, - 0, -res_deg, aligned_lat_max) - + res_deg, 0, aligned_lon_min, 0, -res_deg, aligned_lat_max + ) + global_width = round(360 / res_deg) global_height = round(180 / res_deg) @@ -917,7 +914,7 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): "transform": global_transform, } return target_grid - + @staticmethod def _from_country( country, @@ -927,7 +924,8 @@ def _from_country( total_value=None, reference_year=DEF_REF_YEAR, gpw_version=GPW_VERSION, - data_dir=SYSTEM_DIR): + data_dir=SYSTEM_DIR, + ): """init LitPop exposure object for one single country See docstring of from_countries() for detailled description of parameters. @@ -936,8 +934,8 @@ def _from_country( country : str or int country identifier as iso3alpha, iso3num or name. target_grid : dict - A dictionary containing grid metadata, following the raster grid - specification. + A dictionary containing grid metadata, following the raster grid + specification. exponents : tuple of two integers, optional fin_mode : str, optional total_value : numeric, optional @@ -945,7 +943,7 @@ def _from_country( gpw_version : int, optional data_dir : Path, optional redefines path to input data directory. The default is SYSTEM_DIR. - + Raises ------ ValueError @@ -990,7 +988,7 @@ def _from_country( ) if gdf_tmp is None: LOGGER.debug( - f"Skipping polygon with index {idx} for" + f" country {iso3a}." + "Skipping polygon with index %d for country %s.", idx, iso3a ) continue total_population += meta_tmp["total_population"] @@ -1019,28 +1017,29 @@ def _from_country( # Alias method names for backward compatibility: set_country = set_countries + def _crop_target_grid(target_grid, polygon): """ Crop the target grid to match the bounding box of a given polygon. - This function extracts a subset of the target grid that spatially - aligns with the bounding box of the given polygon. The cropped grid - retains the same resolution and coordinate reference system (CRS) as + This function extracts a subset of the target grid that spatially + aligns with the bounding box of the given polygon. The cropped grid + retains the same resolution and coordinate reference system (CRS) as the original grid but is limited to the extent of the polygon. Parameters ---------- target_grid : dict - A dictionary containing grid metadata, following the raster grid - specification. + A dictionary containing grid metadata, following the raster grid + specification. polygon : shapely.geometry.Polygon The polygon whose bounding box will be used to crop the grid. Returns ------- cropped_grid : dict - A dictionary with the same structure as 'target_grid', but with - updated 'width', 'height', and 'transform' to match the + A dictionary with the same structure as 'target_grid', but with + updated 'width', 'height', and 'transform' to match the polygon's bounding box. """ bbox = polygon.bounds # (minx, miny, maxx, maxy) @@ -1065,11 +1064,16 @@ def _crop_target_grid(target_grid, polygon): cropped_grid["width"] = col_max - col_min cropped_grid["height"] = row_max - row_min cropped_grid["transform"] = Affine( - target_grid["transform"].a, 0, cropped_min_lon, - 0, target_grid["transform"].e, cropped_max_lat + target_grid["transform"].a, + 0, + cropped_min_lon, + 0, + target_grid["transform"].e, + cropped_max_lat, ) return cropped_grid + def _get_litpop_single_polygon( polygon, reference_year, @@ -1104,8 +1108,8 @@ def _get_litpop_single_polygon( Defining power with which lit (nightlights) and pop (gpw) go into LitPop. To get nightlights^3 without population count: (3, 0). To use population count alone: (0, 1). target_grid : dict - A dictionary containing grid metadata, following the raster grid - specification. + A dictionary containing grid metadata, following the raster grid + specification. region_id : int, optional if provided, region_id of gdf is set to value. The default is ``None``, in which case the region is determined automatically for every location of the resulting @@ -1325,8 +1329,8 @@ def reproject_input_data( The meta data with the reference grid used to define the global destination grid should be first in the list, e.g., GPW population data for LitPop. target_grid : dict - A dictionary containing grid metadata, following the raster grid - specification. + A dictionary containing grid metadata, following the raster grid + specification. resampling : resampling function, optional The default is rasterio.warp.Resampling.bilinear conserve : str, optional, either 'mean' or 'sum' @@ -1353,7 +1357,6 @@ def reproject_input_data( data_out_list = [None] * len(data_array_list) meta_out = target_grid.copy() # Copy all existing fields directly - for idx, _ in enumerate(data_array_list): # reproject data grid: dst_bounds = rasterio.transform.array_bounds( @@ -1574,7 +1577,7 @@ def _calc_admin1_one_country( reference_year, gpw_version, data_dir, - target_grid + target_grid, ): """ Calculates the LitPop on admin1 level for provinces/states where such information are @@ -1652,7 +1655,7 @@ def _calc_admin1_one_country( reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, - target_grid=target_grid + target_grid=target_grid, ) ) exp_list[-1].gdf["admin1"] = record["name"] diff --git a/climada/entity/exposures/test/test_litpop.py b/climada/entity/exposures/test/test_litpop.py index d6c799fb91..76f7d166a9 100644 --- a/climada/entity/exposures/test/test_litpop.py +++ b/climada/entity/exposures/test/test_litpop.py @@ -20,10 +20,12 @@ """ import unittest + import numpy as np -from rasterio import Affine import rasterio +from rasterio import Affine from rasterio.crs import CRS + from climada.entity.exposures.litpop import litpop as lp @@ -113,6 +115,7 @@ def data_arrays_resampling_demo(): ] return data_arrays, meta_list + class TestLitPop(unittest.TestCase): """Test LitPop Class methods and functions""" @@ -122,9 +125,7 @@ def test_reproject_input_data_downsample(self): data_in, meta_list = data_arrays_resampling_demo() # data_out, meta_out = lp.reproject_input_data( - data_in, - meta_list, - target_grid=meta_list[0] + data_in, meta_list, target_grid=meta_list[0] ) # test reference data unchanged: np.testing.assert_array_equal(data_in[0], data_out[0]) @@ -143,7 +144,7 @@ def test_reproject_input_data_downsample_conserve_sum(self): data_out, meta_out = lp.reproject_input_data( data_in, meta_list, - target_grid = meta_list[0], + target_grid=meta_list[0], conserve="sum", ) np.testing.assert_array_equal(data_in[0], data_out[0]) @@ -155,15 +156,15 @@ def test_target_grid_alignment(self): orig_res = 1.0 orig_transform = Affine(orig_res, 0, -10.0, 0, -orig_res, 40.0) - + data_in, meta_list = data_arrays_resampling_demo() meta_list[0]["transform"] = orig_transform meta_list[0]["width"] = 3 meta_list[0]["height"] = 3 - target_res = 2.0 + target_res = 2.0 target_transform = Affine(target_res, 0, -10.0, 0, -target_res, 40.0) - + target_grid = { "driver": "GTiff", "dtype": "float32", @@ -179,10 +180,14 @@ def test_target_grid_alignment(self): data_in, meta_list, target_grid=target_grid, - resampling=rasterio.warp.Resampling.bilinear + resampling=rasterio.warp.Resampling.bilinear, ) - self.assertEqual(meta_out["transform"], target_transform, "Transform does not match target grid!") + self.assertEqual( + meta_out["transform"], + target_transform, + "Transform does not match target grid!", + ) self.assertEqual(meta_out["width"], target_grid["width"], "Width mismatch!") self.assertEqual(meta_out["height"], target_grid["height"], "Height mismatch!") diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index ee00293246..5fe6e764bf 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -22,9 +22,9 @@ import unittest import numpy as np -from shapely.geometry import Polygon from rasterio import Affine from rasterio.crs import CRS +from shapely.geometry import Polygon import climada.util.coordinates as u_coord from climada import CONFIG @@ -117,7 +117,6 @@ def test_switzerland30normPop_pass(self): fin_mode=fin_mode, reference_year=2015, ) - # print(cm) self.assertIn("LitPop: Init Exposure for country: CHE", cm.output[0]) self.assertEqual(ent.region_id.min(), 756) self.assertEqual(ent.region_id.max(), 756) @@ -359,10 +358,13 @@ def test_brandenburg(self): admin1_info = admin1_info[country] admin1_shapes = admin1_shapes[country] admin1_names = [record["name"] for record in admin1_info] - print(admin1_names) - for idx, name in enumerate(admin1_names): - if admin1_names[idx] == state_name: - break + + idx = next( + idx + for idx, _name in enumerate(admin1_names) + if admin1_names[idx] == state_name + ) + self.assertEqual(idx, 8) # init LitPop for Brandenburg exp_bra2 = lp.LitPop.from_shape_and_countries( @@ -394,6 +396,7 @@ def test_get_gpw_file_path_pass(self): self.assertIn("gpw_v4_population", str(path)) except FileExistsError as err: self.assertIn("lease download", err.args[0]) + # pylint: disabled=consider-using-f-string self.skipTest("GPW input data for GPW v4.%i not found." % (gpw_version)) def test_load_gpw_pop_shape_pass(self): @@ -413,8 +416,10 @@ def test_load_gpw_pop_shape_pass(self): self.assertEqual(len(data.shape), 2) except FileExistsError as err: self.assertIn("lease download", err.args[0]) + # pylint: disabled=consider-using-f-string self.skipTest("GPW input data for GPW v4.%i not found." % (gpw_version)) + class TestLitPopGridAlignment(unittest.TestCase): """Test grid alignment, geometry validity, and expected structure in LitPop.""" @@ -424,17 +429,22 @@ def setUpClass(cls): # Define the target resolution in degrees (150 arcseconds = 0.0416667°) target_res_deg = 150 / 3600 # Convert arcseconds to degrees (0.0416667°) - #Ensure the grid is aligned exactly at half-cell positions + # Ensure the grid is aligned exactly at half-cell positions aligned_lon_min = -180 + (target_res_deg / 2) # Shift half a cell right - aligned_lat_max = 90 - (target_res_deg / 2) # Shift half a cell down + aligned_lat_max = 90 - (target_res_deg / 2) # Shift half a cell down # Compute grid width & height width = int(360 / target_res_deg) # 360° longitude range height = int(180 / target_res_deg) # 180° latitude range # Define the affine transform for the target grid transform = Affine( - target_res_deg, 0, aligned_lon_min, # X resolution, rotation, X origin - 0, -target_res_deg, aligned_lat_max # Y rotation, Y resolution (negative), Y origin + target_res_deg, + 0, + aligned_lon_min, # X resolution, rotation, X origin + 0, + -target_res_deg, + aligned_lat_max, # Y rotation, Y resolution (negative), + # Y origin ) # Define the target grid metadata @@ -455,15 +465,32 @@ def setUpClass(cls): # Define test cases with correct resolution cls.test_cases = [ ("from_countries", {"countries": cls.test_countries, "res_arcsec": 30}), - ("from_nightlight_intensity", {"countries": cls.test_countries, "res_arcsec": 30}), + ( + "from_nightlight_intensity", + {"countries": cls.test_countries, "res_arcsec": 30}, + ), ("from_population", {"countries": cls.test_countries, "res_arcsec": 30}), - ("from_shape_and_countries", {"shape": cls.test_shape, "countries": cls.test_countries, "res_arcsec": 30}), - ("from_shape", {"shape": cls.test_shape, "total_value": cls.test_total_value, "res_arcsec": 30}), + ( + "from_shape_and_countries", + { + "shape": cls.test_shape, + "countries": cls.test_countries, + "res_arcsec": 30, + }, + ), + ( + "from_shape", + { + "shape": cls.test_shape, + "total_value": cls.test_total_value, + "res_arcsec": 30, + }, + ), ] def check_exposure_geometry(self, exp, method_name): """Validate that the exposure data points align correctly to the expected grid.""" - + # Ensure exposure data exists self.assertIsNotNone(exp, f"{method_name} returned None") self.assertFalse(exp.gdf.empty, f"GeoDataFrame is empty in {method_name}!") @@ -475,40 +502,50 @@ def check_exposure_geometry(self, exp, method_name): # Define the target resolution in degrees (150 arcseconds = 0.0416667°) target_res_deg = 150 / 3600 # Convert arcseconds to degrees (0.0416667°) - #Ensure the grid is aligned exactly at half-cell positions - aligned_lon_min = -180 # Shift half a cell right - aligned_lat_max = 90 + # Ensure the grid is aligned exactly at half-cell positions + aligned_lon_min = -180 # Shift half a cell right + aligned_lat_max = 90 # Compute expected alignment positions - expected_lons = np.arange(aligned_lon_min, aligned_lon_min+360, target_res_deg) - expected_lats = np.arange(aligned_lat_max-180, aligned_lat_max , target_res_deg ) - # Check that all longitudes and latitudes align to expected values + expected_lons = np.arange( + aligned_lon_min, aligned_lon_min + 360, target_res_deg + ) + expected_lats = np.arange( + aligned_lat_max - 180, aligned_lat_max, target_res_deg + ) + # Check that all longitudes and latitudes align to expected values self.assertTrue( np.all(np.isin(np.round(longitudes, 6), np.round(expected_lons, 6))), - f"Longitude coordinates are not aligned to the expected grid in {method_name}!" + f"Longitude coordinates are not aligned to the expected grid in {method_name}!", ) self.assertTrue( np.all(np.isin(np.round(latitudes, 6), np.round(expected_lats, 6))), - f"Latitude coordinates are not aligned to the expected grid in {method_name}!" + f"Latitude coordinates are not aligned to the expected grid in {method_name}!", ) def test_litpop_grid_alignment(self): """Loop through LitPop methods and validate grid alignment.""" for method_name, kwargs in self.test_cases: with self.subTest(method=method_name): - exp = getattr(lp.LitPop, method_name)(**kwargs, target_grid=self.target_grid) + exp = getattr(lp.LitPop, method_name)( + **kwargs, target_grid=self.target_grid + ) self.check_exposure_geometry(exp, method_name) def test_litpop_grid_consistency(self): - """Ensure total exposure values and dataset lengths remain consistent across different grid configurations.""" - + """Ensure total exposure values and dataset lengths remain consistent across different grid + configurations.""" + country = "FRA" # Example country, adjust as needed default_res = 30 # Default resolution in arcseconds test_res = 150 # Test resolution in arcseconds (5x coarser) # Compute exposure with different configurations - exp_default = lp.LitPop.from_countries([country]) # Default settings (30 arcseconds) - exp_res = lp.LitPop.from_countries([country], res_arcsec=test_res) # Using specified resolution (150 arcseconds) - exp_grid = lp.LitPop.from_countries([country], target_grid=self.target_grid) # Using target grid + # Default settings (30 arcseconds): + exp_default = lp.LitPop.from_countries([country]) + # Using specified resolution (150 arcseconds): + exp_res = lp.LitPop.from_countries([country], res_arcsec=test_res) + # Using target grid: + exp_grid = lp.LitPop.from_countries([country], target_grid=self.target_grid) # Extract total exposure values sum_default = exp_default.gdf["value"].sum() @@ -524,14 +561,27 @@ def test_litpop_grid_consistency(self): tolerance = 0.1 * sum_default # 10% of default value # Check that the total exposure values are within tolerance - self.assertAlmostEqual(sum_res, sum_default, delta=tolerance, - msg="Exposure sum differs by more than 10% when using resolution.") - self.assertAlmostEqual(sum_grid, sum_default, delta=tolerance, - msg="Exposure sum differs by more than 10% when using a target grid.") + self.assertAlmostEqual( + sum_res, + sum_default, + delta=tolerance, + msg="Exposure sum differs by more than 10% when using resolution.", + ) + self.assertAlmostEqual( + sum_grid, + sum_default, + delta=tolerance, + msg="Exposure sum differs by more than 10% when using a target grid.", + ) # Check that the number of grid points is similar when using the same resolution - self.assertAlmostEqual(len_res, len_grid, delta=0.05 * len_res, - msg="Dataset lengths differ significantly between resolution and target grid.") + self.assertAlmostEqual( + len_res, + len_grid, + delta=0.05 * len_res, + msg="Dataset lengths differ significantly between resolution and target grid.", + ) + if __name__ == "__main__": TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopExposure) From 7491c0efcc04a3740d99792572f8bee2afff35e2 Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Thu, 8 May 2025 15:21:24 +0200 Subject: [PATCH 10/18] pylint --- climada/entity/exposures/test/test_litpop.py | 30 ++++++++++++-------- climada/test/test_litpop_integr.py | 4 +-- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/climada/entity/exposures/test/test_litpop.py b/climada/entity/exposures/test/test_litpop.py index 76f7d166a9..2e0fd0dca8 100644 --- a/climada/entity/exposures/test/test_litpop.py +++ b/climada/entity/exposures/test/test_litpop.py @@ -32,22 +32,28 @@ def data_arrays_demo(number_of_arrays=2): """init demo data arrays (2d) for LitPop core calculations""" data_arrays = list() + # fmt: off if number_of_arrays > 0: - data_arrays.append(np.array([[0, 1, 2], [3, 4, 5]])) - # array([[0, 1, 2], - # [3, 4, 5]]) + data_arrays.append(np.array([ + [0, 1, 2], + [3, 4, 5] + ])) if number_of_arrays > 1: - data_arrays.append(np.array([[10, 10, 10], [1, 1, 1]])) - # array([[10, 10, 10], - # [1, 1, 1]]) + data_arrays.append(np.array([ + [10, 10, 10], + [ 1, 1, 1] + ])) if number_of_arrays > 2: - data_arrays.append(np.array([[0, 1, 10], [0, 1, 10]])) - # array([[0, 1, 10], - # [0, 1, 10]]) + data_arrays.append(np.array([ + [0, 1, 10], + [0, 1, 10] + ])) if number_of_arrays > 3: - data_arrays.append([[0, 1, 10, 100], [0, 1, 10, 100]]) - # [[0, 1, 10, 100], - # [0, 1, 10, 100]] + data_arrays.append([ + [0, 1, 10, 100], + [0, 1, 10, 100] + ]) + # fmt: on return data_arrays diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index 5fe6e764bf..b02ce5046d 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -396,7 +396,7 @@ def test_get_gpw_file_path_pass(self): self.assertIn("gpw_v4_population", str(path)) except FileExistsError as err: self.assertIn("lease download", err.args[0]) - # pylint: disabled=consider-using-f-string + # pylint: disable=consider-using-f-string self.skipTest("GPW input data for GPW v4.%i not found." % (gpw_version)) def test_load_gpw_pop_shape_pass(self): @@ -416,7 +416,7 @@ def test_load_gpw_pop_shape_pass(self): self.assertEqual(len(data.shape), 2) except FileExistsError as err: self.assertIn("lease download", err.args[0]) - # pylint: disabled=consider-using-f-string + # pylint: disable=consider-using-f-string self.skipTest("GPW input data for GPW v4.%i not found." % (gpw_version)) From dbc5ceb8cfabf74a9bd282a5f4008318c4cb46ac Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Thu, 8 May 2025 16:17:41 +0200 Subject: [PATCH 11/18] name arguments --- climada/entity/exposures/litpop/litpop.py | 90 +++++++++++------------ 1 file changed, 45 insertions(+), 45 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index 811e28c2aa..0f7030084f 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -213,15 +213,15 @@ def from_countries( # GSDP data folder. litpop_list = [ _calc_admin1_one_country( - country, - res_arcsec, - exponents, - fin_mode, - tot_value, - reference_year, - gpw_version, - data_dir, - target_grid, + country=country, + res_arcsec=res_arcsec, + exponents=exponents, + fin_mode=fin_mode, + total_value=tot_value, + reference_year=reference_year, + gpw_version=gpw_version, + data_dir=data_dir, + target_grid=target_grid, ) for tot_value, country in zip(total_values, countries) ] @@ -231,7 +231,7 @@ def from_countries( # within each country and combined at the end. litpop_list = [ cls._from_country( - country, + country=country, exponents=exponents, fin_mode=fin_mode, total_value=total_values[idc], @@ -382,7 +382,7 @@ def from_nightlight_intensity( ) if countries is not None: exp = cls.from_countries( - countries, + countries=countries, res_arcsec=res_arcsec, exponents=(1, 0), fin_mode="none", @@ -393,8 +393,8 @@ def from_nightlight_intensity( ) else: exp = cls.from_shape( - shape, - None, + shape=shape, + total_value=None, res_arcsec=res_arcsec, exponents=(1, 0), value_unit="", @@ -496,7 +496,7 @@ def from_population( ) if countries is not None: exp = cls.from_countries( - countries, + countries=countries, res_arcsec=res_arcsec, exponents=(0, 1), fin_mode="pop", @@ -507,8 +507,8 @@ def from_population( ) else: exp = cls.from_shape( - shape, - None, + shape=shape, + total_value=None, res_arcsec=res_arcsec, exponents=(0, 1), value_unit="people", @@ -617,7 +617,7 @@ def from_shape_and_countries( """ # init countries' exposure: exp = cls.from_countries( - countries, + countries=countries, res_arcsec=res_arcsec, exponents=exponents, fin_mode=fin_mode, @@ -788,13 +788,13 @@ def from_shape( ) litpop_gdf, _ = _get_litpop_single_polygon( - shape, - reference_year, - data_dir, - gpw_version, - exponents, - target_grid, - region_id, + polygon=shape, + reference_year=reference_year, + data_dir=data_dir, + gpw_version=gpw_version, + exponents=exponents, + target_grid=target_grid, + region_id=region_id, ) # disaggregate total value proportional to LitPop values: @@ -977,12 +977,12 @@ def _from_country( gdf_tmp, meta_tmp, ) = _get_litpop_single_polygon( - polygon, - reference_year, - data_dir, - gpw_version, - exponents, - target_grid, + polygon=polygon, + reference_year=reference_year, + data_dir=data_dir, + gpw_version=gpw_version, + exponents=exponents, + target_grid=target_grid, verbose=(idx > 0), region_id=iso3n, ) @@ -1064,12 +1064,12 @@ def _crop_target_grid(target_grid, polygon): cropped_grid["width"] = col_max - col_min cropped_grid["height"] = row_max - row_min cropped_grid["transform"] = Affine( - target_grid["transform"].a, - 0, - cropped_min_lon, - 0, - target_grid["transform"].e, - cropped_max_lat, + a=target_grid["transform"].a, + b=0, + c=cropped_min_lon, + d=0, + e=target_grid["transform"].e, + f=cropped_max_lat, ) return cropped_grid @@ -1133,8 +1133,8 @@ def _get_litpop_single_polygon( offsets = (0, 0) if exponents[1] == 0 else (1, 0) # import population data (2d array), meta data, and global grid info, pop, meta_pop, _ = pop_util.load_gpw_pop_shape( - polygon, - reference_year, + polygon=polygon, + reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, verbose=verbose, @@ -1149,8 +1149,8 @@ def _get_litpop_single_polygon( # reproject Lit and Pop input data to aligned grid with target resolution: try: [pop, nlight], meta_out = reproject_input_data( - [pop, nlight], - [meta_pop, meta_nl], + data_array_list=[pop, nlight], + meta_list=[meta_pop, meta_nl], target_grid=cropped_grid, ) except ValueError as err: @@ -1363,9 +1363,9 @@ def reproject_input_data( dst_height, dst_width, dst_transform ) data_out_list[idx], meta_out["transform"] = u_coord.align_raster_data( - data_array_list[idx], - meta_list[idx]["crs"], - meta_list[idx]["transform"], + source=data_array_list[idx], + src_crs=meta_list[idx]["crs"], + src_transform=meta_list[idx]["transform"], dst_crs=dst_crs, dst_resolution=(res_degree, res_degree), dst_bounds=dst_bounds, @@ -1648,8 +1648,8 @@ def _calc_admin1_one_country( # total value is defined from country multiplied by grp_share: exp_list.append( LitPop.from_shape( - admin1_shapes[idx], - total_value * grp_values[record["name"]], + shape=admin1_shapes[idx], + total_value=total_value * grp_values[record["name"]], res_arcsec=res_arcsec, exponents=exponents, reference_year=reference_year, From b0d61ab1c4bb4380a128d7ba9daeaf4e2dab566e Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Thu, 8 May 2025 16:25:02 +0200 Subject: [PATCH 12/18] cosmetics --- climada/entity/exposures/test/test_litpop.py | 1 - climada/test/test_litpop_integr.py | 3 ++- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/climada/entity/exposures/test/test_litpop.py b/climada/entity/exposures/test/test_litpop.py index 2e0fd0dca8..048819e5ff 100644 --- a/climada/entity/exposures/test/test_litpop.py +++ b/climada/entity/exposures/test/test_litpop.py @@ -429,5 +429,4 @@ def test_get_value_unit_pass(self): if __name__ == "__main__": TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPop) - # TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestUncertainty)) unittest.TextTestRunner(verbosity=2).run(TESTS) diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index b02ce5046d..e19f10c938 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -586,5 +586,6 @@ def test_litpop_grid_consistency(self): if __name__ == "__main__": TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopExposure) TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAdmin1)) - TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopGridAlignment) + TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestGPWPopulation)) + TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestLitPopGridAlignment)) unittest.TextTestRunner(verbosity=2).run(TESTS) From 2415576a94add8e75a2cc287bf02f1af61b37cf9 Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Fri, 9 May 2025 10:22:26 +0200 Subject: [PATCH 13/18] name arguments correctly --- climada/entity/exposures/litpop/litpop.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index 0f7030084f..b45d114e62 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -1133,7 +1133,7 @@ def _get_litpop_single_polygon( offsets = (0, 0) if exponents[1] == 0 else (1, 0) # import population data (2d array), meta data, and global grid info, pop, meta_pop, _ = pop_util.load_gpw_pop_shape( - polygon=polygon, + geometry=polygon, reference_year=reference_year, gpw_version=gpw_version, data_dir=data_dir, From 272ff46339e702756df7a0f54f5e3b6d0d718feb Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Fri, 9 May 2025 11:25:54 +0200 Subject: [PATCH 14/18] move target_grid default into from_shape --- climada/entity/exposures/litpop/litpop.py | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index b45d114e62..51eb9eb574 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -373,13 +373,6 @@ def from_nightlight_intensity( raise ValueError( "Not allowed to set both `countries` and `shape`. Aborting." ) - if target_grid is None: - target_grid = cls._define_target_grid( - reference_year=reference_year, - gpw_version=GPW_VERSION, - data_dir=data_dir, - res_arcsec=res_arcsec, - ) if countries is not None: exp = cls.from_countries( countries=countries, @@ -487,13 +480,6 @@ def from_population( raise ValueError( "Not allowed to set both `countries` and `shape`. Aborting." ) - if target_grid is None: - target_grid = cls._define_target_grid( - reference_year=reference_year, - gpw_version=gpw_version, - data_dir=data_dir, - res_arcsec=res_arcsec, - ) if countries is not None: exp = cls.from_countries( countries=countries, @@ -787,6 +773,14 @@ def from_shape( "GeoSeries. Loop over elements of series outside method." ) + if target_grid is None: + target_grid = cls._define_target_grid( + reference_year=reference_year, + gpw_version=gpw_version, + data_dir=data_dir, + res_arcsec=res_arcsec, + ) + litpop_gdf, _ = _get_litpop_single_polygon( polygon=shape, reference_year=reference_year, From a93aa3f88c3d2f1988ebd79f62a0e0fa1b615a27 Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Fri, 9 May 2025 13:46:48 +0200 Subject: [PATCH 15/18] fix some tests --- climada/test/test_litpop_integr.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index e19f10c938..e408279303 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -167,7 +167,7 @@ def test_from_shape_zurich_pass(self): ent = lp.LitPop.from_shape( shape, total_value, res_arcsec=30, reference_year=2016 ) - self.assertEqual(ent.value.sum(), 1000.0) + self.assertAlmostEqual(ent.value.sum(), 1000.0) self.assertEqual(ent.value.min(), 0.0) self.assertEqual(ent.region_id.min(), 756) self.assertEqual(ent.region_id.max(), 756) @@ -223,7 +223,7 @@ def test_Liechtenstein_15_lit_pass(self): ref_year = 2016 ent = lp.LitPop.from_nightlight_intensity(country_name, reference_year=ref_year) - self.assertEqual(ent.value.sum(), 36469.0) + self.assertAlmostEqual(ent.value.sum(), 36469.0) self.assertEqual(ent.region_id[1], 438) self.assertEqual(ent.value_unit, "") self.assertAlmostEqual(ent.latitude.max(), 47.260416666666664) @@ -330,7 +330,15 @@ def test_calc_admin1(self): resolution = 300 country = "CHE" ent = lp._calc_admin1_one_country( - country, resolution, (2, 1), "pc", None, 2016, lp.GPW_VERSION, SYSTEM_DIR + country, + resolution, + exponents=(2, 1), + fin_mode="pc", + total_value=None, + reference_year=2016, + gpw_version=lp.GPW_VERSION, + data_dir=SYSTEM_DIR, + target_grid=None, ) self.assertEqual(ent.gdf.shape[0], 699) From ce30fa48db177a066e9ce94701fdd9008447eee2 Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Tue, 13 May 2025 14:10:03 +0200 Subject: [PATCH 16/18] distribute parts of _define_target_grid to respective helper modules --- .../entity/exposures/litpop/gpw_population.py | 51 +++++++++++++++++ climada/entity/exposures/litpop/litpop.py | 55 ++++--------------- climada/entity/exposures/litpop/nightlight.py | 17 +++++- climada/test/test_litpop_integr.py | 4 +- 4 files changed, 79 insertions(+), 48 deletions(-) diff --git a/climada/entity/exposures/litpop/gpw_population.py b/climada/entity/exposures/litpop/gpw_population.py index fbb34b4649..6d3f8306cd 100644 --- a/climada/entity/exposures/litpop/gpw_population.py +++ b/climada/entity/exposures/litpop/gpw_population.py @@ -23,6 +23,7 @@ import numpy as np import rasterio +from affine import Affine from climada import CONFIG from climada.util.constants import SYSTEM_DIR @@ -171,3 +172,53 @@ def get_gpw_file_path(gpw_version, reference_year, data_dir=None, verbose=True): f" different folder. The data can be downloaded from {sedac_browse_url}, e.g.," f" {sedac_file_url} (Free NASA Earthdata login required)." ) + + +def grid_aligned_with_gpw(reference_year, gpw_version, res_arcsec, data_dir=SYSTEM_DIR): + """ + Defines a grid based on population metadata. + + Parameters + ---------- + reference_year : int + The reference year for population and nightlight data. + gpw_version : int + Version number of GPW population data. + res_arcsec : int or None + Desired resolution in arcseconds. If None, aligns to population grid. + data_dir : str + Path to input data directory. + + Returns + ------- + grid : dict + A dictionary containing grid metadata, following the raster grid + specification. + """ + res_deg = res_arcsec / 3600 + + file_path = get_gpw_file_path( + gpw_version, reference_year, data_dir=data_dir, verbose=False + ) + with rasterio.open(file_path, "r") as src: + global_crs = src.crs + gpw_transform = src.transform + # Align grid resolution with GPW dataset + aligned_lon_min = -180 + (round((gpw_transform[2] - (-180)) / res_deg) * res_deg) + aligned_lat_max = 90 - (round((90 - gpw_transform[5]) / res_deg) * res_deg) + + global_transform = Affine(res_deg, 0, aligned_lon_min, 0, -res_deg, aligned_lat_max) + + global_width = round(360 / res_deg) + global_height = round(180 / res_deg) + + # Define the target grid using the computed values + return { + "driver": "GTiff", + "dtype": "float32", + "nodata": None, + "crs": global_crs, + "width": global_width, + "height": global_height, + "transform": global_transform, + } diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index 51eb9eb574..b3ec6e7b2b 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -33,10 +33,11 @@ import climada.util.finance as u_fin from climada import CONFIG from climada.entity.exposures.base import DEF_REF_YEAR, INDICATOR_IMPF, Exposures -from climada.entity.exposures.litpop import gpw_population as pop_util -from climada.entity.exposures.litpop import nightlight as nl_util from climada.util.constants import SYSTEM_DIR +from .gpw_population import grid_aligned_with_gpw, load_gpw_pop_shape +from .nightlight import BM_NIGHTLIGHT_GRID, load_nasa_nl_shape + LOGGER = logging.getLogger(__name__) GPW_VERSION = CONFIG.exposures.litpop.gpw_population.gpw_version.int() @@ -868,46 +869,14 @@ def _define_target_grid(reference_year, gpw_version, data_dir, res_arcsec): A dictionary containing grid metadata, following the raster grid specification. """ - res_deg = res_arcsec / 3600 if res_arcsec == 15: - # Black Marble Nightlights Grid - # pylint: disable=c-extension-no-member - global_crs = rasterio.crs.CRS.from_epsg(4326) # WGS84 projection - global_transform = Affine(res_deg, 0, -180, 0, -res_deg, 90) - global_width = round(360 / res_deg) - global_height = round(180 / res_deg) - else: - # GPW Population Grid - file_path = pop_util.get_gpw_file_path( - gpw_version, reference_year, data_dir=data_dir, verbose=False - ) - with rasterio.open(file_path, "r") as src: - global_crs = src.crs - gpw_transform = src.transform - # Align grid resolution with GPW dataset - aligned_lon_min = -180 + ( - round((gpw_transform[2] - (-180)) / res_deg) * res_deg - ) - aligned_lat_max = 90 - (round((90 - gpw_transform[5]) / res_deg) * res_deg) - - global_transform = Affine( - res_deg, 0, aligned_lon_min, 0, -res_deg, aligned_lat_max - ) - - global_width = round(360 / res_deg) - global_height = round(180 / res_deg) - - # Define the target grid using the computed values - target_grid = { - "driver": "GTiff", - "dtype": "float32", - "nodata": None, - "crs": global_crs, - "width": global_width, - "height": global_height, - "transform": global_transform, - } - return target_grid + return BM_NIGHTLIGHT_GRID + return grid_aligned_with_gpw( + reference_year=reference_year, + gpw_version=gpw_version, + res_arcsec=res_arcsec, + data_dir=data_dir, + ) @staticmethod def _from_country( @@ -1126,7 +1095,7 @@ def _get_litpop_single_polygon( # set nightlight offset (delta) to 1 in case n>0, c.f. delta in Eq. 1 of paper: offsets = (0, 0) if exponents[1] == 0 else (1, 0) # import population data (2d array), meta data, and global grid info, - pop, meta_pop, _ = pop_util.load_gpw_pop_shape( + pop, meta_pop, _ = load_gpw_pop_shape( geometry=polygon, reference_year=reference_year, gpw_version=gpw_version, @@ -1135,7 +1104,7 @@ def _get_litpop_single_polygon( ) total_population = pop.sum() # import nightlight data (2d array) and associated meta data: - nlight, meta_nl = nl_util.load_nasa_nl_shape( + nlight, meta_nl = load_nasa_nl_shape( polygon, reference_year, data_dir=data_dir, dtype=float ) diff --git a/climada/entity/exposures/litpop/nightlight.py b/climada/entity/exposures/litpop/nightlight.py index d875b26a9b..c24208507c 100644 --- a/climada/entity/exposures/litpop/nightlight.py +++ b/climada/entity/exposures/litpop/nightlight.py @@ -31,6 +31,7 @@ import numpy as np import rasterio import scipy.sparse as sparse +from affine import Affine from osgeo import gdal from PIL import Image from shapefile import Shape @@ -69,6 +70,17 @@ ] """Nightlight NASA files which generate the whole earth when put together.""" +BM_NIGHTLIGHT_GRID = { + "driver": "GTiff", + "dtype": "float32", + "nodata": None, + "crs": rasterio.crs.CRS.from_epsg(4326), # pylint: disable=c-extension-no-member + "width": round(360 / NASA_RESOLUTION_DEG), + "height": round(180 / NASA_RESOLUTION_DEG), + "transform": Affine(NASA_RESOLUTION_DEG, 0, -180, 0, -NASA_RESOLUTION_DEG, 90), +} +"""Grid aligned with Nightlight NASA files""" + def load_nasa_nl_shape(geometry, year, data_dir=SYSTEM_DIR, dtype="float32"): """Read nightlight data from NASA BlackMarble tiles @@ -146,6 +158,7 @@ def load_nasa_nl_shape(geometry, year, data_dir=SYSTEM_DIR, dtype="float32"): if idx == 0: meta = meta_tmp # set correct CRS from local tile's CRS to global WGS 84: + # pylint: disable=c-extension-no-member meta.update({"crs": rasterio.crs.CRS.from_epsg(4326), "dtype": dtype}) if len(req_files) == 1: # only one tile required: return np.array(out_image, dtype=dtype), meta @@ -264,9 +277,7 @@ def check_nl_local_file_exists(required_files=None, check_path=SYSTEM_DIR, year= required_files = np.ones( np.count_nonzero(BM_FILENAMES), ) - LOGGER.warning( - "The parameter 'required_files' was too short and " "is ignored." - ) + LOGGER.warning("The parameter 'required_files' was too short and is ignored.") if isinstance(check_path, str): check_path = Path(check_path) if not check_path.is_dir(): diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index e408279303..108baaea49 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -44,7 +44,7 @@ ) -class TestLitPopExposure(unittest.TestCase): +class TestLitPopExposures(unittest.TestCase): """Test LitPop exposure data model:""" def test_netherlands150_pass(self): @@ -592,7 +592,7 @@ def test_litpop_grid_consistency(self): if __name__ == "__main__": - TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopExposure) + TESTS = unittest.TestLoader().loadTestsFromTestCase(TestLitPopExposures) TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAdmin1)) TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestGPWPopulation)) TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestLitPopGridAlignment)) From 8944221a53a54066cbd3853ac42315ed4cc318ea Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Mon, 2 Jun 2025 18:32:09 +0200 Subject: [PATCH 17/18] make test_from_shape_zurich_pass pass --- climada/entity/exposures/litpop/litpop.py | 4 ++-- climada/test/test_litpop_integr.py | 14 ++++++++++---- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/climada/entity/exposures/litpop/litpop.py b/climada/entity/exposures/litpop/litpop.py index b3ec6e7b2b..9d25178378 100644 --- a/climada/entity/exposures/litpop/litpop.py +++ b/climada/entity/exposures/litpop/litpop.py @@ -1285,8 +1285,8 @@ def reproject_input_data( ... 'height': 1939, ... 'count': 1, ... 'crs': CRS.from_epsg(4326), - ... 'transform': Affine(0.00833333333333333, 0.0, -18.175000000000068, - ... 0.0, -0.00833333333333333, 43.79999999999993), + ... 'transform': Affine(30/3600, 0.0, -18.175, + ... 0.0, -30/3600, 43.8), ... } The meta data with the reference grid used to define the global destination diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index 108baaea49..617ce0a106 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -169,24 +169,30 @@ def test_from_shape_zurich_pass(self): ) self.assertAlmostEqual(ent.value.sum(), 1000.0) self.assertEqual(ent.value.min(), 0.0) + self.assertAlmostEqual(ent.value.max(), 5.058, places=4) + # Note: up to climada 6.0 this value used to be 5.05792616746308, then 5.058035160766561. + # The discrepancy assumedly caused by rounding differences of the slightly changed + # calculation procedure. self.assertEqual(ent.region_id.min(), 756) self.assertEqual(ent.region_id.max(), 756) - self.assertAlmostEqual(ent.latitude.min(), 47.20416666666661) + self.assertAlmostEqual(ent.latitude.min(), 47.2 + 15 / 3600) # index and coord. of largest value: self.assertEqual( - ent.gdf.loc[ent.gdf["value"] == ent.gdf["value"].max()].index[0], 482 + ent.gdf.loc[ent.gdf["value"] == ent.gdf["value"].max()].index[0], 434 ) + # Note: up to climada 6.0 this index used to be 482, because of its then irregular order. + # Specifically: the dataframe started at 36 and then skipped 71, 107, 143, etc. self.assertAlmostEqual( ent.gdf.loc[ent.gdf["value"] == ent.gdf["value"].max()].geometry.y.values[ 0 ], - 47.34583333333325, + 5681.5 * 30 / 3600, ) self.assertAlmostEqual( ent.gdf.loc[ent.gdf["value"] == ent.gdf["value"].max()].geometry.x.values[ 0 ], - 8.529166666666658, + 1023.5 * 30 / 3600, ) def test_from_shape_and_countries_zurich_pass(self): From a87304fe4d2552c232887769efa6fb2f14bbb2ad Mon Sep 17 00:00:00 2001 From: emanuel-schmid Date: Mon, 2 Jun 2025 18:33:36 +0200 Subject: [PATCH 18/18] test_from_shape_zurich_pass: remove notes about changed excpectations --- climada/test/test_litpop_integr.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/climada/test/test_litpop_integr.py b/climada/test/test_litpop_integr.py index 617ce0a106..6b93b47b42 100644 --- a/climada/test/test_litpop_integr.py +++ b/climada/test/test_litpop_integr.py @@ -170,9 +170,6 @@ def test_from_shape_zurich_pass(self): self.assertAlmostEqual(ent.value.sum(), 1000.0) self.assertEqual(ent.value.min(), 0.0) self.assertAlmostEqual(ent.value.max(), 5.058, places=4) - # Note: up to climada 6.0 this value used to be 5.05792616746308, then 5.058035160766561. - # The discrepancy assumedly caused by rounding differences of the slightly changed - # calculation procedure. self.assertEqual(ent.region_id.min(), 756) self.assertEqual(ent.region_id.max(), 756) self.assertAlmostEqual(ent.latitude.min(), 47.2 + 15 / 3600) @@ -180,8 +177,6 @@ def test_from_shape_zurich_pass(self): self.assertEqual( ent.gdf.loc[ent.gdf["value"] == ent.gdf["value"].max()].index[0], 434 ) - # Note: up to climada 6.0 this index used to be 482, because of its then irregular order. - # Specifically: the dataframe started at 36 and then skipped 71, 107, 143, etc. self.assertAlmostEqual( ent.gdf.loc[ent.gdf["value"] == ent.gdf["value"].max()].geometry.y.values[ 0