diff --git a/.github/workflows/alpine/Dockerfile.ci b/.github/workflows/alpine/Dockerfile.ci
index 2506ff2739da..05d1c09c09d7 100644
--- a/.github/workflows/alpine/Dockerfile.ci
+++ b/.github/workflows/alpine/Dockerfile.ci
@@ -22,6 +22,7 @@ RUN apk add \
kealib-dev \
libaec-dev \
libarchive-dev \
+ libavif-dev \
libdeflate-dev \
libgeotiff-dev \
libheif-dev \
diff --git a/.github/workflows/alpine_32bit/Dockerfile.ci b/.github/workflows/alpine_32bit/Dockerfile.ci
index bd85ef1320cb..7080b05992bb 100644
--- a/.github/workflows/alpine_32bit/Dockerfile.ci
+++ b/.github/workflows/alpine_32bit/Dockerfile.ci
@@ -24,6 +24,7 @@ RUN apk add \
kealib-dev \
libaec-dev \
libarchive-dev \
+ libavif-dev \
libdeflate-dev \
libgeotiff-dev \
libheif-dev \
diff --git a/.github/workflows/cmake_builds.yml b/.github/workflows/cmake_builds.yml
index c9eb17aaaca9..2eee98263c9d 100644
--- a/.github/workflows/cmake_builds.yml
+++ b/.github/workflows/cmake_builds.yml
@@ -325,7 +325,7 @@ jobs:
mingw-w64-x86_64-geos mingw-w64-x86_64-libspatialite mingw-w64-x86_64-proj
mingw-w64-x86_64-cgal mingw-w64-x86_64-libfreexl mingw-w64-x86_64-hdf5 mingw-w64-x86_64-netcdf mingw-w64-x86_64-poppler mingw-w64-x86_64-podofo mingw-w64-x86_64-postgresql
mingw-w64-x86_64-libgeotiff mingw-w64-x86_64-libpng mingw-w64-x86_64-libtiff mingw-w64-x86_64-openjpeg2
- mingw-w64-x86_64-python-pip mingw-w64-x86_64-python-numpy mingw-w64-x86_64-python-pytest mingw-w64-x86_64-python-setuptools mingw-w64-x86_64-python-lxml mingw-w64-x86_64-swig mingw-w64-x86_64-python-psutil mingw-w64-x86_64-blosc
+ mingw-w64-x86_64-python-pip mingw-w64-x86_64-python-numpy mingw-w64-x86_64-python-pytest mingw-w64-x86_64-python-setuptools mingw-w64-x86_64-python-lxml mingw-w64-x86_64-swig mingw-w64-x86_64-python-psutil mingw-w64-x86_64-blosc mingw-w64-x86_64-libavif
- name: Setup cache
uses: actions/cache@0c45773b623bea8c8e75f6c82b208c3cf94ea4f9 # v4.0.2
id: cache
@@ -429,7 +429,7 @@ jobs:
cfitsio freexl geotiff libjpeg-turbo libpq libspatialite libwebp-base pcre pcre2 postgresql \
sqlite tiledb zstd cryptopp cgal doxygen librttopo libkml openssl xz \
openjdk ant qhull armadillo blas blas-devel libblas libcblas liblapack liblapacke blosc libarchive \
- arrow-cpp pyarrow libaec cmake
+ arrow-cpp pyarrow libaec libavif cmake
- name: Check CMake version
shell: bash -l {0}
run: |
diff --git a/.github/workflows/code_checks.yml b/.github/workflows/code_checks.yml
index 8e09edd06b88..905cbdb2a0ab 100644
--- a/.github/workflows/code_checks.yml
+++ b/.github/workflows/code_checks.yml
@@ -101,6 +101,20 @@ jobs:
# SC2129: (style): Consider using { cmd1; cmd2; } >> file instead of individual redirects
run: shellcheck -e SC2086,SC2046,SC2164,SC2054,SC2129 $(find . -name '*.sh' -a -not -name ltmain.sh -a -not -wholename "./autotest/*" -a -not -wholename "./.github/*")
+ binary_files:
+ runs-on: ubuntu-latest
+ steps:
+
+ - name: Install Requirements
+ run: |
+ sudo apt-get install -y python3 coreutils
+
+ - name: Checkout
+ uses: actions/checkout@692973e3d937129bcbf40652eb9f2f61becf3332 # v4.1.7
+
+ - name: Detect binary files
+ run: python3 ./scripts/check_binaries.py
+
linting:
runs-on: ubuntu-latest
steps:
diff --git a/.github/workflows/fedora_rawhide/Dockerfile.ci b/.github/workflows/fedora_rawhide/Dockerfile.ci
index cb764b75ae00..ca6bb876267b 100644
--- a/.github/workflows/fedora_rawhide/Dockerfile.ci
+++ b/.github/workflows/fedora_rawhide/Dockerfile.ci
@@ -20,6 +20,7 @@ RUN dnf install -y clang make diffutils ccache cmake \
armadillo-devel qhull-devel \
hdf-devel hdf5-devel netcdf-devel \
libpq-devel \
+ libavif-devel \
python3-setuptools python3-pip python3-devel python3-lxml swig \
glibc-gconv-extra
diff --git a/.github/workflows/ubuntu_22.04/Dockerfile.ci b/.github/workflows/ubuntu_22.04/Dockerfile.ci
index d6db8850de46..db60a6e92941 100644
--- a/.github/workflows/ubuntu_22.04/Dockerfile.ci
+++ b/.github/workflows/ubuntu_22.04/Dockerfile.ci
@@ -16,6 +16,7 @@ RUN apt-get update && \
g++ \
git \
gpsbabel \
+ libavif-dev \
libblosc-dev \
libboost-dev \
libcairo2-dev \
diff --git a/.github/workflows/ubuntu_22.04/services.sh b/.github/workflows/ubuntu_22.04/services.sh
index d4d9c002aa2c..19fe9c62404a 100755
--- a/.github/workflows/ubuntu_22.04/services.sh
+++ b/.github/workflows/ubuntu_22.04/services.sh
@@ -7,9 +7,9 @@ set -ex
##################
# MSSQL: server side
-docker rm -f gdal-sql1
-docker pull mcr.microsoft.com/mssql/server:2017-latest
-docker run -e 'ACCEPT_EULA=Y' -e 'SA_PASSWORD=DummyPassw0rd' -p 1433:1433 --name gdal-sql1 -d mcr.microsoft.com/mssql/server:2017-latest
+#docker rm -f gdal-sql1
+#docker pull mcr.microsoft.com/mssql/server:2017-latest
+#docker run -e 'ACCEPT_EULA=Y' -e 'SA_PASSWORD=DummyPassw0rd' -p 1433:1433 --name gdal-sql1 -d mcr.microsoft.com/mssql/server:2017-latest
# MySQL 8
docker rm -f gdal-mysql1
@@ -38,7 +38,7 @@ docker run --name gdal-mongo -p 27018:27017 -d mongo:4.4
sleep 10
# MSSQL
-docker exec -t gdal-sql1 /opt/mssql-tools/bin/sqlcmd -l 30 -S localhost -U SA -P DummyPassw0rd -Q "CREATE DATABASE TestDB"
+#docker exec -t gdal-sql1 /opt/mssql-tools/bin/sqlcmd -l 30 -S localhost -U SA -P DummyPassw0rd -Q "CREATE DATABASE TestDB"
# MySQL
docker exec gdal-mysql1 sh -c "echo 'CREATE DATABASE test; SELECT Version()' | mysql -uroot -ppasswd"
diff --git a/.github/workflows/ubuntu_22.04/test.sh b/.github/workflows/ubuntu_22.04/test.sh
index 9af1dd53f7b7..8bcc8b27ba80 100755
--- a/.github/workflows/ubuntu_22.04/test.sh
+++ b/.github/workflows/ubuntu_22.04/test.sh
@@ -31,6 +31,6 @@ AZURE_STORAGE_CONNECTION_STRING=${AZURITE_STORAGE_CONNECTION_STRING} python3 -c
# MongoDB v3
(cd autotest && MONGODBV3_TEST_PORT=27018 MONGODBV3_TEST_HOST=$IP $PYTEST ogr/ogr_mongodbv3.py)
-(cd autotest && OGR_MSSQL_CONNECTION_STRING="MSSQL:server=$IP;database=TestDB;driver=ODBC Driver 17 for SQL Server;UID=SA;PWD=DummyPassw0rd" $PYTEST ogr/ogr_mssqlspatial.py)
+#(cd autotest && OGR_MSSQL_CONNECTION_STRING="MSSQL:server=$IP;database=TestDB;driver=ODBC Driver 17 for SQL Server;UID=SA;PWD=DummyPassw0rd" $PYTEST ogr/ogr_mssqlspatial.py)
(cd autotest && $PYTEST)
diff --git a/.github/workflows/ubuntu_24.04/Dockerfile.ci b/.github/workflows/ubuntu_24.04/Dockerfile.ci
index 6839a077a1b1..213d1a36069a 100644
--- a/.github/workflows/ubuntu_24.04/Dockerfile.ci
+++ b/.github/workflows/ubuntu_24.04/Dockerfile.ci
@@ -16,6 +16,7 @@ RUN apt-get update && \
g++ \
git \
gpsbabel \
+ libavif-dev \
libblosc-dev \
libboost-dev \
libcairo2-dev \
diff --git a/.github/workflows/ubuntu_24.04/services.sh b/.github/workflows/ubuntu_24.04/services.sh
index d4d9c002aa2c..19fe9c62404a 100755
--- a/.github/workflows/ubuntu_24.04/services.sh
+++ b/.github/workflows/ubuntu_24.04/services.sh
@@ -7,9 +7,9 @@ set -ex
##################
# MSSQL: server side
-docker rm -f gdal-sql1
-docker pull mcr.microsoft.com/mssql/server:2017-latest
-docker run -e 'ACCEPT_EULA=Y' -e 'SA_PASSWORD=DummyPassw0rd' -p 1433:1433 --name gdal-sql1 -d mcr.microsoft.com/mssql/server:2017-latest
+#docker rm -f gdal-sql1
+#docker pull mcr.microsoft.com/mssql/server:2017-latest
+#docker run -e 'ACCEPT_EULA=Y' -e 'SA_PASSWORD=DummyPassw0rd' -p 1433:1433 --name gdal-sql1 -d mcr.microsoft.com/mssql/server:2017-latest
# MySQL 8
docker rm -f gdal-mysql1
@@ -38,7 +38,7 @@ docker run --name gdal-mongo -p 27018:27017 -d mongo:4.4
sleep 10
# MSSQL
-docker exec -t gdal-sql1 /opt/mssql-tools/bin/sqlcmd -l 30 -S localhost -U SA -P DummyPassw0rd -Q "CREATE DATABASE TestDB"
+#docker exec -t gdal-sql1 /opt/mssql-tools/bin/sqlcmd -l 30 -S localhost -U SA -P DummyPassw0rd -Q "CREATE DATABASE TestDB"
# MySQL
docker exec gdal-mysql1 sh -c "echo 'CREATE DATABASE test; SELECT Version()' | mysql -uroot -ppasswd"
diff --git a/.github/workflows/ubuntu_24.04/test.sh b/.github/workflows/ubuntu_24.04/test.sh
index 5ce09486593d..f653a1fdd752 100755
--- a/.github/workflows/ubuntu_24.04/test.sh
+++ b/.github/workflows/ubuntu_24.04/test.sh
@@ -34,6 +34,6 @@ AZURE_STORAGE_CONNECTION_STRING=${AZURITE_STORAGE_CONNECTION_STRING} python3 -c
# MongoDB v3
(cd autotest && MONGODBV3_TEST_PORT=27018 MONGODBV3_TEST_HOST=$IP $PYTEST ogr/ogr_mongodbv3.py)
-(cd autotest && OGR_MSSQL_CONNECTION_STRING="MSSQL:server=$IP;database=TestDB;driver=ODBC Driver 17 for SQL Server;UID=SA;PWD=DummyPassw0rd" $PYTEST ogr/ogr_mssqlspatial.py)
+# (cd autotest && OGR_MSSQL_CONNECTION_STRING="MSSQL:server=$IP;database=TestDB;driver=ODBC Driver 17 for SQL Server;UID=SA;PWD=DummyPassw0rd" $PYTEST ogr/ogr_mssqlspatial.py)
(cd autotest && $PYTEST)
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 56a9e0a98db4..58b2285412a8 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -12,7 +12,7 @@ build:
jobs:
post_checkout:
- - (git --no-pager log --pretty="tformat:%s -- %b" -1 | grep -viqP "skip ci|ci skip") || exit 183
+ - (git --no-pager log --pretty="tformat:%s -- %b" -1 | paste -s -d " " | grep -viqP "skip ci|ci skip") || exit 183
pre_build:
- ./doc/rtd/pre_build.sh
- cd doc && make doxygen generated_rst_files
diff --git a/Doxyfile b/Doxyfile
index 87534d9a9456..4bf006b3802a 100644
--- a/Doxyfile
+++ b/Doxyfile
@@ -382,6 +382,7 @@ INPUT = port \
gcore \
frmts/gdalallregister.cpp \
alg \
+ alg/viewshed \
frmts/vrt \
apps \
ogr \
diff --git a/MIGRATION_GUIDE.TXT b/MIGRATION_GUIDE.TXT
index f1155642d493..146ed58ef8ff 100644
--- a/MIGRATION_GUIDE.TXT
+++ b/MIGRATION_GUIDE.TXT
@@ -18,6 +18,9 @@ MIGRATION GUIDE FROM GDAL 3.9 to GDAL 3.10
- Python bindings: Band.GetStatistics() and Band.ComputeStatistics() now
return a None value in case of error (when exceptions are not enabled)
+- New color interpretation (GCI_xxxx) items have been added to the GDALColorInterp
+ enumeration. Code testing color interpretation may need to be adapted.
+
MIGRATION GUIDE FROM GDAL 3.8 to GDAL 3.9
-----------------------------------------
diff --git a/alg/CMakeLists.txt b/alg/CMakeLists.txt
index ca06ca2fafcc..c4f71ee521a1 100644
--- a/alg/CMakeLists.txt
+++ b/alg/CMakeLists.txt
@@ -36,7 +36,12 @@ add_library(
rasterfill.cpp
thinplatespline.cpp
gdal_simplesurf.cpp
- viewshed.cpp
+ viewshed/combiner.cpp
+ viewshed/cumulative.cpp
+ viewshed/progress.cpp
+ viewshed/util.cpp
+ viewshed/viewshed.cpp
+ viewshed/viewshed_executor.cpp
gdalgenericinverse.cpp
gdal_interpolateatpoint.cpp
)
diff --git a/alg/gdaltransformer.cpp b/alg/gdaltransformer.cpp
index 256a79ed6b8b..3eda4c75b31c 100644
--- a/alg/gdaltransformer.cpp
+++ b/alg/gdaltransformer.cpp
@@ -3828,7 +3828,8 @@ static CPLXMLNode *GDALSerializeApproxTransformer(void *pTransformArg)
* @param pBaseTransformArg the callback argument for the high precision
* transformer.
* @param dfMaxError the maximum cartesian error in the "output" space that
- * is to be accepted in the linear approximation.
+ * is to be accepted in the linear approximation, evaluated as a Manhattan
+ * distance.
*
* @return callback pointer suitable for use with GDALApproxTransform(). It
* should be deallocated with GDALDestroyApproxTransformer().
diff --git a/alg/gdalwarper.cpp b/alg/gdalwarper.cpp
index 079211ae4750..3f09b4b2409a 100644
--- a/alg/gdalwarper.cpp
+++ b/alg/gdalwarper.cpp
@@ -1167,6 +1167,20 @@ CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount,
* will be selected, not just those whose center point falls within the
* polygon.
*
+ *
XSCALE: Ratio expressing the resampling factor (number of destination
+ * pixels per source pixel) along the target horizontal axis.
+ * The scale is used to determine the number of source pixels along the x-axis
+ * that are considered by the resampling algorithm.
+ * Equals to one for no resampling, below one for downsampling
+ * and above one for upsampling. This is automatically computed, for each
+ * processing chunk, and may thus vary among them, depending on the
+ * shape of output regions vs input regions. Such variations can be undesired
+ * in some situations. If the resampling factor can be considered as constant
+ * over the warped area, setting a constant value can lead to more reproducible
+ * pixel output.
+ *
+ * YSCALE: Same as XSCALE, but along the horizontal axis.
+ *
* OPTIMIZE_SIZE: This defaults to FALSE, but may be set to TRUE
* typically when writing to a compressed dataset (GeoTIFF with
* COMPRESS creation option set for example) for achieving a smaller
@@ -1176,7 +1190,11 @@ CPLErr GDALWarpDstAlphaMasker(void *pMaskFuncArg, int nBandCount,
* of the file. However sticking to target block size may cause major
* processing slowdown for some particular reprojections. Starting
* with GDAL 3.8, OPTIMIZE_SIZE mode is automatically enabled when it is safe
- * to do so.
+ * to do so.
+ * As this parameter influences the shape of warping chunk, and by default the
+ * XSCALE and YSCALE parameters are computed per warping chunk, this parameter may
+ * influence the pixel output.
+ *
*
* NUM_THREADS: (GDAL >= 1.10) Can be set to a numeric value or ALL_CPUS to
* set the number of threads to use to parallelize the computation part of the
@@ -1536,6 +1554,32 @@ void CPL_STDCALL GDALWarpResolveWorkingDataType(GDALWarpOptions *psOptions)
psOptions->eWorkingDataType = GDT_Byte;
+ // If none of the provided input nodata values can be represented in the
+ // data type of the corresponding source band, ignore them.
+ if (psOptions->hSrcDS && psOptions->padfSrcNoDataReal)
+ {
+ int nCountInvalidSrcNoDataReal = 0;
+ for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
+ {
+ GDALRasterBandH hSrcBand = GDALGetRasterBand(
+ psOptions->hSrcDS, psOptions->panSrcBands[iBand]);
+
+ if (hSrcBand &&
+ !GDALIsValueExactAs(psOptions->padfSrcNoDataReal[iBand],
+ GDALGetRasterDataType(hSrcBand)))
+ {
+ nCountInvalidSrcNoDataReal++;
+ }
+ }
+ if (nCountInvalidSrcNoDataReal == psOptions->nBandCount)
+ {
+ CPLFree(psOptions->padfSrcNoDataReal);
+ psOptions->padfSrcNoDataReal = nullptr;
+ CPLFree(psOptions->padfSrcNoDataImag);
+ psOptions->padfSrcNoDataImag = nullptr;
+ }
+ }
+
for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
{
if (psOptions->hDstDS != nullptr)
diff --git a/alg/gdalwarpoperation.cpp b/alg/gdalwarpoperation.cpp
index 1448de37f95b..d84d2c08b943 100644
--- a/alg/gdalwarpoperation.cpp
+++ b/alg/gdalwarpoperation.cpp
@@ -832,22 +832,6 @@ void GDALDestroyWarpOperation(GDALWarpOperationH hOperation)
/* CollectChunkList() */
/************************************************************************/
-static int OrderWarpChunk(const void *_a, const void *_b)
-{
- const GDALWarpChunk *a = static_cast(_a);
- const GDALWarpChunk *b = static_cast(_b);
- if (a->dy < b->dy)
- return -1;
- else if (a->dy > b->dy)
- return 1;
- else if (a->dx < b->dx)
- return -1;
- else if (a->dx > b->dx)
- return 1;
- else
- return 0;
-}
-
void GDALWarpOperation::CollectChunkList(int nDstXOff, int nDstYOff,
int nDstXSize, int nDstYSize)
@@ -859,10 +843,18 @@ void GDALWarpOperation::CollectChunkList(int nDstXOff, int nDstYOff,
CollectChunkListInternal(nDstXOff, nDstYOff, nDstXSize, nDstYSize);
// Sort chunks from top to bottom, and for equal y, from left to right.
- // TODO(schwehr): Use std::sort.
- if (pasChunkList)
- qsort(pasChunkList, nChunkListCount, sizeof(GDALWarpChunk),
- OrderWarpChunk);
+ if (nChunkListCount > 1)
+ {
+ std::sort(pasChunkList, pasChunkList + nChunkListCount,
+ [](const GDALWarpChunk &a, const GDALWarpChunk &b)
+ {
+ if (a.dy < b.dy)
+ return true;
+ if (a.dy > b.dy)
+ return false;
+ return a.dx < b.dx;
+ });
+ }
/* -------------------------------------------------------------------- */
/* Find the global source window. */
diff --git a/alg/viewshed.cpp b/alg/viewshed.cpp
deleted file mode 100644
index 6fb84ad6f736..000000000000
--- a/alg/viewshed.cpp
+++ /dev/null
@@ -1,1064 +0,0 @@
-/******************************************************************************
- *
- * Project: Viewshed Generation
- * Purpose: Core algorithm implementation for viewshed generation.
- * Author: Tamas Szekeres, szekerest@gmail.com
- *
- ******************************************************************************
- *
- * Permission is hereby granted, free of charge, to any person obtaining a
- * copy of this software and associated documentation files (the "Software"),
- * to deal in the Software without restriction, including without limitation
- * the rights to use, copy, modify, merge, publish, distribute, sublicense,
- * and/or sell copies of the Software, and to permit persons to whom the
- * Software is furnished to do so, subject to the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
- * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
- * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
- * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
- * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
- * DEALINGS IN THE SOFTWARE.
- ****************************************************************************/
-
-#include
-#include
-#include
-#include
-#include
-
-#include "gdal_alg.h"
-#include "gdal_priv_templates.hpp"
-
-#include "viewshed.h"
-
-/************************************************************************/
-/* GDALViewshedGenerate() */
-/************************************************************************/
-
-/**
- * Create viewshed from raster DEM.
- *
- * This algorithm will generate a viewshed raster from an input DEM raster
- * by using a modified algorithm of "Generating Viewsheds without Using
- * Sightlines" published at
- * https://www.asprs.org/wp-content/uploads/pers/2000journal/january/2000_jan_87-90.pdf
- * This appoach provides a relatively fast calculation, since the output raster
- * is generated in a single scan. The gdal/apps/gdal_viewshed.cpp mainline can
- * be used as an example of how to use this function. The output raster will be
- * of type Byte or Float64.
- *
- * \note The algorithm as implemented currently will only output meaningful
- * results if the georeferencing is in a projected coordinate reference system.
- *
- * @param hBand The band to read the DEM data from. Only the part of the raster
- * within the specified maxdistance around the observer point is processed.
- *
- * @param pszDriverName Driver name (GTiff if set to NULL)
- *
- * @param pszTargetRasterName The name of the target raster to be generated.
- * Must not be NULL
- *
- * @param papszCreationOptions creation options.
- *
- * @param dfObserverX observer X value (in SRS units)
- *
- * @param dfObserverY observer Y value (in SRS units)
- *
- * @param dfObserverHeight The height of the observer above the DEM surface.
- *
- * @param dfTargetHeight The height of the target above the DEM surface.
- * (default 0)
- *
- * @param dfVisibleVal pixel value for visibility (default 255)
- *
- * @param dfInvisibleVal pixel value for invisibility (default 0)
- *
- * @param dfOutOfRangeVal The value to be set for the cells that fall outside of
- * the range specified by dfMaxDistance.
- *
- * @param dfNoDataVal The value to be set for the cells that have no data.
- * If set to a negative value, nodata is not set.
- * Note: currently, no special processing of input cells at a
- * nodata value is done (which may result in erroneous results).
- *
- * @param dfCurvCoeff Coefficient to consider the effect of the curvature and
- * refraction. The height of the DEM is corrected according to the following
- * formula: [Height] -= dfCurvCoeff * [Target Distance]^2 / [Earth Diameter] For
- * the effect of the atmospheric refraction we can use 0.85714.
- *
- * @param eMode The mode of the viewshed calculation.
- * Possible values GVM_Diagonal = 1, GVM_Edge = 2 (default), GVM_Max = 3,
- * GVM_Min = 4.
- *
- * @param dfMaxDistance maximum distance range to compute viewshed.
- * It is also used to clamp the extent of the output
- * raster. If set to 0, then unlimited range is assumed, that is to say the
- * computation is performed on the extent of the whole
- * raster.
- *
- * @param pfnProgress A GDALProgressFunc that may be used to report progress
- * to the user, or to interrupt the algorithm. May be NULL if not required.
- *
- * @param pProgressArg The callback data for the pfnProgress function.
- *
- * @param heightMode Type of information contained in output raster. Possible
- * values GVOT_NORMAL = 1 (default), GVOT_MIN_TARGET_HEIGHT_FROM_DEM = 2,
- * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND = 3
- *
- * GVOT_NORMAL returns a raster of type Byte containing
- * visible locations.
- *
- * GVOT_MIN_TARGET_HEIGHT_FROM_DEM and
- * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND will return a raster of type Float64
- * containing the minimum target height for target to be visible from the DEM
- * surface or ground level respectively. Parameters dfTargetHeight, dfVisibleVal
- * and dfInvisibleVal will be ignored.
- *
- *
- * @param papszExtraOptions Future extra options. Must be set to NULL currently.
- *
- * @return not NULL output dataset on success (to be closed with GDALClose()) or
- * NULL if an error occurs.
- *
- * @since GDAL 3.1
- */
-GDALDatasetH GDALViewshedGenerate(
- GDALRasterBandH hBand, const char *pszDriverName,
- const char *pszTargetRasterName, CSLConstList papszCreationOptions,
- double dfObserverX, double dfObserverY, double dfObserverHeight,
- double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal,
- double dfOutOfRangeVal, double dfNoDataVal, double dfCurvCoeff,
- GDALViewshedMode eMode, double dfMaxDistance, GDALProgressFunc pfnProgress,
- void *pProgressArg, GDALViewshedOutputType heightMode,
- [[maybe_unused]] CSLConstList papszExtraOptions)
-{
- using namespace gdal;
-
- Viewshed::Options oOpts;
- oOpts.outputFormat = pszDriverName;
- oOpts.outputFilename = pszTargetRasterName;
- oOpts.creationOpts = papszCreationOptions;
- oOpts.observer.x = dfObserverX;
- oOpts.observer.y = dfObserverY;
- oOpts.observer.z = dfObserverHeight;
- oOpts.targetHeight = dfTargetHeight;
- oOpts.curveCoeff = dfCurvCoeff;
- oOpts.maxDistance = dfMaxDistance;
- oOpts.nodataVal = dfNoDataVal;
-
- switch (eMode)
- {
- case GVM_Edge:
- oOpts.cellMode = Viewshed::CellMode::Edge;
- break;
- case GVM_Diagonal:
- oOpts.cellMode = Viewshed::CellMode::Diagonal;
- break;
- case GVM_Min:
- oOpts.cellMode = Viewshed::CellMode::Min;
- break;
- case GVM_Max:
- oOpts.cellMode = Viewshed::CellMode::Max;
- break;
- }
-
- switch (heightMode)
- {
- case GVOT_MIN_TARGET_HEIGHT_FROM_DEM:
- oOpts.outputMode = Viewshed::OutputMode::DEM;
- break;
- case GVOT_MIN_TARGET_HEIGHT_FROM_GROUND:
- oOpts.outputMode = Viewshed::OutputMode::Ground;
- break;
- case GVOT_NORMAL:
- oOpts.outputMode = Viewshed::OutputMode::Normal;
- break;
- }
-
- if (!GDALIsValueInRange(dfVisibleVal))
- {
- CPLError(CE_Failure, CPLE_AppDefined,
- "dfVisibleVal out of range. Must be [0, 255].");
- return nullptr;
- }
- if (!GDALIsValueInRange(dfInvisibleVal))
- {
- CPLError(CE_Failure, CPLE_AppDefined,
- "dfInvisibleVal out of range. Must be [0, 255].");
- return nullptr;
- }
- if (!GDALIsValueInRange(dfOutOfRangeVal))
- {
- CPLError(CE_Failure, CPLE_AppDefined,
- "dfOutOfRangeVal out of range. Must be [0, 255].");
- return nullptr;
- }
- oOpts.visibleVal = dfVisibleVal;
- oOpts.invisibleVal = dfInvisibleVal;
- oOpts.outOfRangeVal = dfOutOfRangeVal;
-
- gdal::Viewshed v(oOpts);
-
- if (!pfnProgress)
- pfnProgress = GDALDummyProgress;
- v.run(hBand, pfnProgress, pProgressArg);
-
- return GDALDataset::FromHandle(v.output().release());
-}
-
-namespace gdal
-{
-
-namespace
-{
-
-// Calculate the height adjustment factor.
-double CalcHeightAdjFactor(const GDALDataset *poDataset, double dfCurveCoeff)
-{
- const OGRSpatialReference *poDstSRS = poDataset->GetSpatialRef();
-
- if (poDstSRS)
- {
- OGRErr eSRSerr;
-
- // If we can't get a SemiMajor axis from the SRS, it will be SRS_WGS84_SEMIMAJOR
- double dfSemiMajor = poDstSRS->GetSemiMajor(&eSRSerr);
-
- /* If we fetched the axis from the SRS, use it */
- if (eSRSerr != OGRERR_FAILURE)
- return dfCurveCoeff / (dfSemiMajor * 2.0);
-
- CPLDebug("GDALViewshedGenerate",
- "Unable to fetch SemiMajor axis from spatial reference");
- }
- return 0;
-}
-
-/// Calculate the height at nDistance units along a line through the origin given the height
-/// at nDistance - 1 units along the line.
-/// \param nDistance Distance along the line for the target point.
-/// \param Za Height at the line one unit previous to the target point.
-double CalcHeightLine(int nDistance, double Za)
-{
- nDistance = std::abs(nDistance);
- assert(nDistance != 1);
- return Za * nDistance / (nDistance - 1);
-}
-
-// Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
-// and passing through the line connecting (i - 1, j, Za) and (i, j - 1, Zb).
-// In other words, the origin and the two points form a plane and we're calculating Zc
-// of the point (i, j, Zc), also on the plane.
-double CalcHeightDiagonal(int i, int j, double Za, double Zb)
-{
- return (Za * i + Zb * j) / (i + j - 1);
-}
-
-// Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
-// and through the line connecting (i -1, j - 1, Za) and (i - 1, j, Zb). In other words,
-// the origin and the other two points form a plane and we're calculating Zc of the
-// point (i, j, Zc), also on the plane.
-double CalcHeightEdge(int i, int j, double Za, double Zb)
-{
- assert(i != j);
- return (Za * i + Zb * (j - i)) / (j - 1);
-}
-
-} // unnamed namespace
-
-/// Calculate the extent of the output raster in terms of the input raster.
-///
-/// @param nX observer X position in the input raster
-/// @param nY observer Y position in the input raster
-/// @return false on error, true otherwise
-bool Viewshed::calcOutputExtent(int nX, int nY)
-{
- // We start with the assumption that the output size matches the input.
- oOutExtent.xStop = GDALGetRasterBandXSize(pSrcBand);
- oOutExtent.yStop = GDALGetRasterBandYSize(pSrcBand);
-
- if (!oOutExtent.contains(nX, nY))
- CPLError(CE_Warning, CPLE_AppDefined,
- "NOTE: The observer location falls outside of the DEM area");
-
- constexpr double EPSILON = 1e-8;
- if (oOpts.maxDistance > 0)
- {
- //ABELL - This assumes that the transformation is only a scaling. Should be fixed.
- // Find the distance in the direction of the transformed unit vector in the X and Y
- // directions and use those factors to determine the limiting values in the raster space.
- int nXStart = static_cast(
- std::floor(nX - adfInvTransform[1] * oOpts.maxDistance + EPSILON));
- int nXStop = static_cast(
- std::ceil(nX + adfInvTransform[1] * oOpts.maxDistance - EPSILON) +
- 1);
- int nYStart =
- static_cast(std::floor(
- nY - std::fabs(adfInvTransform[5]) * oOpts.maxDistance +
- EPSILON)) -
- (adfInvTransform[5] > 0 ? 1 : 0);
- int nYStop = static_cast(
- std::ceil(nY + std::fabs(adfInvTransform[5]) * oOpts.maxDistance -
- EPSILON) +
- (adfInvTransform[5] < 0 ? 1 : 0));
-
- // If the limits are invalid, set the window size to zero to trigger the error below.
- if (nXStart >= oOutExtent.xStop || nXStop < 0 ||
- nYStart >= oOutExtent.yStop || nYStop < 0)
- {
- oOutExtent = Window();
- }
- else
- {
- oOutExtent.xStart = std::max(nXStart, 0);
- oOutExtent.xStop = std::min(nXStop, oOutExtent.xStop);
-
- oOutExtent.yStart = std::max(nYStart, 0);
- oOutExtent.yStop = std::min(nYStop, oOutExtent.yStop);
- }
- }
-
- if (oOutExtent.xSize() == 0 || oOutExtent.ySize() == 0)
- {
- CPLError(CE_Failure, CPLE_AppDefined,
- "Invalid target raster size due to transform "
- "and/or distance limitation.");
- return false;
- }
- return true;
-}
-
-/// Read a line of raster data.
-///
-/// @param nLine Line number to read.
-/// @param data Pointer to location in which to store data.
-/// @return Success or failure.
-bool Viewshed::readLine(int nLine, double *data)
-{
- std::lock_guard g(iMutex);
-
- if (GDALRasterIO(pSrcBand, GF_Read, oOutExtent.xStart, nLine,
- oOutExtent.xSize(), 1, data, oOutExtent.xSize(), 1,
- GDT_Float64, 0, 0))
- {
- CPLError(CE_Failure, CPLE_AppDefined,
- "RasterIO error when reading DEM at position (%d,%d), "
- "size (%d,%d)",
- oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
- return false;
- }
- return true;
-}
-
-/// Write an output line of either visibility or height data.
-///
-/// @param nLine Line number being written.
-/// @param vResult Result line to write.
-/// @return True on success, false otherwise.
-bool Viewshed::writeLine(int nLine, std::vector &vResult)
-{
- // GDALRasterIO isn't thread-safe.
- std::lock_guard g(oMutex);
-
- if (GDALRasterIO(pDstBand, GF_Write, 0, nLine - oOutExtent.yStart,
- oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
- 1, GDT_Float64, 0, 0))
- {
- CPLError(CE_Failure, CPLE_AppDefined,
- "RasterIO error when writing target raster at position "
- "(%d,%d), size (%d,%d)",
- 0, nLine - oOutExtent.yStart, oOutExtent.xSize(), 1);
- return false;
- }
- return true;
-}
-
-/// Emit progress information saying that a line has been written to output.
-///
-/// @return True on success, false otherwise.
-bool Viewshed::lineProgress()
-{
- if (nLineCount < oCurExtent.ySize())
- nLineCount++;
- return emitProgress(nLineCount / static_cast(oCurExtent.ySize()));
-}
-
-/// Emit progress information saying that a fraction of work has been completed.
-///
-/// @return True on success, false otherwise.
-bool Viewshed::emitProgress(double fraction)
-{
- // Call the progress function.
- if (!oProgress(fraction, ""))
- {
- CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
- return false;
- }
- return true;
-}
-
-/// Adjust the height of the line of data by the observer height and the curvature of the
-/// earth.
-///
-/// @param nYOffset Y offset of the line being adjusted.
-/// @param nX X location of the observer.
-/// @param vThisLineVal Line height data.
-/// @return [left, right) Leftmost and one past the rightmost cell in the line within
-/// the max distance. Indices are limited to the raster extent (right may be just
-/// outside the raster).
-std::pair Viewshed::adjustHeight(int nYOffset, int nX,
- std::vector &vThisLineVal)
-{
- int nLeft = 0;
- int nRight = oCurExtent.xSize();
-
- // Find the starting point in the raster (nX may be outside)
- int nXStart = oCurExtent.clampX(nX);
-
- // If there is a height adjustment factor other than zero or a max distance,
- // calculate the adjusted height of the cell, stopping if we've exceeded the max
- // distance.
- if (static_cast(dfHeightAdjFactor) || dfMaxDistance2 > 0)
- {
- // Hoist invariants from the loops.
- const double dfLineX = adfTransform[2] * nYOffset;
- const double dfLineY = adfTransform[5] * nYOffset;
-
- // Go left
- double *pdfHeight = vThisLineVal.data() + nXStart;
- for (int nXOffset = nXStart - nX; nXOffset >= -nX;
- nXOffset--, pdfHeight--)
- {
- double dfX = adfTransform[1] * nXOffset + dfLineX;
- double dfY = adfTransform[4] * nXOffset + dfLineY;
- double dfR2 = dfX * dfX + dfY * dfY;
- if (dfR2 > dfMaxDistance2)
- {
- nLeft = nXOffset + nX + 1;
- break;
- }
- *pdfHeight -= dfHeightAdjFactor * dfR2 + dfZObserver;
- }
-
- // Go right.
- pdfHeight = vThisLineVal.data() + nXStart + 1;
- for (int nXOffset = nXStart - nX + 1;
- nXOffset < oCurExtent.xSize() - nX; nXOffset++, pdfHeight++)
- {
- double dfX = adfTransform[1] * nXOffset + dfLineX;
- double dfY = adfTransform[4] * nXOffset + dfLineY;
- double dfR2 = dfX * dfX + dfY * dfY;
- if (dfR2 > dfMaxDistance2)
- {
- nRight = nXOffset + nX;
- break;
- }
- *pdfHeight -= dfHeightAdjFactor * dfR2 + dfZObserver;
- }
- }
- else
- {
- // No curvature adjustment. Just normalize for the observer height.
- double *pdfHeight = vThisLineVal.data();
- for (int i = 0; i < oCurExtent.xSize(); ++i)
- {
- *pdfHeight -= dfZObserver;
- pdfHeight++;
- }
- }
- return {nLeft, nRight};
-}
-
-/// Create the output dataset.
-///
-/// @return True on success, false otherwise.
-bool Viewshed::createOutputDataset()
-{
- GDALDriverManager *hMgr = GetGDALDriverManager();
- GDALDriver *hDriver = hMgr->GetDriverByName(oOpts.outputFormat.c_str());
- if (!hDriver)
- {
- CPLError(CE_Failure, CPLE_AppDefined, "Cannot get driver");
- return false;
- }
-
- /* create output raster */
- poDstDS.reset(hDriver->Create(
- oOpts.outputFilename.c_str(), oOutExtent.xSize(), oOutExtent.ySize(), 1,
- oOpts.outputMode == OutputMode::Normal ? GDT_Byte : GDT_Float64,
- const_cast(oOpts.creationOpts.List())));
- if (!poDstDS)
- {
- CPLError(CE_Failure, CPLE_AppDefined, "Cannot create dataset for %s",
- oOpts.outputFilename.c_str());
- return false;
- }
-
- /* copy srs */
- GDALDatasetH hSrcDS = GDALGetBandDataset(pSrcBand);
- if (hSrcDS)
- poDstDS->SetSpatialRef(
- GDALDataset::FromHandle(hSrcDS)->GetSpatialRef());
-
- std::array adfDstTransform;
- adfDstTransform[0] = adfTransform[0] + adfTransform[1] * oOutExtent.xStart +
- adfTransform[2] * oOutExtent.yStart;
- adfDstTransform[1] = adfTransform[1];
- adfDstTransform[2] = adfTransform[2];
- adfDstTransform[3] = adfTransform[3] + adfTransform[4] * oOutExtent.xStart +
- adfTransform[5] * oOutExtent.yStart;
- adfDstTransform[4] = adfTransform[4];
- adfDstTransform[5] = adfTransform[5];
- poDstDS->SetGeoTransform(adfDstTransform.data());
-
- pDstBand = poDstDS->GetRasterBand(1);
- if (!pDstBand)
- {
- CPLError(CE_Failure, CPLE_AppDefined, "Cannot get band for %s",
- oOpts.outputFilename.c_str());
- return false;
- }
-
- if (oOpts.nodataVal >= 0)
- GDALSetRasterNoDataValue(pDstBand, oOpts.nodataVal);
- return true;
-}
-
-namespace
-{
-
-double doDiagonal(int nXOffset, [[maybe_unused]] int nYOffset,
- double dfThisPrev, double dfLast,
- [[maybe_unused]] double dfLastPrev)
-{
- return CalcHeightDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast);
-}
-
-double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
- double dfLastPrev)
-{
- if (nXOffset >= nYOffset)
- return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
- else
- return CalcHeightEdge(nXOffset, nYOffset, dfLastPrev, dfLast);
-}
-
-double doMin(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
- double dfLastPrev)
-{
- double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
- double dfDiagonal =
- doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
- return std::min(dfEdge, dfDiagonal);
-}
-
-double doMax(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
- double dfLastPrev)
-{
- double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
- double dfDiagonal =
- doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
- return std::max(dfEdge, dfDiagonal);
-}
-
-} // unnamed namespace
-
-/// Process the part of the first line to the left of the observer.
-///
-/// @param nX X coordinate of the observer.
-/// @param iStart X coordinate of the first cell to the left of the observer to be procssed.
-/// @param iEnd X coordinate one past the last cell to be processed.
-/// @param vResult Vector in which to store the visibility/height results.
-/// @param vThisLineVal Height of each cell in the line being processed.
-void Viewshed::processFirstLineLeft(int nX, int iStart, int iEnd,
- std::vector &vResult,
- std::vector &vThisLineVal)
-{
- // If end is to the right of start, everything is taken care of by right processing.
- if (iEnd >= iStart)
- return;
-
- iStart = oCurExtent.clampX(iStart);
-
- double *pThis = vThisLineVal.data() + iStart;
-
- // If the start cell is next to the observer, just mark it visible.
- if (iStart + 1 == nX || iStart + 1 == oCurExtent.xStop)
- {
- if (oOpts.outputMode == OutputMode::Normal)
- vResult[iStart] = oOpts.visibleVal;
- else
- setOutput(vResult[iStart], *pThis, *pThis);
- iStart--;
- pThis--;
- }
-
- // Go from the observer to the left, calculating Z as we go.
- for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--)
- {
- int nXOffset = std::abs(iPixel - nX);
- double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
- setOutput(vResult[iPixel], *pThis, dfZ);
- }
- // For cells outside of the [start, end) range, set the outOfRange value.
- std::fill(vResult.begin(), vResult.begin() + iEnd + 1, oOpts.outOfRangeVal);
-}
-
-/// Process the part of the first line to the right of the observer.
-///
-/// @param nX X coordinate of the observer.
-/// @param iStart X coordinate of the first cell to the right of the observer to be processed.
-/// @param iEnd X coordinate one past the last cell to be processed.
-/// @param vResult Vector in which to store the visibility/height results.
-/// @param vThisLineVal Height of each cell in the line being processed.
-void Viewshed::processFirstLineRight(int nX, int iStart, int iEnd,
- std::vector &vResult,
- std::vector &vThisLineVal)
-{
- // If start is to the right of end, everything is taken care of by left processing.
- if (iStart >= iEnd)
- return;
-
- iStart = oCurExtent.clampX(iStart);
-
- double *pThis = vThisLineVal.data() + iStart;
-
- // If the start cell is next to the observer, just mark it visible.
- if (iStart - 1 == nX || iStart == oCurExtent.xStart)
- {
- if (oOpts.outputMode == OutputMode::Normal)
- vResult[iStart] = oOpts.visibleVal;
- else
- setOutput(vResult[iStart], *pThis, *pThis);
- iStart++;
- pThis++;
- }
-
- // Go from the observer to the right, calculating Z as we go.
- for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++)
- {
- int nXOffset = std::abs(iPixel - nX);
- double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
- setOutput(vResult[iPixel], *pThis, dfZ);
- }
- // For cells outside of the [start, end) range, set the outOfRange value.
- std::fill(vResult.begin() + iEnd, vResult.end(), oOpts.outOfRangeVal);
-}
-
-/// Process a line to the left of the observer.
-///
-/// @param nX X coordinate of the observer.
-/// @param nYOffset Offset of the line being processed from the observer
-/// @param iStart X coordinate of the first cell to the left of the observer to be processed.
-/// @param iEnd X coordinate one past the last cell to be processed.
-/// @param vResult Vector in which to store the visibility/height results.
-/// @param vThisLineVal Height of each cell in the line being processed.
-/// @param vLastLineVal Observable height of each cell in the previous line processed.
-void Viewshed::processLineLeft(int nX, int nYOffset, int iStart, int iEnd,
- std::vector &vResult,
- std::vector &vThisLineVal,
- std::vector &vLastLineVal)
-{
- // If start to the left of end, everything is taken care of by processing right.
- if (iStart <= iEnd)
- return;
- iStart = oCurExtent.clampX(iStart);
-
- nYOffset = std::abs(nYOffset);
- double *pThis = vThisLineVal.data() + iStart;
- double *pLast = vLastLineVal.data() + iStart;
-
- // If the observer is to the right of the raster, mark the first cell to the left as
- // visible. This may mark an out-of-range cell with a value, but this will be fixed
- // with the out of range assignment at the end.
- if (iStart == oCurExtent.xStop - 1)
- {
- if (oOpts.outputMode == OutputMode::Normal)
- vResult[iStart] = oOpts.visibleVal;
- else
- setOutput(vResult[iStart], *pThis, *pThis);
- iStart--;
- pThis--;
- pLast--;
- }
-
- // Go from the observer to the left, calculating Z as we go.
- for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--, pLast--)
- {
- int nXOffset = std::abs(iPixel - nX);
- double dfZ;
- if (nXOffset == nYOffset)
- {
- if (nXOffset == 1)
- dfZ = *pThis;
- else
- dfZ = CalcHeightLine(nXOffset, *(pLast + 1));
- }
- else
- dfZ =
- oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
- setOutput(vResult[iPixel], *pThis, dfZ);
- }
-
- // For cells outside of the [start, end) range, set the outOfRange value.
- std::fill(vResult.begin(), vResult.begin() + iEnd + 1, oOpts.outOfRangeVal);
-}
-
-/// Process a line to the right of the observer.
-///
-/// @param nX X coordinate of the observer.
-/// @param nYOffset Offset of the line being processed from the observer
-/// @param iStart X coordinate of the first cell to the right of the observer to be processed.
-/// @param iEnd X coordinate one past the last cell to be processed.
-/// @param vResult Vector in which to store the visibility/height results.
-/// @param vThisLineVal Height of each cell in the line being processed.
-/// @param vLastLineVal Observable height of each cell in the previous line processed.
-void Viewshed::processLineRight(int nX, int nYOffset, int iStart, int iEnd,
- std::vector &vResult,
- std::vector &vThisLineVal,
- std::vector &vLastLineVal)
-{
- // If start is to the right of end, everything is taken care of by processing left.
- if (iStart >= iEnd)
- return;
- iStart = oCurExtent.clampX(iStart);
-
- nYOffset = std::abs(nYOffset);
- double *pThis = vThisLineVal.data() + iStart;
- double *pLast = vLastLineVal.data() + iStart;
-
- // If the observer is to the left of the raster, mark the first cell to the right as
- // visible. This may mark an out-of-range cell with a value, but this will be fixed
- // with the out of range assignment at the end.
- if (iStart == 0)
- {
- if (oOpts.outputMode == OutputMode::Normal)
- vResult[iStart] = oOpts.visibleVal;
- else
- setOutput(vResult[0], *pThis, *pThis);
- iStart++;
- pThis++;
- pLast++;
- }
-
- // Go from the observer to the right, calculating Z as we go.
- for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++, pLast++)
- {
- int nXOffset = std::abs(iPixel - nX);
- double dfZ;
- if (nXOffset == nYOffset)
- {
- if (nXOffset == 1)
- dfZ = *pThis;
- else
- dfZ = CalcHeightLine(nXOffset, *(pLast - 1));
- }
- else
- dfZ =
- oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
- setOutput(vResult[iPixel], *pThis, dfZ);
- }
-
- // For cells outside of the [start, end) range, set the outOfRange value.
- std::fill(vResult.begin() + iEnd, vResult.end(), oOpts.outOfRangeVal);
-}
-
-/// Set the output Z value depending on the observable height and computation mode.
-///
-/// dfResult Reference to the result cell
-/// dfCellVal Reference to the current cell height. Replace with observable height.
-/// dfZ Minimum observable height at cell.
-void Viewshed::setOutput(double &dfResult, double &dfCellVal, double dfZ)
-{
- if (oOpts.outputMode != OutputMode::Normal)
- {
- dfResult += (dfZ - dfCellVal);
- dfResult = std::max(0.0, dfResult);
- }
- else
- dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
- : oOpts.visibleVal;
- dfCellVal = std::max(dfCellVal, dfZ);
-}
-
-/// Process the first line (the one with the Y coordinate the same as the observer).
-///
-/// @param nX X location of the observer
-/// @param nY Y location of the observer
-/// @param vLastLineVal Vector in which to store the read line. Becomes the last line
-/// in further processing.
-/// @return True on success, false otherwise.
-bool Viewshed::processFirstLine(int nX, int nY,
- std::vector &vLastLineVal)
-{
- int nLine = oOutExtent.clampY(nY);
- int nYOffset = nLine - nY;
-
- std::vector vResult(oOutExtent.xSize());
- std::vector vThisLineVal(oOutExtent.xSize());
-
- if (!readLine(nLine, vThisLineVal.data()))
- return false;
-
- // If the observer is outside of the raster, take the specified value as the Z height,
- // otherwise, take it as an offset from the raster height at that location.
- dfZObserver = oOpts.observer.z;
- if (oCurExtent.containsX(nX))
- {
- dfZObserver += vThisLineVal[nX];
- if (oOpts.outputMode == OutputMode::Normal)
- vResult[nX] = oOpts.visibleVal;
- }
- dfHeightAdjFactor = CalcHeightAdjFactor(poDstDS.get(), oOpts.curveCoeff);
-
- // In DEM mode the base is the pre-adjustment value. In ground mode the base is zero.
- if (oOpts.outputMode == OutputMode::DEM)
- vResult = vThisLineVal;
-
- // iLeft and iRight are the processing limits for the line.
- const auto [iLeft, iRight] = adjustHeight(nYOffset, nX, vThisLineVal);
-
- if (!oCurExtent.containsY(nY))
- processFirstLineTopOrBottom(iLeft, iRight, vResult, vThisLineVal);
- else
- {
- auto t1 = std::async(std::launch::async,
- [&, left = iLeft]() {
- processFirstLineLeft(nX, nX - 1, left - 1,
- vResult, vThisLineVal);
- });
-
- auto t2 = std::async(std::launch::async,
- [&, right = iRight]() {
- processFirstLineRight(nX, nX + 1, right,
- vResult, vThisLineVal);
- });
- t1.wait();
- t2.wait();
- }
-
- // Make the current line the last line.
- vLastLineVal = std::move(vThisLineVal);
-
- // Create the output writer.
- if (!writeLine(nLine, vResult))
- return false;
-
- if (!lineProgress())
- return false;
- return true;
-}
-
-// If the observer is above or below the raster, set all cells in the first line near the
-// observer as observable provided they're in range. Mark cells out of range as such.
-/// @param iLeft Leftmost observable raster position in range of the target line.
-/// @param iRight One past the rightmost observable raster position of the target line.
-/// @param vResult Result line.
-/// @param vThisLineVal Heights of the cells in the target line
-void Viewshed::processFirstLineTopOrBottom(int iLeft, int iRight,
- std::vector &vResult,
- std::vector &vThisLineVal)
-{
- double *pResult = vResult.data() + iLeft;
- double *pThis = vThisLineVal.data() + iLeft;
- for (int iPixel = iLeft; iPixel < iRight; ++iPixel, ++pResult, pThis++)
- {
- if (oOpts.outputMode == OutputMode::Normal)
- *pResult = oOpts.visibleVal;
- else
- setOutput(*pResult, *pThis, *pThis);
- }
- std::fill(vResult.begin(), vResult.begin() + iLeft, oOpts.outOfRangeVal);
- std::fill(vResult.begin() + iRight, vResult.begin() + oCurExtent.xStop,
- oOpts.outOfRangeVal);
-}
-
-/// Process a line above or below the observer.
-///
-/// @param nX X location of the observer
-/// @param nY Y location of the observer
-/// @param nLine Line number being processed.
-/// @param vLastLineVal Vector in which to store the read line. Becomes the last line
-/// in further processing.
-/// @return True on success, false otherwise.
-bool Viewshed::processLine(int nX, int nY, int nLine,
- std::vector &vLastLineVal)
-{
- int nYOffset = nLine - nY;
- std::vector vResult(oOutExtent.xSize());
- std::vector vThisLineVal(oOutExtent.xSize());
-
- if (!readLine(nLine, vThisLineVal.data()))
- return false;
-
- // In DEM mode the base is the input DEM value.
- // In ground mode the base is zero.
- if (oOpts.outputMode == OutputMode::DEM)
- vResult = vThisLineVal;
-
- // Adjust height of the read line.
- const auto [iLeft, iRight] = adjustHeight(nYOffset, nX, vThisLineVal);
-
- // Handle the initial position on the line.
- if (oCurExtent.containsX(nX))
- {
- if (iLeft < iRight)
- {
- double dfZ;
- if (std::abs(nYOffset) == 1)
- dfZ = vThisLineVal[nX];
- else
- dfZ = CalcHeightLine(nYOffset, vLastLineVal[nX]);
- setOutput(vResult[nX], vThisLineVal[nX], dfZ);
- }
- else
- vResult[nX] = oOpts.outOfRangeVal;
- }
-
- // process left half then right half of line
- auto t1 =
- std::async(std::launch::async,
- [&, left = iLeft]()
- {
- processLineLeft(nX, nYOffset, nX - 1, left - 1, vResult,
- vThisLineVal, vLastLineVal);
- });
-
- auto t2 =
- std::async(std::launch::async,
- [&, right = iRight]()
- {
- processLineRight(nX, nYOffset, nX + 1, right, vResult,
- vThisLineVal, vLastLineVal);
- });
- t1.wait();
- t2.wait();
-
- // Make the current line the last line.
- vLastLineVal = std::move(vThisLineVal);
-
- if (!writeLine(nLine, vResult))
- return false;
-
- if (!lineProgress())
- return false;
- return true;
-}
-
-/// Compute the viewshed of a raster band.
-///
-/// @param band Pointer to the raster band to be processed.
-/// @param pfnProgress Pointer to the progress function. Can be null.
-/// @param pProgressArg Argument passed to the progress function
-/// @return True on success, false otherwise.
-bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
- void *pProgressArg)
-{
- using namespace std::placeholders;
-
- nLineCount = 0;
- pSrcBand = static_cast(band);
-
- oProgress = std::bind(pfnProgress, _1, _2, pProgressArg);
-
- if (!emitProgress(0))
- return false;
-
- // set up geotransformation
- GDALDatasetH hSrcDS = GDALGetBandDataset(pSrcBand);
- if (hSrcDS != nullptr)
- GDALGetGeoTransform(hSrcDS, adfTransform.data());
-
- if (!GDALInvGeoTransform(adfTransform.data(), adfInvTransform.data()))
- {
- CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
- return false;
- }
-
- // calculate observer position
- double dfX, dfY;
- GDALApplyGeoTransform(adfInvTransform.data(), oOpts.observer.x,
- oOpts.observer.y, &dfX, &dfY);
- if (!GDALIsValueInRange(dfX))
- {
- CPLError(CE_Failure, CPLE_AppDefined, "Observer X value out of range");
- return false;
- }
- if (!GDALIsValueInRange(dfY))
- {
- CPLError(CE_Failure, CPLE_AppDefined, "Observer Y value out of range");
- return false;
- }
- int nX = static_cast(dfX);
- int nY = static_cast(dfY);
-
- // calculate the area of interest
- if (!calcOutputExtent(nX, nY))
- return false;
-
- // normalize horizontal index to [ 0, oOutExtent.xSize() )
- //ABELL - verify this won't underflow.
- oCurExtent = oOutExtent;
- oCurExtent.shiftX(-oOutExtent.xStart);
- nX -= oOutExtent.xStart;
-
- // create the output dataset
- if (!createOutputDataset())
- return false;
-
- std::vector vFirstLineVal(oCurExtent.xSize());
-
- if (!processFirstLine(nX, nY, vFirstLineVal))
- return false;
-
- if (oOpts.cellMode == CellMode::Edge)
- oZcalc = doEdge;
- else if (oOpts.cellMode == CellMode::Diagonal)
- oZcalc = doDiagonal;
- else if (oOpts.cellMode == CellMode::Min)
- oZcalc = doMin;
- else if (oOpts.cellMode == CellMode::Max)
- oZcalc = doMax;
-
- // scan upwards
- int yStart = oCurExtent.clampY(nY);
- std::atomic err(false);
- auto tUp = std::async(std::launch::async,
- [&]()
- {
- std::vector vLastLineVal = vFirstLineVal;
-
- for (int nLine = yStart - 1;
- nLine >= oCurExtent.yStart && !err; nLine--)
- if (!processLine(nX, nY, nLine, vLastLineVal))
- err = true;
- });
-
- // scan downwards
- auto tDown =
- std::async(std::launch::async,
- [&]()
- {
- std::vector vLastLineVal = vFirstLineVal;
-
- for (int nLine = yStart + 1;
- nLine < oCurExtent.yStop && !err; nLine++)
- if (!processLine(nX, nY, nLine, vLastLineVal))
- err = true;
- });
-
- tUp.wait();
- tDown.wait();
-
- if (!emitProgress(1))
- return false;
-
- return true;
-}
-
-} // namespace gdal
diff --git a/alg/viewshed.h b/alg/viewshed.h
deleted file mode 100644
index 13a1ab558091..000000000000
--- a/alg/viewshed.h
+++ /dev/null
@@ -1,262 +0,0 @@
-/******************************************************************************
- *
- * Project: Viewshed Generation
- * Purpose: Core algorithm implementation for viewshed generation.
- * Author: Tamas Szekeres, szekerest@gmail.com
- *
- ******************************************************************************
- *
- * Permission is hereby granted, free of charge, to any person obtaining a
- * copy of this software and associated documentation files (the "Software"),
- * to deal in the Software without restriction, including without limitation
- * the rights to use, copy, modify, merge, publish, distribute, sublicense,
- * and/or sell copies of the Software, and to permit persons to whom the
- * Software is furnished to do so, subject to the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
- * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
- * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
- * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
- * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
- * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
- * DEALINGS IN THE SOFTWARE.
- ****************************************************************************/
-
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-
-#include "cpl_progress.h"
-#include "gdal_priv.h"
-
-namespace gdal
-{
-
-/**
- * Class to support viewshed raster generation.
- */
-class Viewshed
-{
- public:
- /**
- * Raster output mode.
- */
- enum class OutputMode
- {
- Normal, //!< Normal output mode (visibility only)
- DEM, //!< Output height from DEM
- Ground //!< Output height from ground
- };
-
- /**
- * Cell height calculation mode.
- */
- enum class CellMode
- {
- Diagonal, //!< Diagonal Mode
- Edge, //!< Edge Mode
- Max, //!< Maximum value produced by Diagonal and Edge mode
- Min //!< Minimum value produced by Diagonal and Edge mode
- };
-
- /**
- * A point.
- */
- struct Point
- {
- double x; //!< X value
- double y; //!< Y value
- double z; //!< Z value
- };
-
- /**
- * A window in a raster including pixels in [xStart, xStop) and [yStart, yStop).
- */
- struct Window
- {
- int xStart{}; //!< X start position
- int xStop{}; //!< X end position
- int yStart{}; //!< Y start position
- int yStop{}; //!< Y end position
-
- /// \brief Window size in the X direction.
- int xSize() const
- {
- return xStop - xStart;
- }
-
- /// \brief Window size in the Y direction.
- int ySize() const
- {
- return yStop - yStart;
- }
-
- /// \brief Determine if the X window contains the index.
- /// \param nX Index to check
- /// \return True if the index is contained, false otherwise.
- bool containsX(int nX) const
- {
- return nX >= xStart && nX < xStop;
- }
-
- /// \brief Determine if the Y window contains the index.
- /// \param nY Index to check
- /// \return True if the index is contained, false otherwise.
- bool containsY(int nY) const
- {
- return nY >= xStart && nY < yStop;
- }
-
- /// \brief Determine if the window contains the index.
- /// \param nX X coordinate of the index to check
- /// \param nY Y coordinate of the index to check
- /// \return True if the index is contained, false otherwise.
- bool contains(int nX, int nY) const
- {
- return containsX(nX) && containsY(nY);
- }
-
- /// \brief Clamp the argument to be in the window in the X dimension.
- /// \param nX Value to clamp.
- /// \return Clamped value.
- int clampX(int nX) const
- {
- return xSize() ? std::clamp(nX, xStart, xStop - 1) : xStart;
- }
-
- /// \brief Clamp the argument to be in the window in the Y dimension.
- /// \param nY Value to clamp.
- /// \return Clamped value.
- int clampY(int nY) const
- {
- return ySize() ? std::clamp(nY, yStart, yStop - 1) : yStart;
- }
-
- /// \brief Shift the X dimension by nShift.
- /// \param nShift Amount to shift
- void shiftX(int nShift)
- {
- xStart += nShift;
- xStop += nShift;
- }
- };
-
- /**
- * Options for viewshed generation.
- */
- struct Options
- {
- Point observer{0, 0, 0}; //!< x, y, and z of the observer
- double visibleVal{255}; //!< raster output value for visible pixels.
- double invisibleVal{
- 0}; //!< raster output value for non-visible pixels.
- double outOfRangeVal{
- 0}; //!< raster output value for pixels outside of max distance.
- double nodataVal{-1}; //!< raster output value for pixels with no data
- double targetHeight{0.0}; //!< target height above the DEM surface
- double maxDistance{
- 0.0}; //!< maximum distance from observer to compute value
- double curveCoeff{.85714}; //!< coefficient for atmospheric refraction
- OutputMode outputMode{OutputMode::Normal}; //!< Output information.
- //!< Normal, Height from DEM or Height from ground
- std::string outputFormat{}; //!< output raster format
- std::string outputFilename{}; //!< output raster filename
- CPLStringList creationOpts{}; //!< options for output raster creation
- CellMode cellMode{
- CellMode::Edge}; //!< Mode of cell height calculation.
- };
-
- /**
- * Constructor.
- *
- * @param opts Options to use when calculating viewshed.
- */
- CPL_DLL explicit Viewshed(const Options &opts)
- : oOpts{opts}, oOutExtent{}, oCurExtent{},
- dfMaxDistance2{opts.maxDistance * opts.maxDistance},
- dfZObserver{0}, poDstDS{}, pSrcBand{}, pDstBand{},
- dfHeightAdjFactor{0}, nLineCount{0}, adfTransform{0, 1, 0, 0, 0, 1},
- adfInvTransform{}, oProgress{}, oZcalc{}, oMutex{}, iMutex{}
- {
- if (dfMaxDistance2 == 0)
- dfMaxDistance2 = std::numeric_limits::max();
- }
-
- Viewshed(const Viewshed &) = delete;
- Viewshed &operator=(const Viewshed &) = delete;
-
- CPL_DLL bool run(GDALRasterBandH hBand,
- GDALProgressFunc pfnProgress = GDALDummyProgress,
- void *pProgressArg = nullptr);
-
- /**
- * Fetch a pointer to the created raster band.
- *
- * @return Unique pointer to the viewshed dataset.
- */
- CPL_DLL std::unique_ptr output()
- {
- return std::move(poDstDS);
- }
-
- private:
- Options oOpts;
- Window oOutExtent;
- Window oCurExtent;
- double dfMaxDistance2;
- double dfZObserver;
- std::unique_ptr poDstDS;
- GDALRasterBand *pSrcBand;
- GDALRasterBand *pDstBand;
- double dfHeightAdjFactor;
- int nLineCount;
- std::array adfTransform;
- std::array adfInvTransform;
- using ProgressFunc = std::function;
- ProgressFunc oProgress;
- using ZCalc = std::function;
- ZCalc oZcalc;
- std::mutex oMutex;
- std::mutex iMutex;
-
- void setOutput(double &dfResult, double &dfCellVal, double dfZ);
- double calcHeight(double dfZ, double dfZ2);
- bool readLine(int nLine, double *data);
- bool writeLine(int nLine, std::vector &vResult);
- bool processLine(int nX, int nY, int nLine,
- std::vector &vLastLineVal);
- bool processFirstLine(int nX, int nY, std::vector &vLastLineVal);
- void processFirstLineLeft(int nX, int iStart, int iEnd,
- std::vector &vResult,
- std::vector &vThisLineVal);
- void processFirstLineRight(int nX, int iStart, int iEnd,
- std::vector &vResult,
- std::vector &vThisLineVal);
- void processFirstLineTopOrBottom(int iLeft, int iRight,
- std::vector &vResult,
- std::vector &vThisLineVal);
- void processLineLeft(int nX, int nYOffset, int iStart, int iEnd,
- std::vector &vResult,
- std::vector &vThisLineVal,
- std::vector &vLastLineVal);
- void processLineRight(int nX, int nYOffset, int iStart, int iEnd,
- std::vector &vResult,
- std::vector &vThisLineVal,
- std::vector &vLastLineVal);
- std::pair adjustHeight(int iLine, int nX,
- std::vector &thisLineVal);
- bool calcOutputExtent(int nX, int nY);
- bool createOutputDataset();
- bool lineProgress();
- bool emitProgress(double fraction);
-};
-
-} // namespace gdal
diff --git a/alg/viewshed/combiner.cpp b/alg/viewshed/combiner.cpp
new file mode 100644
index 000000000000..74f8b3f6c7bb
--- /dev/null
+++ b/alg/viewshed/combiner.cpp
@@ -0,0 +1,78 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include "combiner.h"
+#include "util.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/// Read viewshed executor output and sum it up in our owned memory raster.
+void Combiner::run()
+{
+ DatasetPtr pTempDataset;
+
+ while (m_inputQueue.pop(pTempDataset))
+ {
+ if (!m_dataset)
+ m_dataset = std::move(pTempDataset);
+ else
+ sum(std::move(pTempDataset));
+ }
+ // Queue remaining summed rasters.
+ queueOutputBuffer();
+}
+
+/// Add the values of the source dataset to those of the owned dataset.
+/// @param src Source dataset.
+void Combiner::sum(DatasetPtr src)
+{
+ if (!m_dataset)
+ {
+ m_dataset = std::move(src);
+ return;
+ }
+ size_t size = bandSize(*m_dataset->GetRasterBand(1));
+
+ uint8_t *dstP =
+ static_cast(m_dataset->GetInternalHandle("MEMORY1"));
+ uint8_t *srcP = static_cast(src->GetInternalHandle("MEMORY1"));
+ for (size_t i = 0; i < size; ++i)
+ *dstP++ += *srcP++;
+ // If we've seen 255 inputs, queue our raster for output and rollup since we might overflow
+ // otherwise.
+ if (++m_count == 255)
+ queueOutputBuffer();
+}
+
+/// Queue the owned buffer as for output.
+void Combiner::queueOutputBuffer()
+{
+ if (m_dataset)
+ m_outputQueue.push(std::move(m_dataset));
+ m_count = 0;
+}
+
+} // namespace viewshed
+} // namespace gdal
diff --git a/alg/viewshed/combiner.h b/alg/viewshed/combiner.h
new file mode 100644
index 000000000000..0e6ce9787a48
--- /dev/null
+++ b/alg/viewshed/combiner.h
@@ -0,0 +1,72 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#ifndef VIEWSHED_COMBINER_H_INCLUDED
+#define VIEWSHED_COMBINER_H_INCLUDED
+
+#include "cumulative.h"
+#include "viewshed_types.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/// Reads completed viewshed rasters and sums them together. When the
+/// summed values may exceed the 8-bit limit, push it on the output
+/// queue.
+class Combiner
+{
+ public:
+ /// Constructor
+ /// @param inputQueue Reference to input queue of datasets
+ /// @param outputQueue Reference to output queue of datasets
+ Combiner(Cumulative::DatasetQueue &inputQueue,
+ Cumulative::DatasetQueue &outputQueue)
+ : m_inputQueue(inputQueue), m_outputQueue(outputQueue)
+ {
+ }
+
+ /// Copy ctor. Allows initialization in a vector of Combiners.
+ /// @param src Source Combiner.
+ // cppcheck-suppress missingMemberCopy
+ Combiner(const Combiner &src)
+ : m_inputQueue(src.m_inputQueue), m_outputQueue(src.m_outputQueue)
+ {
+ }
+
+ void queueOutputBuffer();
+ void run();
+
+ private:
+ Cumulative::DatasetQueue &m_inputQueue;
+ Cumulative::DatasetQueue &m_outputQueue;
+ DatasetPtr m_dataset{};
+ size_t m_count{0};
+
+ void sum(DatasetPtr srcDs);
+};
+
+} // namespace viewshed
+} // namespace gdal
+
+#endif
diff --git a/alg/viewshed/cumulative.cpp b/alg/viewshed/cumulative.cpp
new file mode 100644
index 000000000000..606ca23c8215
--- /dev/null
+++ b/alg/viewshed/cumulative.cpp
@@ -0,0 +1,233 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include
+#include
+#include
+
+#include "cpl_worker_thread_pool.h"
+
+#include "combiner.h"
+#include "cumulative.h"
+#include "notifyqueue.h"
+#include "util.h"
+#include "viewshed_executor.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/// Constructor
+///
+/// @param opts Options for viewshed generation.
+Cumulative::Cumulative(const Options &opts) : m_opts(opts)
+{
+}
+
+/// Destructor
+///
+Cumulative::~Cumulative() = default;
+
+/// Compute the cumulative viewshed of a raster band.
+///
+/// @param srcFilename Source filename.
+/// @param pfnProgress Pointer to the progress function. Can be null.
+/// @param pProgressArg Argument passed to the progress function
+/// @return True on success, false otherwise.
+bool Cumulative::run(const std::string &srcFilename,
+ GDALProgressFunc pfnProgress, void *pProgressArg)
+{
+ // In cumulative mode, we run the executors in normal mode and want "1" where things
+ // are visible.
+ m_opts.outputMode = OutputMode::Normal;
+ m_opts.visibleVal = 1;
+
+ DatasetPtr srcDS(
+ GDALDataset::FromHandle(GDALOpen(srcFilename.c_str(), GA_ReadOnly)));
+ if (!srcDS)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Unable open source file.");
+ return false;
+ }
+
+ GDALRasterBand *pSrcBand = srcDS->GetRasterBand(1);
+
+ // In cumulative mode, the output extent is always the entire source raster.
+ m_extent.xStop = GDALGetRasterBandXSize(pSrcBand);
+ m_extent.yStop = GDALGetRasterBandYSize(pSrcBand);
+
+ // Make a bunch of observer locations based on the spacing and stick them on a queue
+ // to be handled by viewshed executors.
+ for (int x = 0; x < m_extent.xStop; x += m_opts.observerSpacing)
+ for (int y = 0; y < m_extent.yStop; y += m_opts.observerSpacing)
+ m_observerQueue.push({x, y});
+ m_observerQueue.done();
+
+ // Run executors.
+ const int numThreads = m_opts.numJobs;
+ std::atomic err = false;
+ std::atomic running = numThreads;
+ Progress progress(pfnProgress, pProgressArg,
+ m_observerQueue.size() * m_extent.ySize());
+ CPLWorkerThreadPool executorPool(numThreads);
+ for (int i = 0; i < numThreads; ++i)
+ executorPool.SubmitJob(
+ [this, &srcFilename, &progress, &err, &running]
+ { runExecutor(srcFilename, progress, err, running); });
+
+ // Run combiners that create 8-bit sums of executor jobs.
+ CPLWorkerThreadPool combinerPool(numThreads);
+ std::vector combiners(numThreads,
+ Combiner(m_datasetQueue, m_rollupQueue));
+ for (Combiner &c : combiners)
+ combinerPool.SubmitJob([&c] { c.run(); });
+
+ // Run 32-bit rollup job that combines the 8-bit results from the combiners.
+ std::thread sum([this] { rollupRasters(); });
+
+ // When the combiner jobs are done, all the data is in the rollup queue.
+ combinerPool.WaitCompletion();
+ if (m_datasetQueue.isStopped())
+ return false;
+ m_rollupQueue.done();
+
+ // Wait for finalBuf to be fully filled.
+ sum.join();
+ // The executors should exit naturally, but we wait here so that we don't outrun their
+ // completion and exit with outstanding threads.
+ executorPool.WaitCompletion();
+
+ // Scale the data so that we can write an 8-bit raster output.
+ scaleOutput();
+ if (!writeOutput(createOutputDataset(*pSrcBand, m_opts, m_extent)))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Unable to write to output file.");
+ return false;
+ }
+ progress.emit(1);
+
+ return true;
+}
+
+/// Run an executor (single viewshed)
+/// @param srcFilename Source filename
+/// @param progress Progress supporting support.
+/// @param err Shared error flag.
+/// @param running Shared count of number of executors running.
+void Cumulative::runExecutor(const std::string &srcFilename, Progress &progress,
+ std::atomic &err, std::atomic &running)
+{
+ DatasetPtr srcDs(GDALDataset::Open(srcFilename.c_str(), GA_ReadOnly));
+ if (!srcDs)
+ {
+ err = true;
+ }
+ else
+ {
+ Location loc;
+ while (!err && m_observerQueue.pop(loc))
+ {
+ GDALDriver *memDriver =
+ GetGDALDriverManager()->GetDriverByName("MEM");
+ DatasetPtr dstDs(memDriver ? memDriver->Create("", m_extent.xSize(),
+ m_extent.ySize(), 1,
+ GDT_Byte, nullptr)
+ : nullptr);
+ if (!dstDs)
+ {
+ err = true;
+ }
+ else
+ {
+ ViewshedExecutor executor(
+ *srcDs->GetRasterBand(1), *dstDs->GetRasterBand(1), loc.x,
+ loc.y, m_extent, m_extent, m_opts, progress);
+ err = !executor.run();
+ if (!err)
+ m_datasetQueue.push(std::move(dstDs));
+ }
+ }
+ }
+
+ // Job done. Set the output queue state. If all the executor jobs have completed,
+ // set the dataset output queue done.
+ if (err)
+ m_datasetQueue.stop();
+ else
+ {
+ running--;
+ if (!running)
+ m_datasetQueue.done();
+ }
+}
+
+// Add 8-bit rasters into the 32-bit raster buffer.
+void Cumulative::rollupRasters()
+{
+ DatasetPtr pDS;
+
+ m_finalBuf.resize(m_extent.size());
+ while (m_rollupQueue.pop(pDS))
+ {
+ uint8_t *srcP =
+ static_cast(pDS->GetInternalHandle("MEMORY1"));
+ for (size_t i = 0; i < m_extent.size(); ++i)
+ m_finalBuf[i] += srcP[i];
+ }
+}
+
+/// Scale the output so that it's fully spread in 8 bits. Perhaps this shouldn't happen if
+/// the max is less than 255?
+void Cumulative::scaleOutput()
+{
+ uint32_t m = 0; // This gathers all the bits set.
+ for (uint32_t &val : m_finalBuf)
+ m = std::max(val, m);
+
+ if (m == 0)
+ return;
+
+ double factor =
+ std::numeric_limits::max() / static_cast(m);
+ for (uint32_t &val : m_finalBuf)
+ val = static_cast(std::floor(factor * val));
+}
+
+/// Write the output dataset.
+/// @param pDstDS Pointer to the destination dataset.
+/// @return True if the write was successful, false otherwise.
+bool Cumulative::writeOutput(DatasetPtr pDstDS)
+{
+ if (!pDstDS)
+ return false;
+
+ GDALRasterBand *pDstBand = pDstDS->GetRasterBand(1);
+ return (pDstBand->RasterIO(GF_Write, 0, 0, m_extent.xSize(),
+ m_extent.ySize(), m_finalBuf.data(),
+ m_extent.xSize(), m_extent.ySize(), GDT_UInt32,
+ 0, 0, nullptr) == 0);
+}
+
+} // namespace viewshed
+} // namespace gdal
diff --git a/alg/viewshed/cumulative.h b/alg/viewshed/cumulative.h
new file mode 100644
index 000000000000..c9d2f8a3f958
--- /dev/null
+++ b/alg/viewshed/cumulative.h
@@ -0,0 +1,87 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#ifndef VIEWSHED_CUMULATIVE_H_INCLUDED
+#define VIEWSHED_CUMULATIVE_H_INCLUDED
+
+#include
+#include
+
+#include "notifyqueue.h"
+#include "progress.h"
+#include "viewshed_types.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+class Progress;
+
+/// Generates a cumulative viewshed from a matrix of observers.
+class Cumulative
+{
+ public:
+ CPL_DLL explicit Cumulative(const Options &opts);
+ // We define an explicit destructor, whose implementation is in libgdal,
+ // otherwise with gcc 9.4 of Ubuntu 20.04 in debug mode, this would need to
+ // redefinition of the NotifyQueue class in both libgdal and gdal_viewshed,
+ // leading to weird things related to mutex.
+ CPL_DLL ~Cumulative();
+ CPL_DLL bool run(const std::string &srcFilename,
+ GDALProgressFunc pfnProgress = GDALDummyProgress,
+ void *pProgressArg = nullptr);
+
+ private:
+ friend class Combiner; // Provides access to the queue types.
+
+ struct Location
+ {
+ int x;
+ int y;
+ };
+
+ using Buf32 = std::vector;
+ using ObserverQueue = NotifyQueue;
+ using DatasetQueue = NotifyQueue;
+
+ Window m_extent{};
+ Options m_opts;
+ ObserverQueue m_observerQueue{};
+ DatasetQueue m_datasetQueue{};
+ DatasetQueue m_rollupQueue{};
+ Buf32 m_finalBuf{};
+
+ void runExecutor(const std::string &srcFilename, Progress &progress,
+ std::atomic &err, std::atomic &running);
+ void rollupRasters();
+ void scaleOutput();
+ bool writeOutput(DatasetPtr pDstDS);
+
+ Cumulative(const Cumulative &) = delete;
+ Cumulative &operator=(const Cumulative &) = delete;
+};
+
+} // namespace viewshed
+} // namespace gdal
+
+#endif
diff --git a/alg/viewshed/notifyqueue.h b/alg/viewshed/notifyqueue.h
new file mode 100644
index 000000000000..e0e3f9da3725
--- /dev/null
+++ b/alg/viewshed/notifyqueue.h
@@ -0,0 +1,142 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#ifndef VIEWSHED_NOTIFYQUEUE_H_INCLUDED
+#define VIEWSHED_NOTIFYQUEUE_H_INCLUDED
+
+#include "cpl_port.h"
+
+#include
+#include
+#include
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/// This is a thread-safe queue. Things placed in the queue must be move-constructible.
+/// Readers will wait until there is something in the queue or the queue is empty or stopped.
+/// If the queue is stopped (error), it will never be in the done state. If in the
+/// done state (all writers have finished), it will never be in the error state.
+template class NotifyQueue
+{
+ public:
+ /// Destructor
+ ~NotifyQueue()
+ {
+ done();
+ }
+
+ /// Push an object on the queue and notify readers.
+ /// \param t Object to be moved onto the queue.
+ void push(T &&t)
+ {
+ {
+ std::lock_guard lock(m_mutex);
+ m_queue.push(std::move(t));
+ }
+ m_cv.notify_all();
+ }
+
+ /// Get an item from the queue.
+ /// \param t Reference to an item to to which a queued item will be moved.
+ /// \return True if an item was popped. False otherwise. Use isStopped() or isDone()
+ /// to determine the state if you care when false is returned.
+ bool pop(T &t)
+ {
+ std::unique_lock lock(m_mutex);
+ m_cv.wait(lock,
+ [this] { return !m_queue.empty() || m_done || m_stop; });
+
+ if (m_stop)
+ return false;
+
+ if (m_queue.size())
+ {
+ t = std::move(m_queue.front());
+ m_queue.pop();
+ return true;
+ }
+
+ // m_done must be true and the queue is empty.
+ return false;
+ }
+
+ /// When we're done putting things in the queue, set the end condition.
+ void done()
+ {
+ {
+ std::lock_guard lock(m_mutex);
+ m_done = !m_stop; // If we're already stopped, we can't be done.
+ }
+ m_cv.notify_all();
+ }
+
+ /// Unblock all readers regardless of queue state.
+ void stop()
+ {
+ {
+ std::lock_guard lock(m_mutex);
+ m_stop = !m_done; // If we're already done, we can't be stopped.
+ }
+ m_cv.notify_all();
+ }
+
+ /// Determine if the queue was emptied completely. Call after pop() returns false
+ /// to check queue state.
+ /// \return Whether the queue was emptied completely.
+ bool isDone()
+ {
+ std::lock_guard lock(m_mutex);
+ return m_done;
+ }
+
+ /// Determine if the queue was stopped. Call after pop() returns false
+ /// to check queue state.
+ /// \return Whether the queue was stopped.
+ bool isStopped()
+ {
+ std::lock_guard lock(m_mutex);
+ return m_stop;
+ }
+
+ /// Get the current size of the queue.
+ /// \return Current queue size.
+ size_t size() const
+ {
+ std::lock_guard lock(m_mutex);
+ return m_queue.size();
+ }
+
+ private:
+ std::queue m_queue{};
+ mutable std::mutex m_mutex{};
+ std::condition_variable m_cv{};
+ bool m_done{false};
+ bool m_stop{false};
+};
+
+} // namespace viewshed
+} // namespace gdal
+
+#endif
diff --git a/alg/viewshed/progress.cpp b/alg/viewshed/progress.cpp
new file mode 100644
index 000000000000..3f02b324d9f8
--- /dev/null
+++ b/alg/viewshed/progress.cpp
@@ -0,0 +1,81 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include
+
+#include "progress.h"
+
+#include "cpl_error.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/// Constructor
+/// @param pfnProgress Pointer to progress function.
+/// @param pProgressArg Pointer to progress function data.
+/// @param expectedLines Number of lines expected to be processed.
+Progress::Progress(GDALProgressFunc pfnProgress, void *pProgressArg,
+ size_t expectedLines)
+ : m_expectedLines(std::max(expectedLines, static_cast(1)))
+{
+ using namespace std::placeholders;
+
+ // cppcheck-suppress useInitializationList
+ m_cb = std::bind(pfnProgress, _1, _2, pProgressArg);
+}
+
+/// Emit progress information saying that a line has been written to output.
+///
+/// @return True on success, false otherwise.
+bool Progress::lineComplete()
+{
+ double fraction;
+ {
+ std::lock_guard lock(m_mutex);
+
+ if (m_lines < m_expectedLines)
+ m_lines++;
+ fraction = m_lines / static_cast(m_expectedLines);
+ }
+ return emit(fraction);
+}
+
+/// Emit progress information saying that a fraction of work has been completed.
+///
+/// @return True on success, false otherwise.
+bool Progress::emit(double fraction)
+{
+ std::lock_guard lock(m_mutex);
+
+ // Call the progress function.
+ if (!m_cb(fraction, ""))
+ {
+ CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
+ return false;
+ }
+ return true;
+}
+
+} // namespace viewshed
+} // namespace gdal
diff --git a/alg/viewshed/progress.h b/alg/viewshed/progress.h
new file mode 100644
index 000000000000..0eaee7d449c8
--- /dev/null
+++ b/alg/viewshed/progress.h
@@ -0,0 +1,59 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#ifndef VIEWSHED_PROGRESS_H_INCLUDED
+#define VIEWSHED_PROGRESS_H_INCLUDED
+
+#include
+#include
+
+#include "cpl_progress.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/// Support for progress reporting in viewshed construction. Determines the faction of
+/// progress made based on the number of raster lines completed.
+class Progress
+{
+ public:
+ Progress(GDALProgressFunc pfnProgress, void *pProgressArg,
+ size_t expectedLines);
+
+ bool lineComplete();
+ bool emit(double fraction);
+
+ private:
+ using ProgressFunc = std::function;
+
+ size_t m_lines{0}; ///< Number of lines completed.
+ size_t m_expectedLines; ///< Number of lines expected.
+ std::mutex m_mutex{}; ///< Progress function might not be thread-safe.
+ ProgressFunc m_cb{}; ///< Progress callback function.
+};
+
+} // namespace viewshed
+} // namespace gdal
+
+#endif
diff --git a/alg/viewshed/util.cpp b/alg/viewshed/util.cpp
new file mode 100644
index 000000000000..cf935ed3b3c9
--- /dev/null
+++ b/alg/viewshed/util.cpp
@@ -0,0 +1,104 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include
+
+#include "gdal_priv.h"
+#include "util.h"
+#include "viewshed_types.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/// Get the band size
+///
+/// @param band Raster band
+/// @return The raster band size.
+size_t bandSize(GDALRasterBand &band)
+{
+ return static_cast(band.GetXSize()) * band.GetYSize();
+}
+
+/// Create the output dataset.
+///
+/// @param srcBand Source raster band.
+/// @param opts Options.
+/// @param extent Output dataset extent.
+/// @return The output dataset to be filled with data.
+DatasetPtr createOutputDataset(GDALRasterBand &srcBand, const Options &opts,
+ const Window &extent)
+{
+ GDALDriverManager *hMgr = GetGDALDriverManager();
+ GDALDriver *hDriver = hMgr->GetDriverByName(opts.outputFormat.c_str());
+ if (!hDriver)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Cannot get driver");
+ return nullptr;
+ }
+
+ /* create output raster */
+ DatasetPtr dataset(hDriver->Create(
+ opts.outputFilename.c_str(), extent.xSize(), extent.ySize(), 1,
+ opts.outputMode == OutputMode::Normal ? GDT_Byte : GDT_Float64,
+ const_cast(opts.creationOpts.List())));
+ if (!dataset)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Cannot create dataset for %s",
+ opts.outputFilename.c_str());
+ return nullptr;
+ }
+
+ /* copy srs */
+ dataset->SetSpatialRef(srcBand.GetDataset()->GetSpatialRef());
+
+ std::array adfSrcTransform;
+ std::array adfDstTransform;
+ srcBand.GetDataset()->GetGeoTransform(adfSrcTransform.data());
+ adfDstTransform[0] = adfSrcTransform[0] +
+ adfSrcTransform[1] * extent.xStart +
+ adfSrcTransform[2] * extent.yStart;
+ adfDstTransform[1] = adfSrcTransform[1];
+ adfDstTransform[2] = adfSrcTransform[2];
+ adfDstTransform[3] = adfSrcTransform[3] +
+ adfSrcTransform[4] * extent.xStart +
+ adfSrcTransform[5] * extent.yStart;
+ adfDstTransform[4] = adfSrcTransform[4];
+ adfDstTransform[5] = adfSrcTransform[5];
+ dataset->SetGeoTransform(adfDstTransform.data());
+
+ GDALRasterBand *pBand = dataset->GetRasterBand(1);
+ if (!pBand)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Cannot get band for %s",
+ opts.outputFilename.c_str());
+ return nullptr;
+ }
+
+ if (opts.nodataVal >= 0)
+ GDALSetRasterNoDataValue(pBand, opts.nodataVal);
+ return dataset;
+}
+
+} // namespace viewshed
+} // namespace gdal
diff --git a/alg/viewshed/util.h b/alg/viewshed/util.h
new file mode 100644
index 000000000000..3ebf6a9203dc
--- /dev/null
+++ b/alg/viewshed/util.h
@@ -0,0 +1,41 @@
+/******************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#ifndef VIEWSHED_UTIL_H_INCLUDED
+#define VIEWSHED_UTIL_H_INCLUDED
+
+#include "viewshed_types.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+size_t bandSize(GDALRasterBand &band);
+
+DatasetPtr createOutputDataset(GDALRasterBand &srcBand, const Options &opts,
+ const Window &extent);
+
+} // namespace viewshed
+} // namespace gdal
+
+#endif
diff --git a/alg/viewshed/viewshed.cpp b/alg/viewshed/viewshed.cpp
new file mode 100644
index 000000000000..df5638bb70ac
--- /dev/null
+++ b/alg/viewshed/viewshed.cpp
@@ -0,0 +1,365 @@
+/******************************************************************************
+ *
+ * Project: Viewshed Generation
+ * Purpose: Core algorithm implementation for viewshed generation.
+ * Author: Tamas Szekeres, szekerest@gmail.com
+ *
+ * (c) 2024 info@hobu.co
+ *
+ ******************************************************************************
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include
+
+#include "gdal_alg.h"
+#include "gdal_priv_templates.hpp"
+
+#include "progress.h"
+#include "util.h"
+#include "viewshed.h"
+#include "viewshed_executor.h"
+
+/************************************************************************/
+/* GDALViewshedGenerate() */
+/************************************************************************/
+
+/**
+ * Create viewshed from raster DEM.
+ *
+ * This algorithm will generate a viewshed raster from an input DEM raster
+ * by using a modified algorithm of "Generating Viewsheds without Using
+ * Sightlines" published at
+ * https://www.asprs.org/wp-content/uploads/pers/2000journal/january/2000_jan_87-90.pdf
+ * This approach provides a relatively fast calculation, since the output raster
+ * is generated in a single scan. The gdal/apps/gdal_viewshed.cpp mainline can
+ * be used as an example of how to use this function. The output raster will be
+ * of type Byte or Float64.
+ *
+ * \note The algorithm as implemented currently will only output meaningful
+ * results if the georeferencing is in a projected coordinate reference system.
+ *
+ * @param hBand The band to read the DEM data from. Only the part of the raster
+ * within the specified maxdistance around the observer point is processed.
+ *
+ * @param pszDriverName Driver name (GTiff if set to NULL)
+ *
+ * @param pszTargetRasterName The name of the target raster to be generated.
+ * Must not be NULL
+ *
+ * @param papszCreationOptions creation options.
+ *
+ * @param dfObserverX observer X value (in SRS units)
+ *
+ * @param dfObserverY observer Y value (in SRS units)
+ *
+ * @param dfObserverHeight The height of the observer above the DEM surface.
+ *
+ * @param dfTargetHeight The height of the target above the DEM surface.
+ * (default 0)
+ *
+ * @param dfVisibleVal pixel value for visibility (default 255)
+ *
+ * @param dfInvisibleVal pixel value for invisibility (default 0)
+ *
+ * @param dfOutOfRangeVal The value to be set for the cells that fall outside of
+ * the range specified by dfMaxDistance.
+ *
+ * @param dfNoDataVal The value to be set for the cells that have no data.
+ * If set to a negative value, nodata is not set.
+ * Note: currently, no special processing of input cells at a
+ * nodata value is done (which may result in erroneous results).
+ *
+ * @param dfCurvCoeff Coefficient to consider the effect of the curvature and
+ * refraction. The height of the DEM is corrected according to the following
+ * formula: [Height] -= dfCurvCoeff * [Target Distance]^2 / [Earth Diameter] For
+ * the effect of the atmospheric refraction we can use 0.85714.
+ *
+ * @param eMode The mode of the viewshed calculation.
+ * Possible values GVM_Diagonal = 1, GVM_Edge = 2 (default), GVM_Max = 3,
+ * GVM_Min = 4.
+ *
+ * @param dfMaxDistance maximum distance range to compute viewshed.
+ * It is also used to clamp the extent of the output
+ * raster. If set to 0, then unlimited range is assumed, that is to say the
+ * computation is performed on the extent of the whole
+ * raster.
+ *
+ * @param pfnProgress A GDALProgressFunc that may be used to report progress
+ * to the user, or to interrupt the algorithm. May be NULL if not required.
+ *
+ * @param pProgressArg The callback data for the pfnProgress function.
+ *
+ * @param heightMode Type of information contained in output raster. Possible
+ * values GVOT_NORMAL = 1 (default), GVOT_MIN_TARGET_HEIGHT_FROM_DEM = 2,
+ * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND = 3
+ *
+ * GVOT_NORMAL returns a raster of type Byte containing
+ * visible locations.
+ *
+ * GVOT_MIN_TARGET_HEIGHT_FROM_DEM and
+ * GVOT_MIN_TARGET_HEIGHT_FROM_GROUND will return a raster of type Float64
+ * containing the minimum target height for target to be visible from the DEM
+ * surface or ground level respectively. Parameters dfTargetHeight, dfVisibleVal
+ * and dfInvisibleVal will be ignored.
+ *
+ *
+ * @param papszExtraOptions Future extra options. Must be set to NULL currently.
+ *
+ * @return not NULL output dataset on success (to be closed with GDALClose()) or
+ * NULL if an error occurs.
+ *
+ * @since GDAL 3.1
+ */
+GDALDatasetH GDALViewshedGenerate(
+ GDALRasterBandH hBand, const char *pszDriverName,
+ const char *pszTargetRasterName, CSLConstList papszCreationOptions,
+ double dfObserverX, double dfObserverY, double dfObserverHeight,
+ double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal,
+ double dfOutOfRangeVal, double dfNoDataVal, double dfCurvCoeff,
+ GDALViewshedMode eMode, double dfMaxDistance, GDALProgressFunc pfnProgress,
+ void *pProgressArg, GDALViewshedOutputType heightMode,
+ [[maybe_unused]] CSLConstList papszExtraOptions)
+{
+ using namespace gdal;
+
+ viewshed::Options oOpts;
+ oOpts.outputFormat = pszDriverName;
+ oOpts.outputFilename = pszTargetRasterName;
+ oOpts.creationOpts = papszCreationOptions;
+ oOpts.observer.x = dfObserverX;
+ oOpts.observer.y = dfObserverY;
+ oOpts.observer.z = dfObserverHeight;
+ oOpts.targetHeight = dfTargetHeight;
+ oOpts.curveCoeff = dfCurvCoeff;
+ oOpts.maxDistance = dfMaxDistance;
+ oOpts.nodataVal = dfNoDataVal;
+
+ switch (eMode)
+ {
+ case GVM_Edge:
+ oOpts.cellMode = viewshed::CellMode::Edge;
+ break;
+ case GVM_Diagonal:
+ oOpts.cellMode = viewshed::CellMode::Diagonal;
+ break;
+ case GVM_Min:
+ oOpts.cellMode = viewshed::CellMode::Min;
+ break;
+ case GVM_Max:
+ oOpts.cellMode = viewshed::CellMode::Max;
+ break;
+ }
+
+ switch (heightMode)
+ {
+ case GVOT_MIN_TARGET_HEIGHT_FROM_DEM:
+ oOpts.outputMode = viewshed::OutputMode::DEM;
+ break;
+ case GVOT_MIN_TARGET_HEIGHT_FROM_GROUND:
+ oOpts.outputMode = viewshed::OutputMode::Ground;
+ break;
+ case GVOT_NORMAL:
+ oOpts.outputMode = viewshed::OutputMode::Normal;
+ break;
+ }
+
+ if (!GDALIsValueInRange(dfVisibleVal))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "dfVisibleVal out of range. Must be [0, 255].");
+ return nullptr;
+ }
+ if (!GDALIsValueInRange(dfInvisibleVal))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "dfInvisibleVal out of range. Must be [0, 255].");
+ return nullptr;
+ }
+ if (!GDALIsValueInRange(dfOutOfRangeVal))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "dfOutOfRangeVal out of range. Must be [0, 255].");
+ return nullptr;
+ }
+ oOpts.visibleVal = dfVisibleVal;
+ oOpts.invisibleVal = dfInvisibleVal;
+ oOpts.outOfRangeVal = dfOutOfRangeVal;
+
+ gdal::viewshed::Viewshed v(oOpts);
+
+ if (!pfnProgress)
+ pfnProgress = GDALDummyProgress;
+ v.run(hBand, pfnProgress, pProgressArg);
+
+ return GDALDataset::FromHandle(v.output().release());
+}
+
+namespace gdal
+{
+namespace viewshed
+{
+
+namespace
+{
+
+bool getTransforms(GDALRasterBand &band, double *pFwdTransform,
+ double *pRevTransform)
+{
+ band.GetDataset()->GetGeoTransform(pFwdTransform);
+ if (!GDALInvGeoTransform(pFwdTransform, pRevTransform))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
+ return false;
+ }
+ return true;
+}
+
+} // unnamed namespace
+
+Viewshed::Viewshed(const Options &opts) : oOpts{opts}
+{
+}
+
+Viewshed::~Viewshed() = default;
+
+/// Calculate the extent of the output raster in terms of the input raster and
+/// save the input raster extent.
+///
+/// @return false on error, true otherwise
+bool Viewshed::calcExtents(int nX, int nY,
+ const std::array &adfInvTransform)
+{
+ // We start with the assumption that the output size matches the input.
+ oOutExtent.xStop = GDALGetRasterBandXSize(pSrcBand);
+ oOutExtent.yStop = GDALGetRasterBandYSize(pSrcBand);
+
+ if (!oOutExtent.contains(nX, nY))
+ CPLError(CE_Warning, CPLE_AppDefined,
+ "NOTE: The observer location falls outside of the DEM area");
+
+ constexpr double EPSILON = 1e-8;
+ if (oOpts.maxDistance > 0)
+ {
+ //ABELL - This assumes that the transformation is only a scaling. Should be fixed.
+ // Find the distance in the direction of the transformed unit vector in the X and Y
+ // directions and use those factors to determine the limiting values in the raster space.
+ int nXStart = static_cast(
+ std::floor(nX - adfInvTransform[1] * oOpts.maxDistance + EPSILON));
+ int nXStop = static_cast(
+ std::ceil(nX + adfInvTransform[1] * oOpts.maxDistance - EPSILON) +
+ 1);
+ int nYStart =
+ static_cast(std::floor(
+ nY - std::fabs(adfInvTransform[5]) * oOpts.maxDistance +
+ EPSILON)) -
+ (adfInvTransform[5] > 0 ? 1 : 0);
+ int nYStop = static_cast(
+ std::ceil(nY + std::fabs(adfInvTransform[5]) * oOpts.maxDistance -
+ EPSILON) +
+ (adfInvTransform[5] < 0 ? 1 : 0));
+
+ // If the limits are invalid, set the window size to zero to trigger the error below.
+ if (nXStart >= oOutExtent.xStop || nXStop < 0 ||
+ nYStart >= oOutExtent.yStop || nYStop < 0)
+ {
+ oOutExtent = Window();
+ }
+ else
+ {
+ oOutExtent.xStart = std::max(nXStart, 0);
+ oOutExtent.xStop = std::min(nXStop, oOutExtent.xStop);
+
+ oOutExtent.yStart = std::max(nYStart, 0);
+ oOutExtent.yStop = std::min(nYStop, oOutExtent.yStop);
+ }
+ }
+
+ if (oOutExtent.xSize() == 0 || oOutExtent.ySize() == 0)
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "Invalid target raster size due to transform "
+ "and/or distance limitation.");
+ return false;
+ }
+
+ // normalize horizontal index to [ 0, oOutExtent.xSize() )
+ oCurExtent = oOutExtent;
+ oCurExtent.shiftX(-oOutExtent.xStart);
+
+ return true;
+}
+
+/// Compute the viewshed of a raster band.
+///
+/// @param band Pointer to the raster band to be processed.
+/// @param pfnProgress Pointer to the progress function. Can be null.
+/// @param pProgressArg Argument passed to the progress function
+/// @return True on success, false otherwise.
+bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
+ void *pProgressArg)
+{
+ pSrcBand = static_cast(band);
+
+ std::array adfFwdTransform;
+ std::array adfInvTransform;
+ if (!getTransforms(*pSrcBand, adfFwdTransform.data(),
+ adfInvTransform.data()))
+ return false;
+
+ // calculate observer position
+ double dfX, dfY;
+ GDALApplyGeoTransform(adfInvTransform.data(), oOpts.observer.x,
+ oOpts.observer.y, &dfX, &dfY);
+ if (!GDALIsValueInRange(dfX))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Observer X value out of range");
+ return false;
+ }
+ if (!GDALIsValueInRange(dfY))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined, "Observer Y value out of range");
+ return false;
+ }
+ int nX = static_cast(dfX);
+ int nY = static_cast(dfY);
+
+ // Must calculate extents in order to make the output dataset.
+ if (!calcExtents(nX, nY, adfInvTransform))
+ return false;
+
+ poDstDS = createOutputDataset(*pSrcBand, oOpts, oOutExtent);
+ if (!poDstDS)
+ return false;
+
+ // Create the progress reporter.
+ Progress oProgress(pfnProgress, pProgressArg, oOutExtent.ySize());
+
+ // Execute the viewshed algorithm.
+ GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
+ ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
+ oCurExtent, oOpts, oProgress);
+ executor.run();
+ oProgress.emit(1);
+ return static_cast(poDstDS);
+}
+
+} // namespace viewshed
+} // namespace gdal
diff --git a/alg/viewshed/viewshed.h b/alg/viewshed/viewshed.h
new file mode 100644
index 000000000000..996276596dd9
--- /dev/null
+++ b/alg/viewshed/viewshed.h
@@ -0,0 +1,105 @@
+/******************************************************************************
+ *
+ * Project: Viewshed Generation
+ * Purpose: Core algorithm implementation for viewshed generation.
+ * Author: Tamas Szekeres, szekerest@gmail.com
+ *
+ * (c) 2024 info@hobu.co
+ *
+ ******************************************************************************
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#ifndef VIEWSHED_H_INCLUDED
+#define VIEWSHED_H_INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "cpl_progress.h"
+#include "gdal_priv.h"
+#include "viewshed_types.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/**
+ * Class to support viewshed raster generation.
+ */
+class Viewshed
+{
+ public:
+ /**
+ * Constructor.
+ *
+ * @param opts Options to use when calculating viewshed.
+ */
+ CPL_DLL explicit Viewshed(const Options &opts);
+
+ /** Destructor */
+ CPL_DLL ~Viewshed();
+
+ CPL_DLL bool run(GDALRasterBandH hBand,
+ GDALProgressFunc pfnProgress = GDALDummyProgress,
+ void *pProgressArg = nullptr);
+
+ /**
+ * Fetch a pointer to the created raster band.
+ *
+ * @return Unique pointer to the viewshed dataset.
+ */
+ CPL_DLL DatasetPtr output()
+ {
+ return std::move(poDstDS);
+ }
+
+ private:
+ Options oOpts;
+ Window oOutExtent{};
+ Window oCurExtent{};
+ DatasetPtr poDstDS{};
+ GDALRasterBand *pSrcBand = nullptr;
+
+ DatasetPtr execute(int nX, int nY, const std::string &outFilename);
+ void setOutput(double &dfResult, double &dfCellVal, double dfZ);
+ double calcHeight(double dfZ, double dfZ2);
+ bool readLine(int nLine, double *data);
+ std::pair adjustHeight(int iLine, int nX,
+ std::vector &thisLineVal);
+ bool calcExtents(int nX, int nY,
+ const std::array &adfInvTransform);
+
+ Viewshed(const Viewshed &) = delete;
+ Viewshed &operator=(const Viewshed &) = delete;
+};
+
+} // namespace viewshed
+} // namespace gdal
+
+#endif
diff --git a/alg/viewshed/viewshed_executor.cpp b/alg/viewshed/viewshed_executor.cpp
new file mode 100644
index 000000000000..83c9669b9346
--- /dev/null
+++ b/alg/viewshed/viewshed_executor.cpp
@@ -0,0 +1,687 @@
+/******************************************************************************
+ *
+ * Project: Viewshed Generation
+ * Purpose: Core algorithm implementation for viewshed generation.
+ * Author: Tamas Szekeres, szekerest@gmail.com
+ *
+ * (c) 2024 info@hobu.co
+ *
+ ******************************************************************************
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include
+#include
+#include
+#include
+
+#include "viewshed_executor.h"
+#include "progress.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+namespace
+{
+
+/// Calculate the height at nDistance units along a line through the origin given the height
+/// at nDistance - 1 units along the line.
+/// \param nDistance Distance along the line for the target point.
+/// \param Za Height at the line one unit previous to the target point.
+double CalcHeightLine(int nDistance, double Za)
+{
+ nDistance = std::abs(nDistance);
+ assert(nDistance != 1);
+ return Za * nDistance / (nDistance - 1);
+}
+
+// Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
+// and passing through the line connecting (i - 1, j, Za) and (i, j - 1, Zb).
+// In other words, the origin and the two points form a plane and we're calculating Zc
+// of the point (i, j, Zc), also on the plane.
+double CalcHeightDiagonal(int i, int j, double Za, double Zb)
+{
+ return (Za * i + Zb * j) / (i + j - 1);
+}
+
+// Calculate the height Zc of a point (i, j, Zc) given a line through the origin (0, 0, 0)
+// and through the line connecting (i -1, j - 1, Za) and (i - 1, j, Zb). In other words,
+// the origin and the other two points form a plane and we're calculating Zc of the
+// point (i, j, Zc), also on the plane.
+double CalcHeightEdge(int i, int j, double Za, double Zb)
+{
+ assert(i != j);
+ return (Za * i + Zb * (j - i)) / (j - 1);
+}
+
+double doDiagonal(int nXOffset, [[maybe_unused]] int nYOffset,
+ double dfThisPrev, double dfLast,
+ [[maybe_unused]] double dfLastPrev)
+{
+ return CalcHeightDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast);
+}
+
+double doEdge(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
+ double dfLastPrev)
+{
+ if (nXOffset >= nYOffset)
+ return CalcHeightEdge(nYOffset, nXOffset, dfLastPrev, dfThisPrev);
+ else
+ return CalcHeightEdge(nXOffset, nYOffset, dfLastPrev, dfLast);
+}
+
+double doMin(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
+ double dfLastPrev)
+{
+ double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
+ double dfDiagonal =
+ doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
+ return std::min(dfEdge, dfDiagonal);
+}
+
+double doMax(int nXOffset, int nYOffset, double dfThisPrev, double dfLast,
+ double dfLastPrev)
+{
+ double dfEdge = doEdge(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
+ double dfDiagonal =
+ doDiagonal(nXOffset, nYOffset, dfThisPrev, dfLast, dfLastPrev);
+ return std::max(dfEdge, dfDiagonal);
+}
+
+} // unnamed namespace
+
+/// Constructor -- the viewshed algorithm executor
+/// @param srcBand Source raster band
+/// @param dstBand Destination raster band
+/// @param nX X position of observer
+/// @param nY Y position of observer
+/// @param outExtent Extent of output raster (relative to input)
+/// @param curExtent Extent of active raster.
+/// @param opts Configuration options.
+/// @param progress Reference to the progress tracker.
+ViewshedExecutor::ViewshedExecutor(GDALRasterBand &srcBand,
+ GDALRasterBand &dstBand, int nX, int nY,
+ const Window &outExtent,
+ const Window &curExtent, const Options &opts,
+ Progress &progress)
+ : m_pool(4), m_srcBand(srcBand), m_dstBand(dstBand), oOutExtent(outExtent),
+ oCurExtent(curExtent), m_nX(nX - oOutExtent.xStart), m_nY(nY),
+ oOpts(opts), oProgress(progress),
+ m_dfMaxDistance2(opts.maxDistance * opts.maxDistance)
+{
+ if (m_dfMaxDistance2 == 0)
+ m_dfMaxDistance2 = std::numeric_limits::max();
+ m_srcBand.GetDataset()->GetGeoTransform(m_adfTransform.data());
+}
+
+// calculate the height adjustment factor.
+double ViewshedExecutor::calcHeightAdjFactor()
+{
+ std::lock_guard g(oMutex);
+
+ const OGRSpatialReference *poDstSRS =
+ m_dstBand.GetDataset()->GetSpatialRef();
+
+ if (poDstSRS)
+ {
+ OGRErr eSRSerr;
+
+ // If we can't get a SemiMajor axis from the SRS, it will be SRS_WGS84_SEMIMAJOR
+ double dfSemiMajor = poDstSRS->GetSemiMajor(&eSRSerr);
+
+ /* If we fetched the axis from the SRS, use it */
+ if (eSRSerr != OGRERR_FAILURE)
+ return oOpts.curveCoeff / (dfSemiMajor * 2.0);
+
+ CPLDebug("GDALViewshedGenerate",
+ "Unable to fetch SemiMajor axis from spatial reference");
+ }
+ return 0;
+}
+
+/// Set the output Z value depending on the observable height and computation mode.
+///
+/// dfResult Reference to the result cell
+/// dfCellVal Reference to the current cell height. Replace with observable height.
+/// dfZ Minimum observable height at cell.
+void ViewshedExecutor::setOutput(double &dfResult, double &dfCellVal,
+ double dfZ)
+{
+ if (oOpts.outputMode != OutputMode::Normal)
+ {
+ dfResult += (dfZ - dfCellVal);
+ dfResult = std::max(0.0, dfResult);
+ }
+ else
+ dfResult = (dfCellVal + oOpts.targetHeight < dfZ) ? oOpts.invisibleVal
+ : oOpts.visibleVal;
+ dfCellVal = std::max(dfCellVal, dfZ);
+}
+
+/// Read a line of raster data.
+///
+/// @param nLine Line number to read.
+/// @param data Pointer to location in which to store data.
+/// @return Success or failure.
+bool ViewshedExecutor::readLine(int nLine, double *data)
+{
+ std::lock_guard g(iMutex);
+
+ if (GDALRasterIO(&m_srcBand, GF_Read, oOutExtent.xStart, nLine,
+ oOutExtent.xSize(), 1, data, oOutExtent.xSize(), 1,
+ GDT_Float64, 0, 0))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "RasterIO error when reading DEM at position (%d,%d), "
+ "size (%d,%d)",
+ oOutExtent.xStart, nLine, oOutExtent.xSize(), 1);
+ return false;
+ }
+ return true;
+}
+
+/// Write an output line of either visibility or height data.
+///
+/// @param nLine Line number being written.
+/// @param vResult Result line to write.
+/// @return True on success, false otherwise.
+bool ViewshedExecutor::writeLine(int nLine, std::vector &vResult)
+{
+ // GDALRasterIO isn't thread-safe.
+ std::lock_guard g(oMutex);
+
+ if (GDALRasterIO(&m_dstBand, GF_Write, 0, nLine - oOutExtent.yStart,
+ oOutExtent.xSize(), 1, vResult.data(), oOutExtent.xSize(),
+ 1, GDT_Float64, 0, 0))
+ {
+ CPLError(CE_Failure, CPLE_AppDefined,
+ "RasterIO error when writing target raster at position "
+ "(%d,%d), size (%d,%d)",
+ 0, nLine - oOutExtent.yStart, oOutExtent.xSize(), 1);
+ return false;
+ }
+ return true;
+}
+
+/// Adjust the height of the line of data by the observer height and the curvature of the
+/// earth.
+///
+/// @param nYOffset Y offset of the line being adjusted.
+/// @param vThisLineVal Line height data.
+/// @return [left, right) Leftmost and one past the rightmost cell in the line within
+/// the max distance. Indices are limited to the raster extent (right may be just
+/// outside the raster).
+std::pair
+ViewshedExecutor::adjustHeight(int nYOffset, std::vector &vThisLineVal)
+{
+ int nLeft = 0;
+ int nRight = oCurExtent.xSize();
+
+ // Find the starting point in the raster (m_nX may be outside)
+ int nXStart = oCurExtent.clampX(m_nX);
+
+ // If there is a height adjustment factor other than zero or a max distance,
+ // calculate the adjusted height of the cell, stopping if we've exceeded the max
+ // distance.
+ if (static_cast(m_dfHeightAdjFactor) || m_dfMaxDistance2 > 0)
+ {
+ // Hoist invariants from the loops.
+ const double dfLineX = m_adfTransform[2] * nYOffset;
+ const double dfLineY = m_adfTransform[5] * nYOffset;
+
+ // Go left
+ double *pdfHeight = vThisLineVal.data() + nXStart;
+ for (int nXOffset = nXStart - m_nX; nXOffset >= -m_nX;
+ nXOffset--, pdfHeight--)
+ {
+ double dfX = m_adfTransform[1] * nXOffset + dfLineX;
+ double dfY = m_adfTransform[4] * nXOffset + dfLineY;
+ double dfR2 = dfX * dfX + dfY * dfY;
+ if (dfR2 > m_dfMaxDistance2)
+ {
+ nLeft = nXOffset + m_nX + 1;
+ break;
+ }
+ *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
+ }
+
+ // Go right.
+ pdfHeight = vThisLineVal.data() + nXStart + 1;
+ for (int nXOffset = nXStart - m_nX + 1;
+ nXOffset < oCurExtent.xSize() - m_nX; nXOffset++, pdfHeight++)
+ {
+ double dfX = m_adfTransform[1] * nXOffset + dfLineX;
+ double dfY = m_adfTransform[4] * nXOffset + dfLineY;
+ double dfR2 = dfX * dfX + dfY * dfY;
+ if (dfR2 > m_dfMaxDistance2)
+ {
+ nRight = nXOffset + m_nX;
+ break;
+ }
+ *pdfHeight -= m_dfHeightAdjFactor * dfR2 + m_dfZObserver;
+ }
+ }
+ else
+ {
+ // No curvature adjustment. Just normalize for the observer height.
+ double *pdfHeight = vThisLineVal.data();
+ for (int i = 0; i < oCurExtent.xSize(); ++i)
+ {
+ *pdfHeight -= m_dfZObserver;
+ pdfHeight++;
+ }
+ }
+ return {nLeft, nRight};
+}
+
+/// Process the first line (the one with the Y coordinate the same as the observer).
+///
+/// @param vLastLineVal Vector in which to store the read line. Becomes the last line
+/// in further processing.
+/// @return True on success, false otherwise.
+bool ViewshedExecutor::processFirstLine(std::vector &vLastLineVal)
+{
+ int nLine = oOutExtent.clampY(m_nY);
+ int nYOffset = nLine - m_nY;
+
+ std::vector vResult(oOutExtent.xSize());
+ std::vector vThisLineVal(oOutExtent.xSize());
+
+ if (!readLine(nLine, vThisLineVal.data()))
+ return false;
+
+ // If the observer is outside of the raster, take the specified value as the Z height,
+ // otherwise, take it as an offset from the raster height at that location.
+ m_dfZObserver = oOpts.observer.z;
+ if (oCurExtent.containsX(m_nX))
+ {
+ m_dfZObserver += vThisLineVal[m_nX];
+ if (oOpts.outputMode == OutputMode::Normal)
+ vResult[m_nX] = oOpts.visibleVal;
+ }
+ m_dfHeightAdjFactor = calcHeightAdjFactor();
+
+ // In DEM mode the base is the pre-adjustment value. In ground mode the base is zero.
+ if (oOpts.outputMode == OutputMode::DEM)
+ vResult = vThisLineVal;
+
+ // iLeft and iRight are the processing limits for the line.
+ const auto [iLeft, iRight] = adjustHeight(nYOffset, vThisLineVal);
+
+ if (!oCurExtent.containsY(m_nY))
+ processFirstLineTopOrBottom(iLeft, iRight, vResult, vThisLineVal);
+ else
+ {
+ CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
+ pQueue->SubmitJob(
+ [&, left = iLeft]() {
+ processFirstLineLeft(m_nX - 1, left - 1, vResult, vThisLineVal);
+ });
+
+ pQueue->SubmitJob(
+ [&, right = iRight]()
+ { processFirstLineRight(m_nX + 1, right, vResult, vThisLineVal); });
+ pQueue->WaitCompletion();
+ }
+
+ // Make the current line the last line.
+ vLastLineVal = std::move(vThisLineVal);
+
+ // Create the output writer.
+ if (!writeLine(nLine, vResult))
+ return false;
+
+ return oProgress.lineComplete();
+}
+
+// If the observer is above or below the raster, set all cells in the first line near the
+// observer as observable provided they're in range. Mark cells out of range as such.
+/// @param iLeft Leftmost observable raster position in range of the target line.
+/// @param iRight One past the rightmost observable raster position of the target line.
+/// @param vResult Result line.
+/// @param vThisLineVal Heights of the cells in the target line
+void ViewshedExecutor::processFirstLineTopOrBottom(
+ int iLeft, int iRight, std::vector &vResult,
+ std::vector &vThisLineVal)
+{
+ double *pResult = vResult.data() + iLeft;
+ double *pThis = vThisLineVal.data() + iLeft;
+ for (int iPixel = iLeft; iPixel < iRight; ++iPixel, ++pResult, pThis++)
+ {
+ if (oOpts.outputMode == OutputMode::Normal)
+ *pResult = oOpts.visibleVal;
+ else
+ setOutput(*pResult, *pThis, *pThis);
+ }
+ std::fill(vResult.begin(), vResult.begin() + iLeft, oOpts.outOfRangeVal);
+ std::fill(vResult.begin() + iRight, vResult.begin() + oCurExtent.xStop,
+ oOpts.outOfRangeVal);
+}
+
+/// Process the part of the first line to the left of the observer.
+///
+/// @param iStart X coordinate of the first cell to the left of the observer to be procssed.
+/// @param iEnd X coordinate one past the last cell to be processed.
+/// @param vResult Vector in which to store the visibility/height results.
+/// @param vThisLineVal Height of each cell in the line being processed.
+void ViewshedExecutor::processFirstLineLeft(int iStart, int iEnd,
+ std::vector &vResult,
+ std::vector &vThisLineVal)
+{
+ // If end is to the right of start, everything is taken care of by right processing.
+ if (iEnd >= iStart)
+ return;
+
+ iStart = oCurExtent.clampX(iStart);
+
+ double *pThis = vThisLineVal.data() + iStart;
+
+ // If the start cell is next to the observer, just mark it visible.
+ if (iStart + 1 == m_nX || iStart + 1 == oCurExtent.xStop)
+ {
+ if (oOpts.outputMode == OutputMode::Normal)
+ vResult[iStart] = oOpts.visibleVal;
+ else
+ setOutput(vResult[iStart], *pThis, *pThis);
+ iStart--;
+ pThis--;
+ }
+
+ // Go from the observer to the left, calculating Z as we go.
+ for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--)
+ {
+ int nXOffset = std::abs(iPixel - m_nX);
+ double dfZ = CalcHeightLine(nXOffset, *(pThis + 1));
+ setOutput(vResult[iPixel], *pThis, dfZ);
+ }
+ // For cells outside of the [start, end) range, set the outOfRange value.
+ std::fill(vResult.begin(), vResult.begin() + iEnd + 1, oOpts.outOfRangeVal);
+}
+
+/// Process the part of the first line to the right of the observer.
+///
+/// @param iStart X coordinate of the first cell to the right of the observer to be processed.
+/// @param iEnd X coordinate one past the last cell to be processed.
+/// @param vResult Vector in which to store the visibility/height results.
+/// @param vThisLineVal Height of each cell in the line being processed.
+void ViewshedExecutor::processFirstLineRight(int iStart, int iEnd,
+ std::vector &vResult,
+ std::vector &vThisLineVal)
+{
+ // If start is to the right of end, everything is taken care of by left processing.
+ if (iStart >= iEnd)
+ return;
+
+ iStart = oCurExtent.clampX(iStart);
+
+ double *pThis = vThisLineVal.data() + iStart;
+
+ // If the start cell is next to the observer, just mark it visible.
+ if (iStart - 1 == m_nX || iStart == oCurExtent.xStart)
+ {
+ if (oOpts.outputMode == OutputMode::Normal)
+ vResult[iStart] = oOpts.visibleVal;
+ else
+ setOutput(vResult[iStart], *pThis, *pThis);
+ iStart++;
+ pThis++;
+ }
+
+ // Go from the observer to the right, calculating Z as we go.
+ for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++)
+ {
+ int nXOffset = std::abs(iPixel - m_nX);
+ double dfZ = CalcHeightLine(nXOffset, *(pThis - 1));
+ setOutput(vResult[iPixel], *pThis, dfZ);
+ }
+ // For cells outside of the [start, end) range, set the outOfRange value.
+ std::fill(vResult.begin() + iEnd, vResult.end(), oOpts.outOfRangeVal);
+}
+
+/// Process a line to the left of the observer.
+///
+/// @param nYOffset Offset of the line being processed from the observer
+/// @param iStart X coordinate of the first cell to the left of the observer to be processed.
+/// @param iEnd X coordinate one past the last cell to be processed.
+/// @param vResult Vector in which to store the visibility/height results.
+/// @param vThisLineVal Height of each cell in the line being processed.
+/// @param vLastLineVal Observable height of each cell in the previous line processed.
+void ViewshedExecutor::processLineLeft(int nYOffset, int iStart, int iEnd,
+ std::vector &vResult,
+ std::vector &vThisLineVal,
+ std::vector &vLastLineVal)
+{
+ // If start to the left of end, everything is taken care of by processing right.
+ if (iStart <= iEnd)
+ return;
+ iStart = oCurExtent.clampX(iStart);
+
+ nYOffset = std::abs(nYOffset);
+ double *pThis = vThisLineVal.data() + iStart;
+ double *pLast = vLastLineVal.data() + iStart;
+
+ // If the observer is to the right of the raster, mark the first cell to the left as
+ // visible. This may mark an out-of-range cell with a value, but this will be fixed
+ // with the out of range assignment at the end.
+ if (iStart == oCurExtent.xStop - 1)
+ {
+ if (oOpts.outputMode == OutputMode::Normal)
+ vResult[iStart] = oOpts.visibleVal;
+ else
+ setOutput(vResult[iStart], *pThis, *pThis);
+ iStart--;
+ pThis--;
+ pLast--;
+ }
+
+ // Go from the observer to the left, calculating Z as we go.
+ for (int iPixel = iStart; iPixel > iEnd; iPixel--, pThis--, pLast--)
+ {
+ int nXOffset = std::abs(iPixel - m_nX);
+ double dfZ;
+ if (nXOffset == nYOffset)
+ {
+ if (nXOffset == 1)
+ dfZ = *pThis;
+ else
+ dfZ = CalcHeightLine(nXOffset, *(pLast + 1));
+ }
+ else
+ dfZ =
+ oZcalc(nXOffset, nYOffset, *(pThis + 1), *pLast, *(pLast + 1));
+ setOutput(vResult[iPixel], *pThis, dfZ);
+ }
+
+ // For cells outside of the [start, end) range, set the outOfRange value.
+ std::fill(vResult.begin(), vResult.begin() + iEnd + 1, oOpts.outOfRangeVal);
+}
+
+/// Process a line to the right of the observer.
+///
+/// @param nYOffset Offset of the line being processed from the observer
+/// @param iStart X coordinate of the first cell to the right of the observer to be processed.
+/// @param iEnd X coordinate one past the last cell to be processed.
+/// @param vResult Vector in which to store the visibility/height results.
+/// @param vThisLineVal Height of each cell in the line being processed.
+/// @param vLastLineVal Observable height of each cell in the previous line processed.
+void ViewshedExecutor::processLineRight(int nYOffset, int iStart, int iEnd,
+ std::vector &vResult,
+ std::vector &vThisLineVal,
+ std::vector &vLastLineVal)
+{
+ // If start is to the right of end, everything is taken care of by processing left.
+ if (iStart >= iEnd)
+ return;
+ iStart = oCurExtent.clampX(iStart);
+
+ nYOffset = std::abs(nYOffset);
+ double *pThis = vThisLineVal.data() + iStart;
+ double *pLast = vLastLineVal.data() + iStart;
+
+ // If the observer is to the left of the raster, mark the first cell to the right as
+ // visible. This may mark an out-of-range cell with a value, but this will be fixed
+ // with the out of range assignment at the end.
+ if (iStart == 0)
+ {
+ if (oOpts.outputMode == OutputMode::Normal)
+ vResult[iStart] = oOpts.visibleVal;
+ else
+ setOutput(vResult[0], *pThis, *pThis);
+ iStart++;
+ pThis++;
+ pLast++;
+ }
+
+ // Go from the observer to the right, calculating Z as we go.
+ for (int iPixel = iStart; iPixel < iEnd; iPixel++, pThis++, pLast++)
+ {
+ int nXOffset = std::abs(iPixel - m_nX);
+ double dfZ;
+ if (nXOffset == nYOffset)
+ {
+ if (nXOffset == 1)
+ dfZ = *pThis;
+ else
+ dfZ = CalcHeightLine(nXOffset, *(pLast - 1));
+ }
+ else
+ dfZ =
+ oZcalc(nXOffset, nYOffset, *(pThis - 1), *pLast, *(pLast - 1));
+ setOutput(vResult[iPixel], *pThis, dfZ);
+ }
+
+ // For cells outside of the [start, end) range, set the outOfRange value.
+ std::fill(vResult.begin() + iEnd, vResult.end(), oOpts.outOfRangeVal);
+}
+
+/// Process a line above or below the observer.
+///
+/// @param nLine Line number being processed.
+/// @param vLastLineVal Vector in which to store the read line. Becomes the last line
+/// in further processing.
+/// @return True on success, false otherwise.
+bool ViewshedExecutor::processLine(int nLine, std::vector &vLastLineVal)
+{
+ int nYOffset = nLine - m_nY;
+ std::vector vResult(oOutExtent.xSize());
+ std::vector vThisLineVal(oOutExtent.xSize());
+
+ if (!readLine(nLine, vThisLineVal.data()))
+ return false;
+
+ // In DEM mode the base is the input DEM value.
+ // In ground mode the base is zero.
+ if (oOpts.outputMode == OutputMode::DEM)
+ vResult = vThisLineVal;
+
+ // Adjust height of the read line.
+ const auto [iLeft, iRight] = adjustHeight(nYOffset, vThisLineVal);
+
+ // Handle the initial position on the line.
+ if (oCurExtent.containsX(m_nX))
+ {
+ if (iLeft < iRight)
+ {
+ double dfZ;
+ if (std::abs(nYOffset) == 1)
+ dfZ = vThisLineVal[m_nX];
+ else
+ dfZ = CalcHeightLine(nYOffset, vLastLineVal[m_nX]);
+ setOutput(vResult[m_nX], vThisLineVal[m_nX], dfZ);
+ }
+ else
+ vResult[m_nX] = oOpts.outOfRangeVal;
+ }
+
+ // process left half then right half of line
+ CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
+ pQueue->SubmitJob(
+ [&, left = iLeft]()
+ {
+ processLineLeft(nYOffset, m_nX - 1, left - 1, vResult, vThisLineVal,
+ vLastLineVal);
+ });
+ pQueue->SubmitJob(
+ [&, right = iRight]()
+ {
+ processLineRight(nYOffset, m_nX + 1, right, vResult, vThisLineVal,
+ vLastLineVal);
+ });
+ pQueue->WaitCompletion();
+
+ // Make the current line the last line.
+ vLastLineVal = std::move(vThisLineVal);
+
+ if (!writeLine(nLine, vResult))
+ return false;
+
+ return oProgress.lineComplete();
+}
+
+/// Run the viewshed computation
+/// @return Success as true or false.
+bool ViewshedExecutor::run()
+{
+ std::vector vFirstLineVal(oCurExtent.xSize());
+ if (!processFirstLine(vFirstLineVal))
+ return false;
+
+ if (oOpts.cellMode == CellMode::Edge)
+ oZcalc = doEdge;
+ else if (oOpts.cellMode == CellMode::Diagonal)
+ oZcalc = doDiagonal;
+ else if (oOpts.cellMode == CellMode::Min)
+ oZcalc = doMin;
+ else if (oOpts.cellMode == CellMode::Max)
+ oZcalc = doMax;
+
+ // scan upwards
+ int yStart = oCurExtent.clampY(m_nY);
+ std::atomic err(false);
+ CPLJobQueuePtr pQueue = m_pool.CreateJobQueue();
+ pQueue->SubmitJob(
+ [&]()
+ {
+ std::vector vLastLineVal = vFirstLineVal;
+
+ for (int nLine = yStart - 1; nLine >= oCurExtent.yStart && !err;
+ nLine--)
+ if (!processLine(nLine, vLastLineVal))
+ err = true;
+ });
+
+ // scan downwards
+ pQueue->SubmitJob(
+ [&]()
+ {
+ std::vector vLastLineVal = vFirstLineVal;
+
+ for (int nLine = yStart + 1; nLine < oCurExtent.yStop && !err;
+ nLine++)
+ if (!processLine(nLine, vLastLineVal))
+ err = true;
+ });
+ return true;
+}
+
+} // namespace viewshed
+} // namespace gdal
diff --git a/alg/viewshed/viewshed_executor.h b/alg/viewshed/viewshed_executor.h
new file mode 100644
index 000000000000..291581e79b07
--- /dev/null
+++ b/alg/viewshed/viewshed_executor.h
@@ -0,0 +1,103 @@
+/******************************************************************************
+ *
+ * Project: Viewshed Generation
+ * Purpose: Core algorithm implementation for viewshed generation.
+ * Author: Tamas Szekeres, szekerest@gmail.com
+ *
+ * (c) 2024 info@hobu.co
+ *
+ ******************************************************************************
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#pragma once
+
+#include
+#include
+
+#include "gdal_priv.h"
+#include "cpl_worker_thread_pool.h"
+
+#include "viewshed_types.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+class Progress;
+
+/// Executes a viewshed computation on a source band, placing the result
+/// in the destination band.
+class ViewshedExecutor
+{
+ public:
+ ViewshedExecutor(GDALRasterBand &srcBand, GDALRasterBand &dstBand, int nX,
+ int nY, const Window &oOutExtent, const Window &oCurExtent,
+ const Options &opts, Progress &oProgress);
+ bool run();
+
+ private:
+ CPLWorkerThreadPool m_pool;
+ GDALRasterBand &m_srcBand;
+ GDALRasterBand &m_dstBand;
+ const Window oOutExtent;
+ const Window oCurExtent;
+ const int m_nX;
+ const int m_nY;
+ const Options oOpts;
+ Progress &oProgress;
+ double m_dfHeightAdjFactor{0};
+ double m_dfMaxDistance2;
+ double m_dfZObserver{0};
+ std::mutex iMutex{};
+ std::mutex oMutex{};
+ std::array m_adfTransform{0, 1, 0, 0, 0, 1};
+ double (*oZcalc)(int, int, double, double, double){};
+
+ double calcHeightAdjFactor();
+ void setOutput(double &dfResult, double &dfCellVal, double dfZ);
+ bool readLine(int nLine, double *data);
+ bool writeLine(int nLine, std::vector &vResult);
+ bool processLine(int nLine, std::vector &vLastLineVal);
+ bool processFirstLine(std::vector &vLastLineVal);
+ void processFirstLineLeft(int iStart, int iEnd,
+ std::vector &vResult,
+ std::vector &vThisLineVal);
+ void processFirstLineRight(int iStart, int iEnd,
+ std::vector &vResult,
+ std::vector &vThisLineVal);
+ void processFirstLineTopOrBottom(int iLeft, int iRight,
+ std::vector &vResult,
+ std::vector &vThisLineVal);
+ void processLineLeft(int nYOffset, int iStart, int iEnd,
+ std::vector &vResult,
+ std::vector &vThisLineVal,
+ std::vector &vLastLineVal);
+ void processLineRight(int nYOffset, int iStart, int iEnd,
+ std::vector &vResult,
+ std::vector &vThisLineVal,
+ std::vector &vLastLineVal);
+ std::pair adjustHeight(int iLine,
+ std::vector &thisLineVal);
+};
+
+} // namespace viewshed
+} // namespace gdal
diff --git a/alg/viewshed/viewshed_types.h b/alg/viewshed/viewshed_types.h
new file mode 100644
index 000000000000..fa1f0b2b9ead
--- /dev/null
+++ b/alg/viewshed/viewshed_types.h
@@ -0,0 +1,176 @@
+/****************************************************************************
+ * (c) 2024 info@hobu.co
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+#ifndef VIEWSHED_TYPES_H_INCLUDED
+#define VIEWSHED_TYPES_H_INCLUDED
+
+#include
+#include
+
+#include "gdal_priv.h"
+
+namespace gdal
+{
+namespace viewshed
+{
+
+/// Unique pointer to GDAL dataset.
+using DatasetPtr = std::unique_ptr;
+
+/**
+ * Raster output mode.
+ */
+enum class OutputMode
+{
+ Normal, //!< Normal output mode (visibility only)
+ DEM, //!< Output height from DEM
+ Ground, //!< Output height from ground
+ Cumulative //!< Output observability heat map
+};
+
+/**
+ * Cell height calculation mode.
+ */
+enum class CellMode
+{
+ Diagonal, //!< Diagonal Mode
+ Edge, //!< Edge Mode
+ Max, //!< Maximum value produced by Diagonal and Edge mode
+ Min //!< Minimum value produced by Diagonal and Edge mode
+};
+
+/**
+ * A point.
+ */
+struct Point
+{
+ double x; //!< X value
+ double y; //!< Y value
+ double z; //!< Z value
+};
+
+/**
+ * Options for viewshed generation.
+ */
+struct Options
+{
+ Point observer{0, 0, 0}; //!< x, y, and z of the observer
+ double visibleVal{255}; //!< raster output value for visible pixels.
+ double invisibleVal{0}; //!< raster output value for non-visible pixels.
+ double outOfRangeVal{
+ 0}; //!< raster output value for pixels outside of max distance.
+ double nodataVal{-1}; //!< raster output value for pixels with no data
+ double targetHeight{0.0}; //!< target height above the DEM surface
+ double maxDistance{
+ 0.0}; //!< maximum distance from observer to compute value
+ double curveCoeff{.85714}; //!< coefficient for atmospheric refraction
+ OutputMode outputMode{OutputMode::Normal}; //!< Output information.
+ //!< Normal, Height from DEM or Height from ground
+ std::string outputFormat{}; //!< output raster format
+ std::string outputFilename{}; //!< output raster filename
+ CPLStringList creationOpts{}; //!< options for output raster creation
+ CellMode cellMode{CellMode::Edge}; //!< Mode of cell height calculation.
+ int observerSpacing{10}; //!< Observer spacing in cumulative mode.
+ uint8_t numJobs{3}; //!< Relative number of jobs in cumulative mode.
+};
+
+/**
+ * A window in a raster including pixels in [xStart, xStop) and [yStart, yStop).
+ */
+struct Window
+{
+ int xStart{}; //!< X start position
+ int xStop{}; //!< X end position
+ int yStart{}; //!< Y start position
+ int yStop{}; //!< Y end position
+
+ /// \brief Window size in the X direction.
+ int xSize() const
+ {
+ return xStop - xStart;
+ }
+
+ /// \brief Window size in the Y direction.
+ int ySize() const
+ {
+ return yStop - yStart;
+ }
+
+ /// \brief Number of cells.
+ size_t size() const
+ {
+ return static_cast(xSize()) * ySize();
+ }
+
+ /// \brief Determine if the X window contains the index.
+ /// \param nX Index to check
+ /// \return True if the index is contained, false otherwise.
+ bool containsX(int nX) const
+ {
+ return nX >= xStart && nX < xStop;
+ }
+
+ /// \brief Determine if the Y window contains the index.
+ /// \param nY Index to check
+ /// \return True if the index is contained, false otherwise.
+ bool containsY(int nY) const
+ {
+ return nY >= xStart && nY < yStop;
+ }
+
+ /// \brief Determine if the window contains the index.
+ /// \param nX X coordinate of the index to check
+ /// \param nY Y coordinate of the index to check
+ /// \return True if the index is contained, false otherwise.
+ bool contains(int nX, int nY) const
+ {
+ return containsX(nX) && containsY(nY);
+ }
+
+ /// \brief Clamp the argument to be in the window in the X dimension.
+ /// \param nX Value to clamp.
+ /// \return Clamped value.
+ int clampX(int nX) const
+ {
+ return xSize() ? std::clamp(nX, xStart, xStop - 1) : xStart;
+ }
+
+ /// \brief Clamp the argument to be in the window in the Y dimension.
+ /// \param nY Value to clamp.
+ /// \return Clamped value.
+ int clampY(int nY) const
+ {
+ return ySize() ? std::clamp(nY, yStart, yStop - 1) : yStart;
+ }
+
+ /// \brief Shift the X dimension by nShift.
+ /// \param nShift Amount to shift
+ void shiftX(int nShift)
+ {
+ xStart += nShift;
+ xStop += nShift;
+ }
+};
+
+} // namespace viewshed
+} // namespace gdal
+
+#endif
diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt
index 554ca13f1955..4ce116ae1b01 100644
--- a/apps/CMakeLists.txt
+++ b/apps/CMakeLists.txt
@@ -158,6 +158,9 @@ if (BUILD_APPS)
add_executable(gdalasyncread EXCLUDE_FROM_ALL gdalasyncread.cpp)
add_executable(gdalwarpsimple EXCLUDE_FROM_ALL gdalwarpsimple.c)
add_executable(multireadtest EXCLUDE_FROM_ALL multireadtest.cpp)
+ if(NOT MSVC AND CMAKE_THREAD_LIBS_INIT)
+ target_link_libraries(multireadtest PRIVATE ${CMAKE_THREAD_LIBS_INIT})
+ endif()
add_executable(test_ogrsf test_ogrsf.cpp)
add_executable(testreprojmulti EXCLUDE_FROM_ALL testreprojmulti.cpp)
diff --git a/apps/gdal_contour.cpp b/apps/gdal_contour.cpp
index 7dd5f9f9bddb..1248f7b2c190 100644
--- a/apps/gdal_contour.cpp
+++ b/apps/gdal_contour.cpp
@@ -155,13 +155,9 @@ GDALContourAppOptionsGetParser(GDALContourOptions *psOptions)
.help(_("Generate levels on an exponential scale: base ^ k, for k an "
"integer."));
- argParser->add_argument("-fl")
- .metavar("")
- .nargs(argparse::nargs_pattern::at_least_one)
- .scan<'g', double>()
- .action([psOptions](const std::string &s)
- { psOptions->adfFixedLevels.push_back(CPLAtof(s.c_str())); })
- .help(_("Name one or more \"fixed levels\" to extract."));
+ // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
+ argParser->add_argument("-fl").scan<'g', double>().metavar("").help(
+ _("Name one or more \"fixed levels\" to extract."));
argParser->add_argument("-off")
.metavar("")
@@ -246,6 +242,7 @@ MAIN_START(argc, argv)
if (argc < 2)
{
+
try
{
GDALContourOptions sOptions;
@@ -261,11 +258,54 @@ MAIN_START(argc, argv)
}
GDALContourOptions sOptions;
+ CPLStringList aosArgv;
try
{
+ /* -------------------------------------------------------------------- */
+ /* Pre-processing for custom syntax that ArgumentParser does not */
+ /* support. */
+ /* -------------------------------------------------------------------- */
+ for (int i = 1; i < argc && argv != nullptr && argv[i] != nullptr; i++)
+ {
+ // argparser is confused by arguments that have at_least_one
+ // cardinality, if they immediately precede positional arguments.
+ if (EQUAL(argv[i], "-fl") && argv[i + 1])
+ {
+ if (strchr(argv[i + 1], ' '))
+ {
+ const CPLStringList aosTokens(
+ CSLTokenizeString(argv[i + 1]));
+ for (const char *pszToken : aosTokens)
+ {
+ sOptions.adfFixedLevels.push_back(CPLAtof(pszToken));
+ }
+ i += 1;
+ }
+ else
+ {
+ auto isNumeric = [](const char *pszArg) -> bool
+ {
+ char *pszEnd = nullptr;
+ CPLStrtod(pszArg, &pszEnd);
+ return pszEnd != nullptr && pszEnd[0] == '\0';
+ };
+
+ while (i < argc - 1 && isNumeric(argv[i + 1]))
+ {
+ sOptions.adfFixedLevels.push_back(CPLAtof(argv[i + 1]));
+ i += 1;
+ }
+ }
+ }
+ else
+ {
+ aosArgv.AddString(argv[i]);
+ }
+ }
+
auto argParser = GDALContourAppOptionsGetParser(&sOptions);
- argParser->parse_args_without_binary_name(argv + 1);
+ argParser->parse_args_without_binary_name(aosArgv.List());
if (sOptions.dfInterval == 0.0 && sOptions.adfFixedLevels.empty() &&
sOptions.dfExpBase == 0.0)
diff --git a/apps/gdal_rasterize_bin.cpp b/apps/gdal_rasterize_bin.cpp
index 772a5cc9d766..5121cfb3e27d 100644
--- a/apps/gdal_rasterize_bin.cpp
+++ b/apps/gdal_rasterize_bin.cpp
@@ -36,33 +36,11 @@
/* Usage() */
/************************************************************************/
-static void Usage(bool bIsError, const char *pszErrorMsg = nullptr)
+static void Usage()
{
- fprintf(
- bIsError ? stderr : stdout,
- "Usage: gdal_rasterize [--help] [--help-general]\n"
- " [-b ]... [-i] [-at]\n"
- " [-oo =]...\n"
- " {[-burn ]... | [-a ] | [-3d]} [-add]\n"
- " [-l ]... [-where ] "
- "[-sql |@]\n"
- " [-dialect ] [-of ] [-a_srs ] [-to "
- "=]...\n"
- " [-co =]... [-a_nodata ] [-init "
- "]...\n"
- " [-te ] [-tr ] [-tap] "
- "[-ts "
- "]\n"
- " [-ot "
- "{Byte/Int8/Int16/UInt16/UInt32/Int32/UInt64/Int64/Float32/Float64/\n"
- " CInt16/CInt32/CFloat32/CFloat64}] [-optim "
- "{AUTO|VECTOR|RASTER}] [-q]\n"
- " \n");
-
- if (pszErrorMsg != nullptr)
- fprintf(stderr, "\nFAILURE: %s\n", pszErrorMsg);
- exit(bIsError ? 1 : 0);
+ fprintf(stderr, "%s\n", GDALRasterizeAppGetParserUsage().c_str());
+ exit(1);
}
/************************************************************************/
@@ -71,6 +49,7 @@ static void Usage(bool bIsError, const char *pszErrorMsg = nullptr)
MAIN_START(argc, argv)
{
+
/* Check strict compilation and runtime library version as we use C++ API */
if (!GDAL_CHECK_VERSION(argv[0]))
exit(1);
@@ -78,51 +57,36 @@ MAIN_START(argc, argv)
EarlySetConfigOptions(argc, argv);
/* -------------------------------------------------------------------- */
- /* Generic arg processing. */
+ /* Register standard GDAL drivers, and process generic GDAL */
+ /* command options. */
/* -------------------------------------------------------------------- */
GDALAllRegister();
argc = GDALGeneralCmdLineProcessor(argc, &argv, 0);
if (argc < 1)
exit(-argc);
- for (int i = 0; i < argc; i++)
- {
- if (EQUAL(argv[i], "--utility_version"))
- {
- printf("%s was compiled against GDAL %s and "
- "is running against GDAL %s\n",
- argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
- CSLDestroy(argv);
- return 0;
- }
- else if (EQUAL(argv[i], "--help"))
- {
- Usage(false);
- }
- }
+ /* -------------------------------------------------------------------- */
+ /* Generic arg processing. */
+ /* -------------------------------------------------------------------- */
GDALRasterizeOptionsForBinary sOptionsForBinary;
- // coverity[tainted_data]
- GDALRasterizeOptions *psOptions =
- GDALRasterizeOptionsNew(argv + 1, &sOptionsForBinary);
+ std::unique_ptr
+ psOptions{GDALRasterizeOptionsNew(argv + 1, &sOptionsForBinary),
+ GDALRasterizeOptionsFree};
+
CSLDestroy(argv);
- if (psOptions == nullptr)
+ if (!psOptions)
{
- Usage(true);
+ Usage();
}
if (!(sOptionsForBinary.bQuiet))
{
- GDALRasterizeOptionsSetProgress(psOptions, GDALTermProgress, nullptr);
+ GDALRasterizeOptionsSetProgress(psOptions.get(), GDALTermProgress,
+ nullptr);
}
- if (sOptionsForBinary.osSource.empty())
- Usage(true, "No input file specified.");
-
- if (!sOptionsForBinary.bDestSpecified)
- Usage(true, "No output file specified.");
-
/* -------------------------------------------------------------------- */
/* Open input file. */
/* -------------------------------------------------------------------- */
@@ -186,16 +150,19 @@ MAIN_START(argc, argv)
}
int bUsageError = FALSE;
- GDALDatasetH hRetDS = GDALRasterize(sOptionsForBinary.osDest.c_str(),
- hDstDS, hInDS, psOptions, &bUsageError);
+ GDALDatasetH hRetDS =
+ GDALRasterize(sOptionsForBinary.osDest.c_str(), hDstDS, hInDS,
+ psOptions.get(), &bUsageError);
+
if (bUsageError == TRUE)
- Usage(true);
+ Usage();
+
int nRetCode = hRetDS ? 0 : 1;
GDALClose(hInDS);
+
if (GDALClose(hRetDS) != CE_None)
nRetCode = 1;
- GDALRasterizeOptionsFree(psOptions);
GDALDestroyDriverManager();
diff --git a/apps/gdal_rasterize_lib.cpp b/apps/gdal_rasterize_lib.cpp
index 493922ee1068..e3b7d346bb0a 100644
--- a/apps/gdal_rasterize_lib.cpp
+++ b/apps/gdal_rasterize_lib.cpp
@@ -50,17 +50,314 @@
#include "ogr_api.h"
#include "ogr_core.h"
#include "ogr_srs_api.h"
+#include "gdalargumentparser.h"
/************************************************************************/
-/* ArgIsNumericRasterize() */
+/* GDALRasterizeOptions() */
/************************************************************************/
-static bool ArgIsNumericRasterize(const char *pszArg)
+struct GDALRasterizeOptions
+{
+ std::vector anBandList{};
+ std::vector adfBurnValues{};
+ bool bInverse = false;
+ std::string osFormat{};
+ bool b3D = false;
+ GDALProgressFunc pfnProgress = GDALDummyProgress;
+ void *pProgressData = nullptr;
+ std::vector aosLayers{};
+ std::string osSQL{};
+ std::string osDialect{};
+ std::string osBurnAttribute{};
+ std::string osWHERE{};
+ CPLStringList aosRasterizeOptions{};
+ CPLStringList aosTO{};
+ double dfXRes = 0;
+ double dfYRes = 0;
+ CPLStringList aosCreationOptions{};
+ GDALDataType eOutputType = GDT_Unknown;
+ std::vector adfInitVals{};
+ std::string osNoData{};
+ OGREnvelope sEnvelop{};
+ int nXSize = 0;
+ int nYSize = 0;
+ OGRSpatialReference oOutputSRS{};
+
+ bool bTargetAlignedPixels = false;
+ bool bCreateOutput = false;
+};
+
+/************************************************************************/
+/* GDALRasterizeOptionsGetParser() */
+/************************************************************************/
+static std::unique_ptr
+GDALRasterizeOptionsGetParser(GDALRasterizeOptions *psOptions,
+ GDALRasterizeOptionsForBinary *psOptionsForBinary)
{
- char *pszEnd = nullptr;
- CPLStrtod(pszArg, &pszEnd);
- return pszEnd != nullptr && pszEnd[0] == '\0';
+ auto argParser = std::make_unique(
+ "gdal_rasterize", /* bForBinary=*/psOptionsForBinary != nullptr);
+
+ argParser->add_description(_("Burns vector geometries into a raster."));
+
+ argParser->add_epilog(
+ _("This program burns vector geometries (points, lines, and polygons) "
+ "into the raster band(s) of a raster image."));
+
+ // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
+ argParser->add_argument("-b")
+ .metavar("")
+ .append()
+ .scan<'i', int>()
+ //.nargs(argparse::nargs_pattern::at_least_one)
+ .help(_("The band(s) to burn values into."));
+
+ argParser->add_argument("-i")
+ .flag()
+ .store_into(psOptions->bInverse)
+ .help(_("Invert rasterization."));
+
+ argParser->add_argument("-at")
+ .flag()
+ .action(
+ [psOptions](const std::string &) {
+ psOptions->aosRasterizeOptions.SetNameValue("ALL_TOUCHED",
+ "TRUE");
+ })
+ .help(_("Enables the ALL_TOUCHED rasterization option."));
+
+ // Mutually exclusive options: -burn, -3d, -a
+ {
+ // Required if options for binary
+ auto &group = argParser->add_mutually_exclusive_group(
+ psOptionsForBinary != nullptr);
+
+ // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
+ group.add_argument("-burn")
+ .metavar("")
+ .scan<'g', double>()
+ .append()
+ //.nargs(argparse::nargs_pattern::at_least_one)
+ .help(_("A fixed value to burn into the raster band(s)."));
+
+ group.add_argument("-a")
+ .metavar("")
+ .store_into(psOptions->osBurnAttribute)
+ .help(_("Name of the field in the input layer to get the burn "
+ "values from."));
+
+ group.add_argument("-3d")
+ .flag()
+ .store_into(psOptions->b3D)
+ .action(
+ [psOptions](const std::string &) {
+ psOptions->aosRasterizeOptions.SetNameValue(
+ "BURN_VALUE_FROM", "Z");
+ })
+ .help(_("Indicates that a burn value should be extracted from the "
+ "\"Z\" values of the feature."));
+ }
+
+ argParser->add_argument("-add")
+ .flag()
+ .action(
+ [psOptions](const std::string &) {
+ psOptions->aosRasterizeOptions.SetNameValue("MERGE_ALG", "ADD");
+ })
+ .help(_("Instead of burning a new value, this adds the new value to "
+ "the existing raster."));
+
+ // Undocumented
+ argParser->add_argument("-chunkysize")
+ .flag()
+ .hidden()
+ .action(
+ [psOptions](const std::string &s) {
+ psOptions->aosRasterizeOptions.SetNameValue("CHUNKYSIZE",
+ s.c_str());
+ });
+
+ // Mutually exclusive -l, -sql
+ {
+ auto &group = argParser->add_mutually_exclusive_group(false);
+
+ group.add_argument("-l")
+ .metavar("")
+ .append()
+ .store_into(psOptions->aosLayers)
+ .help(_("Name of the layer(s) to process."));
+
+ group.add_argument("-sql")
+ .metavar("