From 2048869867239d574fff25d5c6649adfccf20a44 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 10 Mar 2024 17:45:11 +0100 Subject: [PATCH 1/2] gdal_viewshed: add support for south-up source datasets (fixes #9432) --- alg/viewshed.cpp | 34 ++++--- autotest/utilities/test_gdal_viewshed.py | 118 +++++++++++++++++++++++ 2 files changed, 140 insertions(+), 12 deletions(-) diff --git a/alg/viewshed.cpp b/alg/viewshed.cpp index 20d530d70569..e5cf43f2f12f 100644 --- a/alg/viewshed.cpp +++ b/alg/viewshed.cpp @@ -282,29 +282,39 @@ GDALDatasetH GDALViewshedGenerate( } /* calculate the area of interest */ + constexpr double EPSILON = 1e-8; int nXStart = dfMaxDistance > 0 - ? (std::max)(0, static_cast(std::floor( - nX - adfInvGeoTransform[1] * dfMaxDistance))) + ? (std::max)( + 0, static_cast(std::floor( + nX - adfInvGeoTransform[1] * dfMaxDistance + EPSILON))) : 0; int nXStop = dfMaxDistance > 0 - ? (std::min)(nXSize, - static_cast(std::ceil(nX + adfInvGeoTransform[1] * - dfMaxDistance) + - 1)) + ? (std::min)( + nXSize, + static_cast( + std::ceil(nX + adfInvGeoTransform[1] * dfMaxDistance - + EPSILON) + + 1)) : nXSize; int nYStart = dfMaxDistance > 0 - ? (std::max)(0, static_cast(std::floor( - nY + adfInvGeoTransform[5] * dfMaxDistance))) + ? (std::max)( + 0, static_cast(std::floor( + nY - std::fabs(adfInvGeoTransform[5]) * dfMaxDistance + + EPSILON)) - + (adfInvGeoTransform[5] > 0 ? 1 : 0)) : 0; int nYStop = dfMaxDistance > 0 - ? (std::min)(nYSize, - static_cast(std::ceil(nY - adfInvGeoTransform[5] * - dfMaxDistance) + - 1)) + ? (std::min)( + nYSize, + static_cast(std::ceil(nY + + std::fabs(adfInvGeoTransform[5]) * + dfMaxDistance - + EPSILON) + + (adfInvGeoTransform[5] < 0 ? 1 : 0))) : nYSize; /* normalize horizontal index (0 - nXSize) */ diff --git a/autotest/utilities/test_gdal_viewshed.py b/autotest/utilities/test_gdal_viewshed.py index 8e9a4755df81..01ea6688e1d1 100755 --- a/autotest/utilities/test_gdal_viewshed.py +++ b/autotest/utilities/test_gdal_viewshed.py @@ -29,6 +29,8 @@ # DEALINGS IN THE SOFTWARE. ############################################################################### +import struct + import gdaltest import pytest import test_cli_utilities @@ -304,3 +306,119 @@ def test_gdal_viewshed_invalid_output_filename(gdal_viewshed_path): + " -ox -79.5 -oy 43.5 ../gdrivers/data/n43.tif i/do_not/exist.tif" ) assert "Cannot create dataset" in err + + +############################################################################### +# Test bug fix for https://github.com/OSGeo/gdal/issues/9432 + + +def test_gdal_viewshed_south_up(gdal_viewshed_path, tmp_path, viewshed_input): + + width = 7 + height = 5 + res = 1 + left_x = 1000 + top_y = 2000 + + # "Reference" case with north-up dataset + src_ds_north_up_filename = str(tmp_path / "test_gdal_viewshed_src_ds_north_up.tif") + src_ds_north_up = gdal.GetDriverByName("GTiff").Create( + src_ds_north_up_filename, width, height + ) + src_ds_north_up.SetGeoTransform([left_x, res, 0, top_y, 0, -res]) + expected_gt = src_ds_north_up.GetGeoTransform() + src_ds_north_up.GetRasterBand(1).WriteRaster(width // 2, height // 2, 1, 1, b"\x80") + src_ds_north_up.Close() + + viewshed_out = str(tmp_path / "test_gdal_viewshed_north_up_out.tif") + + _, err = gdaltest.runexternal_out_and_err( + gdal_viewshed_path + + " -oz {} -ox {} -oy {} {} {}".format( + 130, + left_x + float(width) / 2 * res, + top_y - float(height) / 2 * res, + src_ds_north_up_filename, + viewshed_out, + ) + ) + assert err is None or err == "" + ds = gdal.Open(viewshed_out) + assert ds + assert ds.RasterXSize == width + assert ds.RasterYSize == height + assert ds.GetGeoTransform() == pytest.approx(expected_gt) + expected_data = ( + 255, + 255, + 255, + 255, + 255, + 255, + 255, # end of line + 255, + 255, + 0, + 0, + 0, + 255, + 255, # end of line + 255, + 255, + 255, + 255, + 255, + 255, + 255, # end of line + 255, + 255, + 0, + 0, + 0, + 255, + 255, # end of line + 255, + 255, + 255, + 255, + 255, + 255, + 255, + ) + assert ( + struct.unpack("B" * (width * height), ds.GetRasterBand(1).ReadRaster()) + == expected_data + ) + + # Tested case with south-up dataset + src_ds_south_up_filename = str(tmp_path / "test_gdal_viewshed_src_ds_south_up.tif") + src_ds_south_up = gdal.GetDriverByName("GTiff").Create( + src_ds_south_up_filename, width, height + ) + src_ds_south_up.SetGeoTransform([left_x, res, 0, top_y - res * height, 0, res]) + expected_gt = src_ds_south_up.GetGeoTransform() + src_ds_south_up.GetRasterBand(1).WriteRaster(width // 2, height // 2, 1, 1, b"\x80") + src_ds_south_up.Close() + + viewshed_out = str(tmp_path / "test_gdal_viewshed_south_up_out.tif") + + _, err = gdaltest.runexternal_out_and_err( + gdal_viewshed_path + + " -oz {} -ox {} -oy {} {} {}".format( + 130, + left_x + float(width) / 2 * res, + top_y - float(height) / 2 * res, + src_ds_south_up_filename, + viewshed_out, + ) + ) + assert err is None or err == "" + ds = gdal.Open(viewshed_out) + assert ds + assert ds.RasterXSize == width + assert ds.RasterYSize == height + assert ds.GetGeoTransform() == pytest.approx(expected_gt) + assert ( + struct.unpack("B" * (width * height), ds.GetRasterBand(1).ReadRaster()) + == expected_data + ) From a1d0a75e5e1aff25f0c7aab170f924d6d8774e10 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 11 Mar 2024 12:07:38 +0100 Subject: [PATCH 2/2] viewshed.cpp: remove useless parentheses around std::min/max --- alg/viewshed.cpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/alg/viewshed.cpp b/alg/viewshed.cpp index e5cf43f2f12f..234fae858594 100644 --- a/alg/viewshed.cpp +++ b/alg/viewshed.cpp @@ -110,9 +110,9 @@ inline static double CalcHeight(double dfZ, double dfZ2, GDALViewshedMode eMode) if (eMode == GVM_Edge) return dfZ2; else if (eMode == GVM_Max) - return (std::max)(dfZ, dfZ2); + return std::max(dfZ, dfZ2); else if (eMode == GVM_Min) - return (std::min)(dfZ, dfZ2); + return std::min(dfZ, dfZ2); else return dfZ; } @@ -285,13 +285,13 @@ GDALDatasetH GDALViewshedGenerate( constexpr double EPSILON = 1e-8; int nXStart = dfMaxDistance > 0 - ? (std::max)( + ? std::max( 0, static_cast(std::floor( nX - adfInvGeoTransform[1] * dfMaxDistance + EPSILON))) : 0; int nXStop = dfMaxDistance > 0 - ? (std::min)( + ? std::min( nXSize, static_cast( std::ceil(nX + adfInvGeoTransform[1] * dfMaxDistance - @@ -300,7 +300,7 @@ GDALDatasetH GDALViewshedGenerate( : nXSize; int nYStart = dfMaxDistance > 0 - ? (std::max)( + ? std::max( 0, static_cast(std::floor( nY - std::fabs(adfInvGeoTransform[5]) * dfMaxDistance + EPSILON)) - @@ -308,9 +308,8 @@ GDALDatasetH GDALViewshedGenerate( : 0; int nYStop = dfMaxDistance > 0 - ? (std::min)( - nYSize, - static_cast(std::ceil(nY + + ? std::min(nYSize, static_cast( + std::ceil(nY + std::fabs(adfInvGeoTransform[5]) * dfMaxDistance - EPSILON) +