OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.78k stars 2.5k forks source link

gdal_viewshed inaccurate/wrong results #10588

Open leonidb555 opened 1 month ago

leonidb555 commented 1 month ago

What is the bug?

Problem description

The gdal_viewshed tool sometimes produces inaccurate results. This issue becomes evident when comparing its output with other viewshed engines or when performing a line-of-sight (LOS) analysis over specific areas using a DEM profile.

Here is an example of gdal_viewshed calculation compared with other engines and with profile LOS.

Example

DEM Source : SRTM 3sec Observer Coordinates: 588900, 3504400 Observer Height: 20m above ground Target Height: 1m above ground Distance: 100000m Curvature Coefficient: 0.85714 LOS Target Coordinates: 551014.9, 3596874.7 (for LOS check only) Project CRS: _+proj=tmerc +lat_0=0 +lon_0=34.2 +k=0.9996 +x_0=500000 +y_0=0 +datum=WGS84 +units=m +nodefs +type=crs

Comparison Results

  1. GDAL vs other engines Viewshed Calculation:

    image

  2. Line-of-Sight (LOS) Check:

    image

Observations

In the LOS check, it is clear that there is an obstruction at approximately 11 km from the observer. The target point should not be visible from the observer's location due to this obstruction, yet gdal_viewshed indicates otherwise.

Expected Behavior

The gdal_viewshed tool should accurately reflect the obstructions in the line of sight and produce results consistent with other viewshed engines and LOS analysis.

Steps to reproduce the issue

gdal_viewshed -ox 588900.0 -oy 3504400.0 -oz 20.0 -tz 1.0 -md 100000.0 -f GTiff -cc 0.85714 srtm_3.tiff out.tiff

Versions and provenance

Computed with GDAL Plugin on QGIS on Windows 11 machine

QGIS version: 3.36.3-Maidenhead QGIS code revision: 2df96554 Qt version: 5.15.13 Python version: 3.12.3 GDAL version: 3.9.0 GEOS version: 3.12.1-CAPI-1.18.1 PROJ version: Rel. 9.4.0, March 1st, 2024 PDAL version: 2.6.3 (git-version: b5523a)

Additional context

Please let me know if further details or data are required to address this issue. Thanks!

rouault commented 1 month ago

CC @abellgithub

abellgithub commented 1 month ago

The GDAL implementation does not use sightlines. It uses the the Wang method described in the documentation. I am not sure how you've generated the line-of-sight analysis so I can't judge whether it is relevant to the GDAL viewshed algorithm or not.

There are two cell modes that can be used (diagonal and normal) but these can only be accessed through the API, not the command line. I have no idea whether the other mode might make a difference, but you could try.

leonidb555 commented 1 month ago

Thank you for your response.

The line-of-sight in my example is a straight line above height profile calculated directly from DEM pixels without interpolation. I mentioned it to provide additional evidence of potential issues with the reliability of the gdal_viewshed algorithm, beyond comparisons with other engines.

I have tried running the viewshed with different mode settings (GVM_Diagonal, GVM_Edge, GVM_Max, GVM_Min), but these settings don’t seem to have a significant impact on the results.

Regarding the Wang algorithm, you are correct that it does not use line-of-sight calculations. However, I would expect that the results should still be consistent with a simple line-of-sight check.

abellgithub commented 1 month ago

You might use DEM mode to see the actual heights used during the calculation.

The algorithm does not assume that the surface is an array of rectangular solids. If your own line-of-sight calculation does this, the results will be different. There is no single correct answer to the viewshed problem.

leonidb555 commented 4 weeks ago

Hi @abellgithub,

I've checked the DEM mode, and it produces consistent results with the NORMAL mode, but the same inaccuracy persists.

I understand that there isn't a single "correct" answer for viewshed analysis since DEMs are approximations of the real world and different algorithms interpret DEM heights in various ways. However, if there is a visible obstruction that prevents direct line of sight between two points, this obstruction should be accurately reflected in the viewshed results. In other words, if a hill in the DEM blocks the view, the area behind that hill should be marked as non-visible.

Additionally, I've conducted further tests and found that densifying the DEM grid (by simple oversampling without altering pixel values—e.g., dividing each pixel into 4, 9, 16, etc.) results in viewshed outputs that are more reliable and consistent with other engines. This leads me to suspect that the issue might be related to the "skipping" of DEM pixels in the original calculation.

Given these findings, it appears there may be an underlying issue affecting the accuracy of the gdal_viewshed tool. I would appreciate it if this could be reviewed to ensure that the tool provides results that accurately reflect the visible areas considering DEM obstructions.

Thank you for your attention to this matter.

jratike80 commented 4 weeks ago

Do you know if the Wang algorithm is supposed to give same results than line-of-sight? This article https://www.sciencedirect.com/science/article/pii/S0098300422001625 says about the reference plane approach

the reference plane (RP) algorithm (Wang et al., 1996) adopt the approximate calculation method to reduce the computation complexity to O(r2) by sacrificing accuracy

Can you find as a reference other utilities that implement the RP algorithm but give different results than gdal_viewshed?

rouault commented 4 weeks ago

Note that we have now GDALIsLineOfSightVisible that could potentially be used to build a viewshed alternative based on it (although it would likely be much slower, probably O(r^3))

abellgithub commented 4 weeks ago

The Wang algorithm does interpolation between points that might block visibility. Choice of those points is determined by the cell mode (mentioned earlier). In your test, when you make cells smaller without more data you are reducing the ability of the algorithm to interpolate between known data.

If you draw a line between observer and target and assume that any cell whose height is greater than the line at that cell location will block observability (which is what I believe the mentioned Breshenham algorithm to be doing), you will get different (not necessarily worse or even less accurate) results than computed with the Wang algorithm. This is because the Wang algorithm uses an interpolated value between blocking cells to determine the visibility of a successive cell.

There is only so much data available in a raster. If a line does not exactly intersect a raster location, one must decide how to determine an obstruction since the terrain at that location is unknown. Even if a line does pass directly through a raster point, it's likely that the value of the raster is itself the result of some data interpolation, leading to inaccuracies.

The fact that the results you see are different than those from some other algorithm isn't sufficient to indicate a bug. Wang isn't a complicated algorithm and if you can make a simple test case that exhibits an error I'd be happy to spend some time looking at it, but that error needs to be within the context of the Wang algorithm. If you're not happy with the Wang algorithm, I would suggest you continue to use other software.

abellgithub commented 4 weeks ago

One more comment: I've looked at the GDAL line of sight algorithm. It uses a distance to a pixel (which the line probably doesn't exactly intersect) to determine a height along the line from observer to target. This computation may well yield a Z value for a position on the line that isn't in the subject cell at all. There isn't anything wrong with this, but understand that determinations about observability are based on estimates of various sorts because a raster is not a continuous surface.

I can't comment on any other algorithms because I haven't studied them, but they necessarily are making decisions based on estimates.