Update GHSL processing to produce height and proportion residential

Author: collijk

src/rra_building_density/constants.py

@@ -95,6 +95,8 @@ def process_resources(self, resolution: str) -> tuple[str, str]:
     ),
 }
 
+LATEST_MICROSOFT_VERSION = MICROSOFT_VERSIONS["6"]
+
 
 class GHSLVersion(BuiltVersion):
     measure_map: ClassVar[dict[str, tuple[str, str]]] = {

src/rra_building_density/process/ghsl.py

@@ -1,92 +1,112 @@
 import click
 import numpy as np
+import rasterra as rt
 from rra_tools import jobmon
 
 from rra_building_density import cli_options as clio
 from rra_building_density import constants as bdc
-from rra_building_density import utils
 from rra_building_density.data import BuildingDensityData
+from rra_building_density.process import utils
 
 
 def format_ghsl_main(
     block_key: str,
-    measure: str,
     time_point: str,
     resolution: int | str,
     output_dir: str,
 ) -> None:
-    ghsl_version = bdc.GHSL_VERSIONS["r2023a"]
-    ghsl_measure = ghsl_version.prefix_and_measure(measure)[1]
+    version = bdc.GHSL_VERSIONS["r2023a"]
     crs = bdc.CRSES["wgs84"].to_pyproj()
-
     bd_data = BuildingDensityData(output_dir)
+
+    print("Loading the reference block")
+    mask_version = bdc.LATEST_MICROSOFT_VERSION
+    reference_block = bd_data.load_tile(
+        resolution,
+        provider=mask_version.name,
+        block_key=block_key,
+        # This will be the transition time point
+        time_point=mask_version.time_points[0],
+        measure="density",
+    )
+
     print("Loading the tile index")
     tile_index = bd_data.load_tile_index(resolution)
-    tile_index_info = bd_data.load_tile_index_info(resolution)
-
-    print("Building template")
     block_index = tile_index[tile_index.block_key == block_key]
+
     block_poly_series = block_index.dissolve("block_key").geometry
-    block_poly = block_poly_series.iloc[0]
     block_poly_ghsl = (
         utils.bbox_safe_buffer(block_poly_series, 5000).to_crs(crs).iloc[0]
     )
 
-    block_template = utils.make_raster_template(
-        block_poly,
-        resolution=tile_index_info.tile_resolution,
-        crs=bdc.CRSES["equal_area"],
+    print("Loading the GHSL data")
+    density_arr, raw_volume_arr, raw_nonresidential_density_arr = (
+        utils.load_and_resample_ghsl_data(  # noqa: SLF001
+            measure=measure,
+            time_point=time_point,
+            ghsl_version=version,
+            bounds=block_poly_ghsl,
+            reference_block=reference_block,
+            bd_data=bd_data,
+        )._ndarray
+        for measure in ["density", "volume", "nonresidential_density"]
     )
-
-    print("Loading GHSL data")
-    raw_tile = bd_data.load_provider_tile(
-        ghsl_version,
-        bounds=block_poly_ghsl,
-        measure=ghsl_measure,
-        time_point=time_point,
-        year=time_point[:4],
+    print("Generating height and proportion residential arrays")
+    height_arr = utils.generate_height_array(density_arr, raw_volume_arr)  # type: ignore[arg-type]
+    proportion_residential_arr = utils.generate_proportion_residential_array(
+        density_arr,  # type: ignore[arg-type]
+        raw_nonresidential_density_arr,  # type: ignore[arg-type]
     )
-    raw_tile = raw_tile.astype(np.float32) / 10000.0
 
-    print("Resampling")
-    tile = raw_tile.set_no_data_value(np.nan).resample_to(block_template, "average")
-    tile = utils.suppress_noise(tile)
-    print("Saving")
-    bd_data.save_tile(
-        tile,
-        resolution,
-        provider="ghsl_r2023a",
-        block_key=block_key,
-        time_point=time_point,
-        measure=measure,
-    )
+    print("Generating and saving output arrays")
+    out_ops = {
+        "height": lambda _, h, __: h,
+        "proportion_residential": lambda _, __, p: p,
+        "density": lambda d, _, __: d,
+        "residential_density": lambda d, _, p: d * p,
+        "nonresidential_density": lambda d, _, p: d * (1 - p),
+        "volume": lambda d, h, _: h * d,
+        "residential_volume": lambda d, h, p: h * d * p,
+        "nonresidential_volume": lambda d, h, p: h * d * (1 - p),
+    }
+    for measure, op in out_ops.items():
+        out = rt.RasterArray(
+            data=op(density_arr, height_arr, proportion_residential_arr),  # type: ignore[no-untyped-call]
+            transform=reference_block.transform,
+            crs=reference_block.crs,
+            no_data_value=np.nan,
+        )
+        bd_data.save_tile(
+            out,
+            resolution,
+            provider="ghsl_r2023a",
+            block_key=block_key,
+            time_point=time_point,
+            measure=measure,
+        )
 
 
 @click.command()
-@clio.with_measure(bdc.GHSLVersion.measure_map)
 @clio.with_block_key()
 @clio.with_time_point()
 @clio.with_resolution(bdc.RESOLUTIONS)
 @clio.with_output_directory(bdc.MODEL_ROOT)
 def format_ghsl_task(
     block_key: str,
-    measure: str,
     time_point: str,
     resolution: str,
     output_dir: str,
 ) -> None:
     """Build predictors for a given tile and time point."""
-    format_ghsl_main(block_key, measure, time_point, resolution, output_dir)
+    format_ghsl_main(block_key, time_point, resolution, output_dir)
 
 
 @click.command()
-@clio.with_measure(bdc.GHSLVersion.measure_map, allow_all=True)
 @clio.with_time_point(allow_all=True)
 @clio.with_resolution(bdc.RESOLUTIONS)
 @clio.with_output_directory(bdc.MODEL_ROOT)
 @clio.with_queue()
 def format_ghsl(
-    measure: list[str],
     time_point: str,
     resolution: str,
     output_dir: str,
@@ -100,7 +120,7 @@ def format_ghsl(
     print("Loading the tile index")
     tile_index = bd_data.load_tile_index(resolution)
     block_keys = tile_index.block_key.unique().tolist()
-    njobs = len(block_keys) * len(time_point) * len(measure)
+    njobs = len(block_keys) * len(time_points)
     print(f"Formating building density for {njobs} block-times")
 
     memory, runtime = ghsl_version.process_resources(resolution)
@@ -114,7 +134,6 @@ def format_ghsl(
         },
         node_args={
             "block-key": block_keys,
-            "measure": measure,
             "time-point": time_points,
         },
         task_resources={

src/rra_building_density/process/microsoft.py

@@ -1,15 +1,12 @@
 import click
 import numpy as np
 import rasterra as rt
-from affine import Affine
 from rra_tools import jobmon
 
 from rra_building_density import cli_options as clio
 from rra_building_density import constants as bdc
-from rra_building_density import utils
 from rra_building_density.data import BuildingDensityData
-
-USE_WATER_MASK = False
+from rra_building_density.process import utils
 
 
 def format_microsoft_main(
@@ -25,28 +22,19 @@ def format_microsoft_main(
     tile_index_info = bd_data.load_tile_index_info(resolution)
     msft_index = bd_data.load_provider_index(msft_version, "union")
 
-    block_index = tile_index[tile_index.block_key == block_key]
-    block_poly_series = block_index.dissolve("block_key").geometry
-    block_poly = block_poly_series.iloc[0]
-    block_poly_msft = block_poly_series.to_crs(msft_index.crs).iloc[0]
+    block_poly, block_poly_msft = utils.get_block_polys(
+        tile_index[tile_index.block_key == block_key], msft_index.crs
+    )
 
     block_template = utils.make_raster_template(
         block_poly,
         resolution=tile_index_info.tile_resolution,
         crs=bdc.CRSES["equal_area"],
     )
 
-    overlapping = msft_index.loc[
-        msft_index.intersects(block_poly_msft), "quad_name"
-    ].tolist()
-
-    msft_tile_keys = [
-        tile_key
-        for tile_key in overlapping
-        if bd_data.provider_tile_exists(
-            msft_version, tile_key=tile_key, time_point=time_point
-        )
-    ]
+    msft_tile_keys = utils.get_provider_tile_keys(
+        msft_index, block_poly_msft, msft_version, bd_data, time_point=time_point
+    )
 
     if not msft_tile_keys:
         print("No overlapping building tiles, likely open ocean.")
@@ -65,31 +53,8 @@ def format_microsoft_main(
         bd_tile = bd_data.load_provider_tile(
             msft_version, tile_key=tile_key, time_point=time_point
         )
-
-        # The resolution of the MSFT tiles has too many decimal points.
-        # This causes tiles slightly west of the antimeridian to cross
-        # over and really mucks up reprojection. We'll clip the values
-        # here to 5 decimal places (ie to 100 microns), explicitly
-        # rounding down. This reduces the width of the tile by
-        # 512*0.0001 = 0.05m or 50cm, enough to fix roundoff issues.
-        x_res, y_res = bd_tile.resolution
-        xmin, xmax, ymin, ymax = bd_tile.bounds
-        bd_tile._transform = Affine(  # noqa: SLF001
-            a=utils.precise_floor(x_res, 4),
-            b=0.0,
-            c=xmin,
-            d=0.0,
-            e=-utils.precise_floor(-y_res, 4),
-            f=ymax,
-        )
+        bd_tile = utils.fix_microsoft_tile(bd_tile)
         bd_tile = bd_tile.unset_no_data_value().set_no_data_value(np.nan)
-        if USE_WATER_MASK:
-            # mask out water
-            mask_version = bdc.MICROSOFT_VERSIONS["water_mask"]
-            mask = bd_data.load_provider_tile(
-                mask_version, tile_key=tile_key
-            ).to_numpy()
-            bd_tile._ndarray[mask] = np.nan  # noqa: SLF001
 
         reprojected_tile = bd_tile.reproject(
             dst_resolution=block_template.x_resolution,