Urban morphology

def urban_parameters_from_morphology(
    lcz_map: xr.DataArray,
    usm_driver: xr.Dataset,
    coarsen: int | None = None,
    pad: Dict[str, float] | None = None,
) -> xr.Dataset:
    zones = np.unique(lcz_map)
    params = {}
    dx_usm = float(usm_driver.x[1] - usm_driver.x[0])
    dy_usm = float(usm_driver.y[1] - usm_driver.y[0])
    nx_usm = usm_driver.x.size
    ny_usm = usm_driver.y.size
    extent_x_usm = dx_usm * nx_usm
    extent_y_usm = dy_usm * ny_usm

    for zone in zones:
        params[zone] = {}

        total_gridpoints = (lcz_map == zone).sum()
        total_area = total_gridpoints * dx_usm * dy_usm

        lcz_mask = lcz_map == zone
        building_mask = lcz_mask & (usm_driver["building_type"] > 0)
        street_mask = lcz_mask & (usm_driver["street_type"] > 0)

        # Mean building height
        params[zone]["building_height"] = float(
            usm_driver["buildings_2d"].where(building_mask).mean()
        )

        # Urban fraction
        params[zone]["urban_fraction"] = float((
            building_mask | street_mask
        ).sum() / total_gridpoints)

        urban_area = total_area * params[zone]["urban_fraction"]

        # Plan area fraction
        params[zone]["building_plan_area_fraction"] = float((
            building_mask.sum() / total_gridpoints
        ) / params[zone]["urban_fraction"] )

        # Compute total wall area
        wall_area = 0.0
        for i in range(building_mask.x.size):
            for j in range(building_mask.y.size):
                if building_mask[i, j]:
                    # Check neighbors
                    if i == 0 or not building_mask[i - 1, j]:
                        wall_area += (
                            usm_driver["buildings_2d"][i, j] * dy_usm
                        )  # Top edge
                    if (
                        i == building_mask.x.size - 1
                        or not building_mask[i + 1, j]
                    ):
                        wall_area += (
                            usm_driver["buildings_2d"][i, j] * dy_usm
                        )  # Bottom edge
                    if j == 0 or not building_mask[i, j - 1]:
                        wall_area += (
                            usm_driver["buildings_2d"][i, j] * dx_usm
                        )  # Left edge
                    if (
                        j == building_mask.y.size - 1
                        or not building_mask[i, j + 1]
                    ):
                        wall_area += (
                            usm_driver["buildings_2d"][i, j] * dx_usm
                        )  # Right edge

        wall_area_density = wall_area / urban_area

        # Street canyon aspect ratio
        params[zone]["street_canyon_aspect_ratio"] = float(
            0.5
            * wall_area_density
            / (1.0 - params[zone]["building_plan_area_fraction"])
        )

        # Frontal area fraction
        params[zone]["building_frontal_area_fraction"] = float(
            wall_area_density / 4
        )

    return params