[napattan]

Changes proposed in this PR:

• Vectorize Impact.impact_at_reg in climada/engine/impact.py by replacing the per-region Python list comprehension ([self.imp_mat[:, np.where(agg_regions == reg)[0]].sum(1) for reg in agg_reg_unique]) with np.unique(agg_regions, return_inverse=True) and a single CSR sparse indicator matrix multiplication (self.imp_mat @ reg_ind).toarray().
• Vectorize the bounding-box construction inside points_to_raster in climada/util/coordinates.py (apply_box) by replacing row-wise df_exp.apply(lambda r: ..., axis=1) (Point(...).buffer(res / 2).envelope) with shapely.box(lon - res / 2, lat - res / 2, lon + res / 2, lat + res / 2).
• Remove the redundant df_geom test block in climada/util/test/test_coordinates.py and replace the duplicated list-comprehension check in climada/engine/test/test_impact.py (TestImpactReg.test_agg_regions) with direct per-event assertions against self.imp.at_event and self.imp.imp_mat.toarray().

Benchmark Verification Script & Results

Below is a self-contained script comparing develop vs. this PR across all 4 rows of the benchmark table (Impact.impact_at_reg end-to-end, points_to_raster apply_box with latitude/longitude columns, points_to_raster apply_box with geometry column, and full end-to-end points_to_raster including df_poly.to_crs and rasterio.features.rasterize):

import time
from unittest.mock import patch
import geopandas as gpd
import numpy as np
import pandas as pd
from scipy import sparse
import shapely
from shapely.geometry import Point

from climada.engine.impact import Impact
import climada.util.coordinates as u_coord


def bench_impact_at_reg(n_runs=10):
    rng = np.random.default_rng(42)
    n_events, n_exp, n_regions = 2_000, 50_000, 100
    imp_mat = sparse.random(
        n_events, n_exp, density=0.005, format="csr", random_state=rng, dtype=np.float64
    )
    agg_regions = rng.choice([f"REG_{i:03d}" for i in range(n_regions)], size=n_exp)

    imp = Impact(
        event_id=np.arange(1, n_events + 1),
        event_name=[f"ev_{i}" for i in range(1, n_events + 1)],
        date=np.arange(1, n_events + 1),
        frequency=np.ones(n_events) / n_events,
        coord_exp=np.zeros((n_exp, 2)),
        crs="EPSG:4326",
        eai_exp=np.zeros(n_exp),
        at_event=np.asarray(imp_mat.sum(axis=1)).ravel(),
        tot_value=1.0,
        aai_agg=0.0,
        unit="USD",
        imp_mat=imp_mat,
    )

    def develop_impact_at_reg(self_imp, agg_regs):
        agg_regs = np.asanyarray(agg_regs)
        agg_reg_unique = np.unique(agg_regs)
        at_reg_event = np.hstack(
            [
                self_imp.imp_mat[:, np.where(agg_regs == reg)[0]].sum(1)
                for reg in agg_reg_unique
            ]
        )
        return pd.DataFrame(
            at_reg_event, columns=agg_reg_unique, index=self_imp.event_id
        )

    df_old = develop_impact_at_reg(imp, agg_regions)
    df_new = imp.impact_at_reg(agg_regions)
    pd.testing.assert_frame_equal(df_old, df_new, rtol=1e-13, atol=1e-13)

    t_old = []
    for _ in range(n_runs):
        t0 = time.perf_counter()
        develop_impact_at_reg(imp, agg_regions)
        t_old.append((time.perf_counter() - t0) * 1000.0)

    t_new = []
    for _ in range(n_runs):
        t0 = time.perf_counter()
        imp.impact_at_reg(agg_regions)
        t_new.append((time.perf_counter() - t0) * 1000.0)

    print(
        f"1. Impact.impact_at_reg ({n_events}x{n_exp}, {n_regions} regions): "
        f"develop median={np.median(t_old):.2f} ms (mean={np.mean(t_old):.2f}±{np.std(t_old):.2f} ms) -> "
        f"PR #1350 median={np.median(t_new):.2f} ms (mean={np.mean(t_new):.2f}±{np.std(t_new):.2f} ms) "
        f"[{np.median(t_old) / np.median(t_new):.1f}x speedup]"
    )


