2

My goal: Using a .TIF DEM, I would like to conduct a viewshed analysis. This is to be used as one feature for a larger program written in Python that will significantly improve work quality in my profession. Ultimately, I am looking for a binary output raster that I can manipulate with GeoPandas or Shapely, then plot as polygons in Folium.

Raster data viewshed analysis with Python was extremely helpful, but I have been unable to implement the suggested methods successfully. Below is my closest successful attempt using GDAL, but I have not found a specific example of the required inputs in the documentation.

from osgeo import gdal

band = gdal.Open('elevation.TIF').GetRasterBand(1)

gdal.ViewshedGenerate(
    srcBand = band,
    driverName = 'GTiff',
    targetRasterName = 'viewtest.TIF',
    creationOptions = 'COMPRESS=NONE',
    observerX = test_center[1],
    observerY = test_center[0],
    observerHeight = 1.9,
    targetHeight = 2,
    visibleVal = 1,
    invisibleVal = 0,
    outOfRangeVal = 0,
    noDataVal = 0,
    dfCurvCoeff = 1 - 1/7,
    mode = 'NORMAL',
    maxDistance = 1000)

I get the following error, related to the fourth line: creationOptions

File ~\anaconda3\envs\ARES\lib\site-packages\osgeo\gdal.py:4084, in ViewshedGenerate(*args, **kwargs)
   4082 def ViewshedGenerate(*args, **kwargs):
   4083     r"""ViewshedGenerate(Band srcBand, char const * driverName, char const * targetRasterName, char ** creationOptions, double observerX, double observerY, double observerHeight, double targetHeight, double visibleVal, double invisibleVal, double outOfRangeVal, double noDataVal, double dfCurvCoeff, GDALViewshedMode mode, double maxDistance, GDALProgressFunc callback=0, void * callback_data=None, GDALViewshedOutputType heightMode=GVOT_NORMAL, char ** papszOptions=None) -> Dataset"""
-> 4084     return _gdal.ViewshedGenerate(*args, **kwargs)

TypeError: in method 'ViewshedGenerate', argument 4 of type 'char **'

Here is the documentation I found related to creationOptions for drivers, but I was not able to discern how to implement it specifically in Python. I cannot find anywhere to describe what format needs to be passed into the function.

I am not invested in any specific method/GIS- I have also tried QGIS and GRASS APIs to no avail, and this package but I was unable to install it successfully.

2 Answers 2

3

When I set creationOptions=None I was able to avoid that particular error. However, I suspect you will have a problem with mode='NORMAL as that is what I'm currently stuck on.

EDIT: From the gdal_alg.h file contained somewhere in the package, I discovered that the GDALViewshedMode is an enumerated type with the following enumerations:

GVM_Diagonal = 1,
GVM_Edge = 2,
GVM_Max = 3,
GVM_Min = 4

So if you do fix the creationOptions issue, replace mode=Normal with a number 1 to 4 to get that option to work correctly.

3

The Python API documentation for the ViewshedGenerate algorithm is at https://gdal.org/api/python/osgeo.gdal.html#osgeo.gdal.ViewshedGenerate

osgeo.gdal.ViewshedGenerate(Band srcBand, char const * driverName, char const * targetRasterName, char ** creationOptions, double observerX, double observerY, double observerHeight, double targetHeight, double visibleVal, double invisibleVal, double outOfRangeVal, double noDataVal, double dfCurvCoeff, GDALViewshedMode mode, double maxDistance, GDALProgressFunc callback=0, void * callback_data=None, GDALViewshedOutputType heightMode=GVOT_NORMAL, char ** options=None)→ Dataset

A complete Python example can be found from the GDAL autotest suite https://github.com/OSGeo/gdal/blob/master/autotest/utilities/test_gdal_viewshed.py#L161

Here as slightly modified so it can be run with copy-paste after the input parameters are edited to match the real data.

from osgeo import gdal

ox=[292294]
oy=[6571861]
oz=[2]
src_ds = gdal.Open('dem.tif')
ds = gdal.ViewshedGenerate(
    src_ds.GetRasterBand(1),
    "GTiff",
    "outtest.tif",
    ["INTERLEAVE=BAND"],
    ox[0],
    oy[0],
    oz[0],
    0,  # targetHeight
    255,  # visibleVal
    0,  # invisibleVal
    0,  # outOfRangeVal
    -1.0,  # noDataVal,
    0.85714,  # dfCurvCoeff
    gdal.GVM_Edge,
    0,  # maxDistance
    heightMode=gdal.GVOT_MIN_TARGET_HEIGHT_FROM_GROUND,
    options=["UNUSED=YES"],
)
2
  • what are the units of ox=[292294] and oy=[6571861]. If my raster is in geographic CS, must I reproject my raster and my location ox and oy before performing viewshed or can I pass them directly as lon and lat and let raster as it is?
    – Khaled
    Jan 13 at 22:06
  • From the documentation gdal.org/programs/gdal_viewshed.html: "The algorithm as implemented currently will only output meaningful results if the georeferencing is in a projected coordinate reference system." Reproject your DEM into some suitable projected system and give ox and oy in the units of that system
    – user30184
    Jan 14 at 10:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.