diff --git a/lib/raster/R.h b/lib/raster/R.h index 098ccd0895a..cd61a26164a 100644 --- a/lib/raster/R.h +++ b/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 */ diff --git a/lib/raster/get_row.c b/lib/raster/get_row.c index a05277ae07a..0f34b23dd91 100644 --- a/lib/raster/get_row.c +++ b/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); } diff --git a/lib/raster/open.c b/lib/raster/open.c index 487e6095eb4..b058e33525f 100644 --- a/lib/raster/open.c +++ b/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; 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..d3c74041ab6 --- /dev/null +++ b/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 diff --git a/lib/raster/window_map.c b/lib/raster/window_map.c index 7ed622e7c42..b888c342cd2 100644 --- a/lib/raster/window_map.c +++ b/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 =