From 254a5ecba55238abbecc08f612c596566ccf5f06 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 2 Sep 2026 14:06:14 +0200 Subject: [PATCH 1/8] raster: Restrict GDAL-linked reads to the region's column range --- lib/raster/R.h | 11 ++++++++--- lib/raster/get_row.c | 30 ++++++++++++++++++++++-------- lib/raster/window_map.c | 19 +++++++++++++++++++ 3 files changed, 49 insertions(+), 11 deletions(-) diff --git a/lib/raster/R.h b/lib/raster/R.h index 098ccd0895a..a33838c05bb 100644 --- a/lib/raster/R.h +++ b/lib/raster/R.h @@ -53,9 +53,14 @@ 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 (0-based, inclusive) actually needed by + * the current region for GDAL-linked, non-hflip'ed maps, derived + * from col_map. -1 if not applicable (native full-width read). */ + 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 */ diff --git a/lib/raster/get_row.c b/lib/raster/get_row.c index 5fddb183900..d27b7415679 100644 --- a/lib/raster/get_row.c +++ b/lib/raster/get_row.c @@ -207,26 +207,40 @@ static void read_data_gdal(int fd, int row, unsigned char *data_buf, struct fileinfo *fcb = &R__.fileinfo[fd]; unsigned char *buf; CPLErr err; + /* Restrict the read to the native columns actually needed by the + * region (computed once in Rast__create_window_mapping()) instead + * of always reading the full native row width. This matters a lot + * for maps much wider than the region, such as country-wide + * mosaics built from many source tiles: a full-width read touches + * every one of those source tiles on every row, even if the + * region only overlaps a couple of them. hflip'ed maps keep the + * original full-width behavior, to avoid also having to mirror + * the needed column range. */ + int col_off = 0; + int ncols = fcb->cellhd.cols; + + if (!fcb->gdal->hflip && fcb->gdal_min_col >= 0) { + col_off = fcb->gdal_min_col; + ncols = fcb->gdal_max_col - fcb->gdal_min_col + 1; + } *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++) + for (i = 0; i < ncols; i++) memcpy(data_buf + i * fcb->cur_nbytes, - buf + (fcb->cellhd.cols - 1 - i) * fcb->cur_nbytes, - fcb->cur_nbytes); + buf + (ncols - 1 - i) * fcb->cur_nbytes, fcb->cur_nbytes); G_free(buf); } diff --git a/lib/raster/window_map.c b/lib/raster/window_map.c index 513a00b5270..0a94ca48cae 100644 --- a/lib/raster/window_map.c +++ b/lib/raster/window_map.c @@ -108,6 +108,25 @@ void Rast__create_window_mapping(int fd) fprintf(stderr, "\n"); */ + /* For GDAL-linked, non-hflip'ed maps, find the range of native + * columns actually needed by the current region, so that + * read_data_gdal() can avoid reading (and thus touching every + * source file of, for a mosaic) the full native row width when + * the map is much wider than the region. */ + fcb->gdal_min_col = -1; + fcb->gdal_max_col = -1; + if (fcb->gdal && !fcb->gdal->hflip) { + 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 = From 4515c8c11f198e9b1eb2ba70db782e2f0035743e Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 2 Sep 2026 14:13:43 +0200 Subject: [PATCH 2/8] raster: shorten comments --- lib/raster/get_row.c | 6 +----- lib/raster/window_map.c | 3 +-- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/lib/raster/get_row.c b/lib/raster/get_row.c index d27b7415679..0f0c369be46 100644 --- a/lib/raster/get_row.c +++ b/lib/raster/get_row.c @@ -209,11 +209,7 @@ static void read_data_gdal(int fd, int row, unsigned char *data_buf, CPLErr err; /* Restrict the read to the native columns actually needed by the * region (computed once in Rast__create_window_mapping()) instead - * of always reading the full native row width. This matters a lot - * for maps much wider than the region, such as country-wide - * mosaics built from many source tiles: a full-width read touches - * every one of those source tiles on every row, even if the - * region only overlaps a couple of them. hflip'ed maps keep the + * of always reading the full native row width. hflip'ed maps keep the * original full-width behavior, to avoid also having to mirror * the needed column range. */ int col_off = 0; diff --git a/lib/raster/window_map.c b/lib/raster/window_map.c index 0a94ca48cae..f6867064452 100644 --- a/lib/raster/window_map.c +++ b/lib/raster/window_map.c @@ -110,8 +110,7 @@ void Rast__create_window_mapping(int fd) /* For GDAL-linked, non-hflip'ed maps, find the range of native * columns actually needed by the current region, so that - * read_data_gdal() can avoid reading (and thus touching every - * source file of, for a mosaic) the full native row width when + * read_data_gdal() can avoid reading the full native row width when * the map is much wider than the region. */ fcb->gdal_min_col = -1; fcb->gdal_max_col = -1; From c552f7526a901d73f87209be4f36505e230e261f Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 2 Sep 2026 14:14:37 +0200 Subject: [PATCH 3/8] add testsuite for limiting GDAL linked data reads --- .../tests/lib_raster_gdal_link_window_test.py | 106 ++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 lib/raster/tests/lib_raster_gdal_link_window_test.py diff --git a/lib/raster/tests/lib_raster_gdal_link_window_test.py b/lib/raster/tests/lib_raster_gdal_link_window_test.py new file mode 100644 index 00000000000..927a9443fe4 --- /dev/null +++ b/lib/raster/tests/lib_raster_gdal_link_window_test.py @@ -0,0 +1,106 @@ +"""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. + +These tests generate a raster with a value unique to each cell, link it +back with r.external, and check that reading it through various regions +(fully inside, partially outside, and fully outside the file's extent) +returns the expected values. +""" + +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 + + +@pytest.fixture(scope="module") +def linked_session(tmp_path_factory): + """Module-scoped session with a GeoTIFF exported and linked as 'linked' + + The source raster has a value of (row - 1) * 1000 + (col - 1) at its + 1-based row/col, so a cell's expected value can be computed from its + position 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") + 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 + + +def expected_values(row_offset, col_offset, rows, cols): + """Expected 'linked' values for a region shifted by row/col_offset cells""" + row_values = (row_offset + np.arange(rows)) * 1000 + col_values = col_offset + np.arange(cols) + return row_values[:, None] + col_values[None, :] + + +def test_region_fully_inside_source_extent(session): + """A region fully inside the file reads the correct sub-window""" + # Source extent is n=20, s=0, w=0, e=30. This region's origin is shifted + # by (20 - 15) rows and (12 - 0) columns into the source raster. + Tools(session=session).g_region(n=15, s=8, w=12, e=25, res=1) + arr = np.array(garray.array("linked", null=NULL, env=session.env)) + assert np.array_equal(arr, expected_values(5, 12, *arr.shape)) + + +def test_region_partially_outside_source_extent(session): + """Columns outside the file's extent read as null, the rest as data""" + # w=-5 puts the first 5 columns of the region outside the source + # raster's extent (west=0); columns 5 and up still overlap it. + Tools(session=session).g_region(n=10, s=5, w=-5, e=10, res=1) + arr = np.array(garray.array("linked", 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)) + + +def test_region_fully_outside_source_extent(session): + """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("linked", null=NULL, env=session.env)) + assert np.all(arr == NULL) + + +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)) From 87d3010f83bbe4bcdaf115a87c93ab5bd2b055ba Mon Sep 17 00:00:00 2001 From: ninsbl Date: Thu, 3 Sep 2026 12:59:42 +0200 Subject: [PATCH 4/8] shorten comments --- lib/raster/R.h | 5 ++--- lib/raster/get_row.c | 5 +---- lib/raster/window_map.c | 4 +--- 3 files changed, 4 insertions(+), 10 deletions(-) diff --git a/lib/raster/R.h b/lib/raster/R.h index a33838c05bb..cd61a26164a 100644 --- a/lib/raster/R.h +++ b/lib/raster/R.h @@ -56,9 +56,8 @@ struct fileinfo /* Information for opened cell files */ 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 (0-based, inclusive) actually needed by - * the current region for GDAL-linked, non-hflip'ed maps, derived - * from col_map. -1 if not applicable (native full-width read). */ + /* 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 */ diff --git a/lib/raster/get_row.c b/lib/raster/get_row.c index 0f0c369be46..7f88e9412de 100644 --- a/lib/raster/get_row.c +++ b/lib/raster/get_row.c @@ -208,10 +208,7 @@ static void read_data_gdal(int fd, int row, unsigned char *data_buf, unsigned char *buf; CPLErr err; /* Restrict the read to the native columns actually needed by the - * region (computed once in Rast__create_window_mapping()) instead - * of always reading the full native row width. hflip'ed maps keep the - * original full-width behavior, to avoid also having to mirror - * the needed column range. */ + * region (except for hflip'ed maps). */ int col_off = 0; int ncols = fcb->cellhd.cols; diff --git a/lib/raster/window_map.c b/lib/raster/window_map.c index f6867064452..b9a93f606d0 100644 --- a/lib/raster/window_map.c +++ b/lib/raster/window_map.c @@ -109,9 +109,7 @@ void Rast__create_window_mapping(int fd) */ /* For GDAL-linked, non-hflip'ed maps, find the range of native - * columns actually needed by the current region, so that - * read_data_gdal() can avoid reading the full native row width when - * the map is much wider than the region. */ + * columns needed by the current region. */ fcb->gdal_min_col = -1; fcb->gdal_max_col = -1; if (fcb->gdal && !fcb->gdal->hflip) { From d43814817449c0cebe241e4c7cc19cabbeed7ff8 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Thu, 3 Sep 2026 13:04:54 +0200 Subject: [PATCH 5/8] add tests for OpenMP and latlon cases --- .../tests/lib_raster_gdal_link_window_test.py | 165 +++++++++++++++--- 1 file changed, 144 insertions(+), 21 deletions(-) diff --git a/lib/raster/tests/lib_raster_gdal_link_window_test.py b/lib/raster/tests/lib_raster_gdal_link_window_test.py index 927a9443fe4..6e877ee031d 100644 --- a/lib/raster/tests/lib_raster_gdal_link_window_test.py +++ b/lib/raster/tests/lib_raster_gdal_link_window_test.py @@ -1,13 +1,8 @@ -"""Tests for column-windowed reading of GDAL-linked (r.external) raster maps +"""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. - -These tests generate a raster with a value unique to each cell, link it -back with r.external, and check that reading it through various regions -(fully inside, partially outside, and fully outside the file's extent) -returns the expected values. """ import os @@ -23,15 +18,19 @@ 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): - """Module-scoped session with a GeoTIFF exported and linked as 'linked' + """Session with a GeoTIFF exported and linked as 'linked'. - The source raster has a value of (row - 1) * 1000 + (col - 1) at its - 1-based row/col, so a cell's expected value can be computed from its - position without needing a second, independent read of the file. + 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) @@ -49,7 +48,7 @@ def linked_session(tmp_path_factory): @pytest.fixture def session(linked_session): - """A session in its own temporary mapset, so each test has its own region""" + """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 @@ -57,31 +56,105 @@ def session(linked_session): @pytest.fixture def source_tif(linked_session): - """Path to the GeoTIFF file linked as 'linked' in the session fixture""" + """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 expected_values(row_offset, col_offset, rows, cols): - """Expected 'linked' values for a region shifted by row/col_offset cells""" + """Expected 'linked' values for a region shifted by row/col_offset cells.""" row_values = (row_offset + np.arange(rows)) * 1000 col_values = col_offset + np.arange(cols) return row_values[:, None] + col_values[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): + """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 native_rows, native_cols + + +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 + + def test_region_fully_inside_source_extent(session): - """A region fully inside the file reads the correct sub-window""" - # Source extent is n=20, s=0, w=0, e=30. This region's origin is shifted - # by (20 - 15) rows and (12 - 0) columns into the source raster. + """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("linked", null=NULL, env=session.env)) assert np.array_equal(arr, expected_values(5, 12, *arr.shape)) def test_region_partially_outside_source_extent(session): - """Columns outside the file's extent read as null, the rest as data""" - # w=-5 puts the first 5 columns of the region outside the source - # raster's extent (west=0); columns 5 and up still overlap it. + """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("linked", null=NULL, env=session.env)) assert np.all(arr[:, :5] == NULL) @@ -89,14 +162,30 @@ def test_region_partially_outside_source_extent(session): def test_region_fully_outside_source_extent(session): - """A region with no overlap at all reads back as entirely null""" + """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("linked", null=NULL, env=session.env)) assert np.all(arr == NULL) +def test_region_coarser_than_source_resolution(session): + """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("linked", null=NULL, env=session.env)) + native_rows, native_cols = native_indices_for_region(16, 10, 2, *arr.shape) + assert np.array_equal(arr, native_rows[:, None] * 1000 + native_cols[None, :]) + + +def test_region_finer_than_source_resolution(session): + """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("linked", null=NULL, env=session.env)) + native_rows, native_cols = native_indices_for_region(16, 10, 0.5, *arr.shape) + 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""" + """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") @@ -104,3 +193,37 @@ def test_r_in_gdal_ignores_region(session, source_tif): 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 From 76c2e7273868a5c4c066e979fd0160ca025367d9 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 9 Sep 2026 19:10:32 +0200 Subject: [PATCH 6/8] add tests for h- and v-flipped raster maps --- .../tests/lib_raster_gdal_link_window_test.py | 92 ++++++++++++++----- 1 file changed, 70 insertions(+), 22 deletions(-) diff --git a/lib/raster/tests/lib_raster_gdal_link_window_test.py b/lib/raster/tests/lib_raster_gdal_link_window_test.py index 6e877ee031d..d61895731aa 100644 --- a/lib/raster/tests/lib_raster_gdal_link_window_test.py +++ b/lib/raster/tests/lib_raster_gdal_link_window_test.py @@ -2,7 +2,11 @@ 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. +current region instead of always reading the full native row width. That +column restriction is skipped for maps linked with a horizontal flip, so +the window tests below are parametrized over hflip/vflip, reusing the same +source GeoTIFF linked with r.external's -h/-v flags instead of writing out +an actually mirrored file. """ import os @@ -43,6 +47,10 @@ def linked_session(tmp_path_factory): 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 @@ -89,11 +97,25 @@ def latlon_mapset(latlon_session): yield mapset_session -def expected_values(row_offset, col_offset, rows, cols): - """Expected 'linked' values for a region shifted by row/col_offset cells.""" - row_values = (row_offset + np.arange(rows)) * 1000 - col_values = col_offset + np.arange(cols) - return row_values[:, None] + col_values[None, :] +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): @@ -105,7 +127,7 @@ def nearest_native_index(offset, step, count): return np.floor(offset + step * np.arange(count)).astype(int) -def native_indices_for_region(north, west, res, rows, cols): +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( @@ -114,7 +136,7 @@ def native_indices_for_region(north, west, res, rows, cols): native_rows = nearest_native_index( (FILE_NORTH - north + res / 2.0) / FILE_RES, step, rows ) - return native_rows, native_cols + return apply_flip(native_rows, native_cols, hflip, vflip) def wrapped_native_col_indices(region_west, region_east, res, file_west, file_cols): @@ -146,41 +168,67 @@ def native_for(west): return native -def test_region_fully_inside_source_extent(session): +# 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("linked", null=NULL, env=session.env)) - assert np.array_equal(arr, expected_values(5, 12, *arr.shape)) + 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) + ) -def test_region_partially_outside_source_extent(session): +@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("linked", null=NULL, env=session.env)) + 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)) + assert np.array_equal( + arr[:, 5:], expected_values(10, 0, arr.shape[0], 10, hflip=hflip, vflip=vflip) + ) -def test_region_fully_outside_source_extent(session): +@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("linked", null=NULL, env=session.env)) + arr = np.array(garray.array(raster_name, null=NULL, env=session.env)) assert np.all(arr == NULL) -def test_region_coarser_than_source_resolution(session): +@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("linked", null=NULL, env=session.env)) - native_rows, native_cols = native_indices_for_region(16, 10, 2, *arr.shape) + 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, :]) -def test_region_finer_than_source_resolution(session): +@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("linked", null=NULL, env=session.env)) - native_rows, native_cols = native_indices_for_region(16, 10, 0.5, *arr.shape) + 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, :]) From 0ff63cf3c5e6fb2482d7ec2f0c72ea95f241cf71 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 9 Sep 2026 22:06:48 +0200 Subject: [PATCH 7/8] address code review --- lib/raster/open.c | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lib/raster/open.c b/lib/raster/open.c index bf9b0fb1bb8..f2c025abd41 100644 --- a/lib/raster/open.c +++ b/lib/raster/open.c @@ -40,6 +40,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; } @@ -54,6 +56,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; From dded5ba5055d81702a8d4f7fdb5efa2da4e28bb9 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 9 Sep 2026 22:37:10 +0200 Subject: [PATCH 8/8] implement window limiting for flipped raster maps --- lib/raster/get_row.c | 20 +++++++++---------- .../tests/lib_raster_gdal_link_window_test.py | 6 +----- lib/raster/window_map.c | 7 ++++--- 3 files changed, 15 insertions(+), 18 deletions(-) diff --git a/lib/raster/get_row.c b/lib/raster/get_row.c index 7f88e9412de..c35ef1c5687 100644 --- a/lib/raster/get_row.c +++ b/lib/raster/get_row.c @@ -207,15 +207,15 @@ static void read_data_gdal(int fd, int row, unsigned char *data_buf, struct fileinfo *fcb = &R__.fileinfo[fd]; unsigned char *buf; CPLErr err; - /* Restrict the read to the native columns actually needed by the - * region (except for hflip'ed maps). */ - int col_off = 0; - int ncols = fcb->cellhd.cols; - - if (!fcb->gdal->hflip && fcb->gdal_min_col >= 0) { - col_off = fcb->gdal_min_col; - ncols = fcb->gdal_max_col - fcb->gdal_min_col + 1; - } + /* 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; @@ -232,7 +232,7 @@ static void read_data_gdal(int fd, int row, unsigned char *data_buf, int i; for (i = 0; i < ncols; i++) - memcpy(data_buf + i * fcb->cur_nbytes, + memcpy(data_buf + (min_col + i) * fcb->cur_nbytes, buf + (ncols - 1 - i) * fcb->cur_nbytes, fcb->cur_nbytes); G_free(buf); } diff --git a/lib/raster/tests/lib_raster_gdal_link_window_test.py b/lib/raster/tests/lib_raster_gdal_link_window_test.py index d61895731aa..d3c74041ab6 100644 --- a/lib/raster/tests/lib_raster_gdal_link_window_test.py +++ b/lib/raster/tests/lib_raster_gdal_link_window_test.py @@ -2,11 +2,7 @@ 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. That -column restriction is skipped for maps linked with a horizontal flip, so -the window tests below are parametrized over hflip/vflip, reusing the same -source GeoTIFF linked with r.external's -h/-v flags instead of writing out -an actually mirrored file. +current region instead of always reading the full native row width. """ import os diff --git a/lib/raster/window_map.c b/lib/raster/window_map.c index b9a93f606d0..2d6ac306462 100644 --- a/lib/raster/window_map.c +++ b/lib/raster/window_map.c @@ -108,11 +108,12 @@ void Rast__create_window_mapping(int fd) fprintf(stderr, "\n"); */ - /* For GDAL-linked, non-hflip'ed maps, find the range of native - * columns needed by the current region. */ + /* 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 && !fcb->gdal->hflip) { + if (fcb->gdal) { for (i = 0; i < R__.rd_window.cols; i++) { if (!fcb->col_map[i]) continue;