r.proj: size parallel bands from a precomputed footprint grid

include/grass/defs/gprojects.h

@@ -9,6 +9,8 @@ int GPJ_transform(const struct pj_info *, const struct pj_info *,
 int GPJ_transform_array(const struct pj_info *, const struct pj_info *,
                         const struct pj_info *, int, double *, double *,
                         double *, int);
+void GPJ_clone_transform(const struct pj_info *, struct gpj_transform_clone *);
+void GPJ_free_transform_clone(struct gpj_transform_clone *);
 
 /* old API, to be removed */
 int pj_do_proj(double *, double *, const struct pj_info *,

include/grass/gprojects.h

@@ -48,6 +48,15 @@ struct pj_info {
     char *wkt;
 };
 
+/* Per-thread clone of a transform, filled by GPJ_clone_transform() and
+ * released by GPJ_free_transform_clone(). Bundles the cloned transform with
+ * the private PROJ context it was cloned into, so ownership is a single unit.
+ */
+struct gpj_transform_clone {
+    struct pj_info info;
+    PJ_CONTEXT *ctx;
+};
+
 struct gpj_datum {
     char *name, *longname, *ellps;
     double dx, dy, dz;

lib/proj/do_proj.c

@@ -1414,3 +1414,46 @@ int pj_do_transform(int count, double *x, double *y, double *h,
     }
     return ok;
 }
+
+/*!
+ * \brief Clone a transform into a fresh per-thread PROJ context
+ *
+ * PROJ transformation objects are not safe for concurrent use, so each thread
+ * needs its own. This fills \p clone with a copy of \p src whose PJ is cloned
+ * into a new private context. Release it with GPJ_free_transform_clone().
+ *
+ * Safe to call concurrently from multiple threads with the same \p src,
+ * provided \p src is not modified during the calls. Each call clones into its
+ * own new context and touches no shared mutable state.
+ *
+ * \param src source transform (as set up by GPJ_init_transform())
+ * \param[out] clone receives the per-thread clone (info plus private context)
+ */
+void GPJ_clone_transform(const struct pj_info *src,
+                         struct gpj_transform_clone *clone)
+{
+    clone->ctx = proj_context_create();
+    /* A failed context leaves the thread with no usable transform, so this
+     * aborts the run. */
+    if (clone->ctx == NULL)
+        G_fatal_error(_("proj_context_create() failed for a per-thread "
+                        "transform clone"));
+    clone->info = *src;
+    clone->info.pj = proj_clone(clone->ctx, src->pj);
+    if (clone->info.pj == NULL)
+        G_fatal_error(_("proj_clone() failed for a per-thread transform "
+                        "clone"));
+}
+
+/*!
+ * \brief Free a per-thread transform clone and its context
+ *
+ * \param clone clone filled by GPJ_clone_transform(); its cloned PJ is set to
+ *              NULL after release
+ */
+void GPJ_free_transform_clone(struct gpj_transform_clone *clone)
+{
+    proj_destroy(clone->info.pj);
+    proj_context_destroy(clone->ctx);
+    clone->info.pj = NULL;
+}

raster/r.proj/Makefile

@@ -2,7 +2,7 @@ MODULE_TOPDIR = ../..
 
 PGM = r.proj
 
-LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PARSONLIB)
+LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PARSONLIB) $(PROJLIB)
 DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP)
 
 EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB)

raster/r.proj/benchmark/benchmark_r_proj.py

