distribute parts of _define_target_grid to respective helper modules

Author: emanuel-schmid

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,
+    }

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
     )
 

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():

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))