def bench_points_to_raster(n_runs=7):
    gx, gy = np.meshgrid(np.linspace(5.0, 14.95, 200), np.linspace(45.0, 49.95, 100))
    rng = np.random.default_rng(42)
    vals = rng.uniform(1.0, 100.0, size=gx.size)
    df_val = gpd.GeoDataFrame(
        {"latitude": gy.ravel(), "longitude": gx.ravel(), "value": vals}
    )
    gdf_geom = gpd.GeoDataFrame(
        {"value": vals},
        geometry=gpd.points_from_xy(gx.ravel(), gy.ravel()),
        crs="EPSG:4326",
    )
    res = 0.05

    # 2a. Isolated apply_box step (latitude/longitude columns)
    fun_ll = lambda r: Point(r["longitude"], r["latitude"]).buffer(res / 2).envelope
    t_box_old_ll, t_box_new_ll = [], []
    for _ in range(n_runs):
        t0 = time.perf_counter()
        s_old_ll = df_val.apply(fun_ll, axis=1)
        t_box_old_ll.append((time.perf_counter() - t0) * 1000.0)
        t0 = time.perf_counter()
        lon, lat = df_val["longitude"], df_val["latitude"]
        s_new_ll = shapely.box(lon - res / 2, lat - res / 2, lon + res / 2, lat + res / 2)
        t_box_new_ll.append((time.perf_counter() - t0) * 1000.0)
    assert shapely.equals(s_old_ll.to_numpy(), s_new_ll.to_numpy()).all()
    print(
        f"2. apply_box lat/lon ({len(df_val)} pts): "
        f"develop median={np.median(t_box_old_ll):.2f} ms (mean={np.mean(t_box_old_ll):.2f}±{np.std(t_box_old_ll):.2f} ms) -> "
        f"PR #1350 median={np.median(t_box_new_ll):.2f} ms (mean={np.mean(t_box_new_ll):.2f}±{np.std(t_box_new_ll):.2f} ms) "
        f"[{np.median(t_box_old_ll) / np.median(t_box_new_ll):.1f}x speedup]"
    )

    # 2b. Isolated apply_box step (geometry column)
    fun_g = lambda r: r.geometry.buffer(res / 2).envelope
    t_box_old_g, t_box_new_g = [], []
    for _ in range(n_runs):
        t0 = time.perf_counter()
        s_old_g = gdf_geom.apply(fun_g, axis=1)
        t_box_old_g.append((time.perf_counter() - t0) * 1000.0)
        t0 = time.perf_counter()
        lon, lat = gdf_geom.geometry.x, gdf_geom.geometry.y
        s_new_g = shapely.box(lon - res / 2, lat - res / 2, lon + res / 2, lat + res / 2)
        t_box_new_g.append((time.perf_counter() - t0) * 1000.0)
    assert shapely.equals(s_old_g.to_numpy(), s_new_g.to_numpy()).all()
    print(
        f"3. apply_box geometry ({len(gdf_geom)} pts): "
        f"develop median={np.median(t_box_old_g):.2f} ms (mean={np.mean(t_box_old_g):.2f}±{np.std(t_box_old_g):.2f} ms) -> "
        f"PR #1350 median={np.median(t_box_new_g):.2f} ms (mean={np.mean(t_box_new_g):.2f}±{np.std(t_box_new_g):.2f} ms) "
        f"[{np.median(t_box_old_g) / np.median(t_box_new_g):.1f}x speedup]"
    )

    # 2c. Full end-to-end u_coord.points_to_raster call (includes apply_box + to_crs + rasterize)
    t_e2e_old, t_e2e_new = [], []
    with patch.object(
        u_coord.shapely,
        "box",
        side_effect=lambda xmin, ymin, xmax, ymax: df_val.apply(fun_ll, axis=1),
    ):
        r_old, _ = u_coord.points_to_raster(df_val, val_names=["value"], res=res, raster_res=res)
        for _ in range(n_runs):
            t0 = time.perf_counter()
            u_coord.points_to_raster(df_val, val_names=["value"], res=res, raster_res=res)
            t_e2e_old.append((time.perf_counter() - t0) * 1000.0)

    r_new, _ = u_coord.points_to_raster(df_val, val_names=["value"], res=res, raster_res=res)
    np.testing.assert_array_equal(r_old, r_new)
    for _ in range(n_runs):
        t0 = time.perf_counter()
        u_coord.points_to_raster(df_val, val_names=["value"], res=res, raster_res=res)
        t_e2e_new.append((time.perf_counter() - t0) * 1000.0)

    print(
        f"4. points_to_raster end-to-end ({len(df_val)} pts, res={res}): "
        f"develop median={np.median(t_e2e_old):.2f} ms (mean={np.mean(t_e2e_old):.2f}±{np.std(t_e2e_old):.2f} ms) -> "
        f"PR #1350 median={np.median(t_e2e_new):.2f} ms (mean={np.mean(t_e2e_new):.2f}±{np.std(t_e2e_new):.2f} ms) "
        f"[{np.median(t_e2e_old) / np.median(t_e2e_new):.1f}x speedup]"
    )


