raster: Restrict GDAL-linked raster reads to the needed column range

lib/raster/R.h

@@ -53,9 +53,13 @@ struct fileinfo /* Information for opened cell files */
     struct Range range;      /* Range structure              */
     struct FPRange fp_range; /* float Range structure        */
     int want_histogram;
-    int reclass_flag;         /* Automatic reclass flag       */
-    off_t *row_ptr;           /* File row addresses           */
-    COLUMN_MAPPING *col_map;  /* Data to window col mapping   */
+    int reclass_flag;        /* Automatic reclass flag       */
+    off_t *row_ptr;          /* File row addresses           */
+    COLUMN_MAPPING *col_map; /* Data to window col mapping   */
+    /* Range of native columns of GDAL-linked maps needed with
+     * the current region. */
+    COLUMN_MAPPING gdal_min_col;
+    COLUMN_MAPPING gdal_max_col;
     double C1, C2;            /* Data to window row constants */
     int cur_row;              /* Current data row in memory   */
     int null_cur_row;         /* Current null row in memory   */

lib/raster/get_row.c

@@ -205,26 +205,33 @@ static void read_data_gdal(int fd, int row, unsigned char *data_buf,
     struct fileinfo *fcb = &R__.fileinfo[fd];
     unsigned char *buf;
     CPLErr err;
+    /* Logical (pre-flip) column range actually needed by the region;
+     * unrestricted (full row) if the window mapping left it unset. */
+    int min_col = fcb->gdal_min_col >= 0 ? fcb->gdal_min_col : 0;
+    int max_col =
+        fcb->gdal_min_col >= 0 ? fcb->gdal_max_col : fcb->cellhd.cols - 1;
+    int ncols = max_col - min_col + 1;
+    /* hflip'ed maps store columns mirrored, so the logical range read
+     * from disk is the physical range at the opposite end of the row. */
+    int col_off = fcb->gdal->hflip ? fcb->cellhd.cols - 1 - max_col : min_col;
 
     *nbytes = fcb->nbytes;
 
     if (fcb->gdal->vflip)
         row = fcb->cellhd.rows - 1 - row;
 
-    buf = fcb->gdal->hflip ? G_malloc(fcb->cellhd.cols * fcb->cur_nbytes)
-                           : data_buf;
+    buf = fcb->gdal->hflip ? G_malloc((size_t)ncols * fcb->cur_nbytes)
+                           : data_buf + (size_t)col_off * fcb->cur_nbytes;
 
-    err =
-        Rast_gdal_raster_IO(fcb->gdal->band, GF_Read, 0, row, fcb->cellhd.cols,
-                            1, buf, fcb->cellhd.cols, 1, fcb->gdal->type, 0, 0);
+    err = Rast_gdal_raster_IO(fcb->gdal->band, GF_Read, col_off, row, ncols, 1,
+                              buf, ncols, 1, fcb->gdal->type, 0, 0);
 
     if (fcb->gdal->hflip) {
         int i;
 
-        for (i = 0; i < fcb->cellhd.cols; i++)
-            memcpy(data_buf + i * fcb->cur_nbytes,
-                   buf + (fcb->cellhd.cols - 1 - i) * fcb->cur_nbytes,
-                   fcb->cur_nbytes);
+        for (i = 0; i < ncols; i++)
+            memcpy(data_buf + (min_col + i) * fcb->cur_nbytes,
+                   buf + (ncols - 1 - i) * fcb->cur_nbytes, fcb->cur_nbytes);
         G_free(buf);
     }
 

lib/raster/open.c

@@ -37,6 +37,8 @@ static int new_fileinfo(void)
         if (R__.fileinfo[i].open_mode <= 0) {
             memset(&R__.fileinfo[i], 0, sizeof(struct fileinfo));
             R__.fileinfo[i].open_mode = -1;
+            R__.fileinfo[i].gdal_min_col = -1;
+            R__.fileinfo[i].gdal_max_col = -1;
             return i;
         }
 
@@ -51,6 +53,8 @@ static int new_fileinfo(void)
     for (i = oldsize; i < newsize; i++) {
         memset(&R__.fileinfo[i], 0, sizeof(struct fileinfo));
         R__.fileinfo[i].open_mode = -1;
+        R__.fileinfo[i].gdal_min_col = -1;
+        R__.fileinfo[i].gdal_max_col = -1;
     }
 
     R__.fileinfo_count = newsize;

lib/raster/tests/lib_raster_gdal_link_window_test.py

@@ -0,0 +1,273 @@
+"""Tests for column-windowed reading of GDAL-linked (r.external) raster maps.
+
+Rast_get_row() reads GDAL-linked maps through read_data_gdal(), which
+restricts the GDAL read to the range of native columns that overlap the
+current region instead of always reading the full native row width.
+"""
+
+import os
+
+import numpy as np
+import pytest
+
+import grass.script as gs
+from grass.experimental import TemporaryMapsetSession
+from grass.script import array as garray
+from grass.tools import Tools
+
+ROWS = 20
+COLS = 30
+NULL = -999999
+# Geometry of the source file, matching the g.region call in linked_session.
+FILE_NORTH = ROWS
+FILE_WEST = 0
+FILE_RES = 1
+
+
+@pytest.fixture(scope="module")
+def linked_session(tmp_path_factory):
+    """Session with a GeoTIFF exported and linked as 'linked'.
+
+    The source raster has cell values that can be computed from their
+    row and column positions without needing a second, independent
+    read of the file.
+    """
+    project = tmp_path_factory.mktemp("gdal_link_window") / "project"
+    gs.create_project(project, epsg=3358)
+    tif_path = tmp_path_factory.mktemp("gdal_link_window_data") / "source.tif"
+    with gs.setup.init(project, env=os.environ.copy()) as session:
+        tools = Tools(session=session)
+        tools.g_region(n=ROWS, s=0, w=0, e=COLS, res=1)
+        tools.r_mapcalc(expression="source = (row() - 1) * 1000 + (col() - 1)")
+        tools.r_out_gdal(
+            input="source", output=str(tif_path), format="GTiff", type="Int32"
+        )
+        tools.r_external(input=str(tif_path), output="linked")
+        # Link also as flipped raster maps
+        tools.r_external(input=str(tif_path), output="linked_h", flags="h")
+        tools.r_external(input=str(tif_path), output="linked_v", flags="v")
+        tools.r_external(input=str(tif_path), output="linked_hv", flags="hv")
+        yield session, tif_path
+
+
+@pytest.fixture
+def session(linked_session):
+    """A session in its own temporary mapset, so each test has its own region."""
+    session, _ = linked_session
+    with TemporaryMapsetSession(env=session.env) as mapset_session:
+        yield mapset_session
+
+
+@pytest.fixture
+def source_tif(linked_session):
+    """Path to the GeoTIFF file linked as 'linked' in the session fixture."""
+    _, tif_path = linked_session
+    return tif_path
+
+
+@pytest.fixture(scope="module")
+def latlon_session(tmp_path_factory):
+    """WGS84 session with a full-longitude GeoTIFF linked as 'latlon'.
+
+    The region covers a -180 to 180 longitude range at 1 degree resolution,
+    to test wrapping of lat/lon in Rast__create_window_mapping() (window_map.c).
+    """
+    project = tmp_path_factory.mktemp("gdal_link_window_ll") / "project"
+    gs.create_project(project, epsg=4326)
+    tif_path = tmp_path_factory.mktemp("gdal_link_window_ll_data") / "source_ll.tif"
+    with gs.setup.init(project, env=os.environ.copy()) as session:
+        tools = Tools(session=session)
+        tools.g_region(n=5, s=-5, w=-180, e=180, res=1)
+        tools.r_mapcalc(expression="source_ll = (row() - 1) * 1000 + (col() - 1)")
+        tools.r_out_gdal(
+            input="source_ll", output=str(tif_path), format="GTiff", type="Int32"
+        )
+        tools.r_external(input=str(tif_path), output="latlon")
+        yield session
+
+
+@pytest.fixture
+def latlon_mapset(latlon_session):
+    """A session in its own temporary mapset, in the WGS84 project"""
+    with TemporaryMapsetSession(env=latlon_session.env) as mapset_session:
+        yield mapset_session
+
+
+def apply_flip(rows_idx, cols_idx, hflip, vflip):
+    """Map native (row, col) indices through a GDAL link's hflip/vflip.
+
+    Mirrors the row/column reversal read_data_gdal() (get_row.c) applies,
+    over the file's full row/column range, for maps linked with a flip.
+    """
+    if vflip:
+        rows_idx = ROWS - 1 - rows_idx
+    if hflip:
+        cols_idx = COLS - 1 - cols_idx
+    return rows_idx, cols_idx
+
+
+def expected_values(row_offset, col_offset, rows, cols, hflip=False, vflip=False):
+    """Expected values for a region shifted by row/col_offset cells."""
+    row_idx, col_idx = apply_flip(
+        row_offset + np.arange(rows), col_offset + np.arange(cols), hflip, vflip
+    )
+    return row_idx[:, None] * 1000 + col_idx[None, :]
+
+
+def nearest_native_index(offset, step, count):
+    """offset + i * step for i in range(count), floored.
+
+    Reproduces the nearest-neighbor mapping from a region cell to
+    a native file cell in Rast__create_window_mapping() (window_map.c).
+    """
+    return np.floor(offset + step * np.arange(count)).astype(int)
+
+
+def native_indices_for_region(north, west, res, rows, cols, hflip=False, vflip=False):
+    """Native (row, col) indices 'linked' resolves to for a region."""
+    step = res / FILE_RES
+    native_cols = nearest_native_index(
+        (west - FILE_WEST + res / 2.0) / FILE_RES, step, cols
+    )
+    native_rows = nearest_native_index(
+        (FILE_NORTH - north + res / 2.0) / FILE_RES, step, rows
+    )
+    return apply_flip(native_rows, native_cols, hflip, vflip)
+
+
+def wrapped_native_col_indices(region_west, region_east, res, file_west, file_cols):
+    """Native (0-based) column indices of GRASS's lat/lon wraparound mapping.
+
+    Mirrors Rast__create_window_mapping() (window_map.c).
+    """
+    west, east = region_west, region_east
+    while west > file_west + 360.0:
+        west -= 360.0
+        east -= 360.0
+    while west < file_west:
+        west += 360.0
+        east += 360.0
+
+    cols = round((region_east - region_west) / res)
+
+    def native_for(west):
+        x = np.floor((west - file_west + res / 2.0) / res + np.arange(cols))
+        x[(x < 0) | (x >= file_cols)] = -1
+        return x.astype(int)
+
+    native = native_for(west)
+    while east - 360.0 > file_west:
+        east -= 360.0
+        west -= 360.0
+        unresolved = native < 0
+        native[unresolved] = native_for(west)[unresolved]
+    return native
+
+
+# Configure parametrization
+FLIP_CASES = pytest.mark.parametrize(
+    ("hflip", "vflip", "raster_name"),
+    [
+        (False, False, "linked"),
+        (True, False, "linked_h"),
+        (False, True, "linked_v"),
+        (True, True, "linked_hv"),
+    ],
+    ids=["noflip", "hflip", "vflip", "hvflip"],
+)
+
+
+@FLIP_CASES
+def test_region_fully_inside_source_extent(session, hflip, vflip, raster_name):
+    """A region fully inside the file reads the correct sub-window."""
+    Tools(session=session).g_region(n=15, s=8, w=12, e=25, res=1)
+    arr = np.array(garray.array(raster_name, null=NULL, env=session.env))
+    assert np.array_equal(
+        arr, expected_values(5, 12, *arr.shape, hflip=hflip, vflip=vflip)
+    )
+
+
+@FLIP_CASES
+def test_region_partially_outside_source_extent(session, hflip, vflip, raster_name):
+    """Columns outside the file's extent read as null, the rest as data."""
+    Tools(session=session).g_region(n=10, s=5, w=-5, e=10, res=1)
+    arr = np.array(garray.array(raster_name, null=NULL, env=session.env))
+    assert np.all(arr[:, :5] == NULL)
+    assert np.array_equal(
+        arr[:, 5:], expected_values(10, 0, arr.shape[0], 10, hflip=hflip, vflip=vflip)
+    )
+
+
+@FLIP_CASES
+def test_region_fully_outside_source_extent(session, hflip, vflip, raster_name):
+    """A region with no overlap at all reads back as entirely null."""
+    Tools(session=session).g_region(n=10, s=5, w=-50, e=-40, res=1)
+    arr = np.array(garray.array(raster_name, null=NULL, env=session.env))
+    assert np.all(arr == NULL)
+
+
+@FLIP_CASES
+def test_region_coarser_than_source_resolution(session, hflip, vflip, raster_name):
+    """A region coarser than the file's resolution reads the nearest cell."""
+    Tools(session=session).g_region(n=16, s=6, w=10, e=24, res=2)
+    arr = np.array(garray.array(raster_name, null=NULL, env=session.env))
+    native_rows, native_cols = native_indices_for_region(
+        16, 10, 2, *arr.shape, hflip=hflip, vflip=vflip
+    )
+    assert np.array_equal(arr, native_rows[:, None] * 1000 + native_cols[None, :])
+
+
+@FLIP_CASES
+def test_region_finer_than_source_resolution(session, hflip, vflip, raster_name):
+    """A region finer than the file's resolution duplicates the nearest cell."""
+    Tools(session=session).g_region(n=16, s=11, w=10, e=15, res=0.5)
+    arr = np.array(garray.array(raster_name, null=NULL, env=session.env))
+    native_rows, native_cols = native_indices_for_region(
+        16, 10, 0.5, *arr.shape, hflip=hflip, vflip=vflip
+    )
+    assert np.array_equal(arr, native_rows[:, None] * 1000 + native_cols[None, :])
+
+
+def test_r_in_gdal_ignores_region(session, source_tif):
+    """r.in.gdal imports the full file regardless of the current region."""
+    tools = Tools(session=session)
+    tools.g_region(n=15, s=8, w=12, e=25, res=1)
+    tools.r_in_gdal(input=str(source_tif), output="imported")
+    tools.g_region(raster="imported")
+    arr = np.array(garray.array("imported", env=session.env))
+    assert arr.shape == (ROWS, COLS)
+    assert np.array_equal(arr, expected_values(0, 0, ROWS, COLS))
+
+
+def test_region_wraps_across_antimeridian(latlon_mapset):
+    """A region crossing the antimeridian reads correctly wrapped columns."""
+    session = latlon_mapset
+    Tools(session=session).g_region(n=5, s=-5, w=170, e=190, res=1)
+    arr = np.array(garray.array("latlon", null=NULL, env=session.env))
+    native_cols = wrapped_native_col_indices(
+        region_west=170, region_east=190, res=1, file_west=-180, file_cols=360
+    )
+    assert np.all(native_cols >= 0)
+    native_rows = np.arange(arr.shape[0])  # n=5, s=-5, res=1 matches the file
+    assert np.array_equal(arr, native_rows[:, None] * 1000 + native_cols[None, :])
+
+
+def test_r_mapcalc_parallel_matches_serial(session):
+    """Results of r.mapcalc with linked data and nprocs = 1 and > 1 match."""
+    tools = Tools(session=session)
+    tools.g_region(n=15, s=8, w=12, e=25, res=1)
+    tools.r_mapcalc(expression="serial = linked", nprocs=1)
+    tools.r_mapcalc(expression="parallel = linked", nprocs=2)
+    serial = np.array(garray.array("serial", env=session.env))
+    parallel = np.array(garray.array("parallel", env=session.env))
+    assert np.array_equal(parallel, expected_values(5, 12, *parallel.shape))
+    assert np.array_equal(parallel, serial)
+
+
+def test_r_univar_parallel_matches_serial(session):
+    """Results of r.univar with  linked data and nprocs = 1 and > 1 match."""
+    tools = Tools(session=session)
+    tools.g_region(n=15, s=8, w=12, e=25, res=1)
+    serial = tools.r_univar(map="linked", nprocs=1, format="json").json
+    parallel = tools.r_univar(map="linked", nprocs=2, format="json").json
+    assert parallel == serial

lib/raster/window_map.c

@@ -106,6 +106,23 @@ void Rast__create_window_mapping(int fd)
        fprintf(stderr, "\n");
      */
 
+    /* For GDAL-linked maps, find the range of logical (pre-flip) columns
+     * needed by the current region. read_data_gdal() (get_row.c) mirrors
+     * this range to the physical columns needed for hflip'ed maps. */
+    fcb->gdal_min_col = -1;
+    fcb->gdal_max_col = -1;
+    if (fcb->gdal) {
+        for (i = 0; i < R__.rd_window.cols; i++) {
+            if (!fcb->col_map[i])
+                continue;
+            if (fcb->gdal_min_col < 0 ||
+                fcb->col_map[i] - 1 < fcb->gdal_min_col)
+                fcb->gdal_min_col = fcb->col_map[i] - 1;
+            if (fcb->col_map[i] - 1 > fcb->gdal_max_col)
+                fcb->gdal_max_col = fcb->col_map[i] - 1;
+        }
+    }
+
     /* compute C1,C2 for row window mapping */
     fcb->C1 = R__.rd_window.ns_res / fcb->cellhd.ns_res;
     fcb->C2 =