@@ -0,0 +1,128 @@
+"""This is a benchmark script for r.proj thread scaling with grass.benchmark.
+
+This script sweeps through a raster size at a fixed memory setting, and then
+memory at a fixed size. It then plots time, speedup, and efficiency. Creates
+its own source raster and projects, so it runs standalone with
+grass --exec python benchmark_r_proj.py."""
+
+import os
+import tempfile
+
+from grass.exceptions import CalledModuleError, GrassError
+from grass.pygrass.modules import Module
+import grass.script as gs
+import grass.benchmark as bm
+
+# Baselines held fixed while one dimension is swept.
+BASE_MAPSIZE = 50e6  # cells
+BASE_MEMORY = 300  # MB
+MAPSIZES = [10e6, 50e6, 100e6]
+MEMORIES = [50, 100, 300, 1000]
+METRICS = ["time", "speedup", "efficiency"]
+MAX_NPROCS = 8
+REPEAT = 3
+
+SRC_PROJECT = "src4326"
+DST_PROJECT = "dst3857"
+INPUT = "benchmark_r_proj_reference"
+OUTPUT = "benchmark_r_proj"
+
+
+def main():
+    gisdbase = tempfile.mkdtemp(prefix="bench_r_proj_")
+    gs.create_project(os.path.join(gisdbase, SRC_PROJECT), epsg="4326")
+    gs.create_project(os.path.join(gisdbase, DST_PROJECT), epsg="3857")
+
+    # Sweep raster size at the baseline memory.
+    results = []
+    for mapsize in MAPSIZES:
+        benchmark(
+            gisdbase,
+            size=int(mapsize**0.5),
+            memory=BASE_MEMORY,
+            label=f"r.proj_{int(mapsize / 1e6)}M",
+            results=results,
+        )
+    plot(results, "rastersize")
+
+    # Sweep memory at the baseline raster size.
+    results = []
+    for memory in MEMORIES:
+        benchmark(
+            gisdbase,
+            size=int(BASE_MAPSIZE**0.5),
+            memory=memory,
+            label=f"r.proj_memory_{memory}MB",
+            results=results,
+        )
+    plot(results, "memory")
+
+
+def benchmark(gisdbase, size, memory, label, results):
+    generate_input(gisdbase, size)
+    with gs.setup.init(
+        os.path.join(gisdbase, DST_PROJECT), env=os.environ.copy()
+    ) as session:
+        env = session.env
+        # Output region from r.proj's own suggested bounds for this input.
+        text = gs.read_command(
+            "r.proj",
+            project=SRC_PROJECT,
+            mapset="PERMANENT",
+            dbase=gisdbase,
+            input=INPUT,
+            method="nearest",
+            flags="g",
+            env=env,
+        )
+        region = dict(token.split("=") for token in text.split())
+        gs.run_command("g.region", env=env, **region)
+
+        module = Module(
+            "r.proj",
+            project=SRC_PROJECT,
+            mapset="PERMANENT",
+            dbase=gisdbase,
+            input=INPUT,
+            output=OUTPUT,
+            method="nearest",
+            memory=memory,
+            env_=env,
+            run_=False,
+            overwrite=True,
+        )
+        results.append(
+            bm.benchmark_nprocs(
+                module, label=label, max_nprocs=MAX_NPROCS, repeat=REPEAT
+            )
+        )
+
+
+def generate_input(gisdbase, size):
+    """Generate the source raster in the EPSG:4326 project. Uses
+    r.surf.fractal, or r.random.surface when FFTW is unavailable."""
+    with gs.setup.init(
+        os.path.join(gisdbase, SRC_PROJECT), env=os.environ.copy()
+    ) as session:
+        env = session.env
+        gs.run_command(
+            "g.region", n=50, s=40, w=-110, e=-90, rows=size, cols=size, env=env
+        )
+        try:
+            Module("r.surf.fractal", output=INPUT, overwrite=True, env_=env)
+        except (CalledModuleError, GrassError):
+            Module("r.random.surface", output=INPUT, overwrite=True, env_=env)
+
+
+def plot(results, sweep):
+    for metric in METRICS:
+        bm.nprocs_plot(
+            results,
+            filename=f"r_proj_{sweep}_{metric}.svg",
+            title=f"r.proj {sweep} {metric}",
+            metric=metric,
+        )
+
+
+if __name__ == "__main__":
+    main()