if __name__ == "__main__":
    bench_impact_at_reg()
    bench_points_to_raster()

Measured results (Python 3.13, NumPy 2.4, SciPy 1.17, Shapely 2.1, GeoPandas 1.1):

Target	Scope	develop (median / mean ± std)	PR #1350 (median / mean ± std)	Speedup (median)
Impact.impact_at_reg	End-to-end (2,000 events × 50,000 exposures × 100 regions, CSR density=0.005)	138.23 ms (138.79 ± 5.89 ms)	11.78 ms (12.08 ± 1.01 ms)	11.7x
points_to_raster (apply_box)	Isolated bounding-box creation (20,000 points, lat/lon columns)	510.35 ms (509.87 ± 13.36 ms)	7.30 ms (7.47 ± 0.86 ms)	69.9x
points_to_raster (apply_box)	Isolated bounding-box creation (20,000 points, geometry column)	451.23 ms (447.10 ± 12.45 ms)	10.28 ms (10.10 ± 1.87 ms)	43.9x
points_to_raster (full function)	End-to-end (apply_box + to_crs + rasterio.features.rasterize, 20,000 points)	1001.34 ms (1003.06 ± 21.27 ms)	483.83 ms (482.17 ± 27.66 ms)	2.1x

This PR fixes # (N/A — performance optimization)

PR Author Checklist

• Read the Contribution Guide and AI Policy
• Correct target branch selected (if unsure, select develop)
• Descriptive pull request title added
• Source branch up-to-date with target branch
• Documentation updated
• Tests updated
• Tests passing
• No new linter issues
• Changelog updated

PR Reviewer Checklist

• Read the Contribution Guide and AI Policy
• CLIMADA Reviewer Checklist passed
• Tests passing
• No new linter issues

[peanutfun]

@napattan Thank you for your contribution! Please provide Python scripts to verify your claims on the speedup. Please also use the provided pull request description template (https://github.com/CLIMADA-project/climada_python/blob/main/.github/pull_request_template.md) for a smoother review. Finally, you added a test that does not verify anything (new). Please revise before we can merge

[peanutfun]

This test does not make any sense. df_geom contains the same data as df_val. Then you call points_to_raster as above and check against the previous result.

[napattan]

Removed this redundant df_geom test block in the latest commit. Thank you for noticing.

[napattan]

Hi @peanutfun, thank you for taking the time to review.

My background is in landscape architecture, and I have been trying on integrating computational workflow and spatial analysis into my landscape and urban fields. Also I have used AI aided in my coding workflow and I read CLIMADA's AI policy markdown now, I take full responsibility for every line of code in this PR. I am still very new to contributing in open-source projects, I truly apologize for missing description template, issuing very unclearly, and for redundant test.

I have updated the PR description above using .github/pull_request_template.md along with the Python benchmark script and a breakdown of both the isolated apply_box speedup (69.9x) and the full end-to-end points_to_raster (2.1x) and Impact.impact_at_reg (11.7x) speedups, and pushed a commit removing the redundant df_geom test in test_coordinates.py.

If you checked and think this is unhelpful, please feel free to close this issue.
Thank you for your time. I promise to be more considerate from now on.

[peanutfun]

@napattan All good, I appreciate your frankness and I'll walk you through it. Thanks for providing the benchmark code. Could you simplify it by only including the end-to-end calculations with timings? The script does not have to provide the speedup, just the run time. Then one can simply switch between branches, re-run the script, and see the timings for the actual branch code.

The speedup for impact_at_reg looks impressive, but the speedup for the end-to-end points_to_raster is significantly lower than the isolated apply_box. I guess this is because the "boxing algorithm" is not the actual bottleneck of the function?