Calculating limits of visible horizon of location in QGIS?

Question

I am an amateur at GIS. I installed QGIS Las Palmas version in Windows, I searched YouTube and also in Google but I didn't find this question.

There is a website that one can input the coordinates of a location and its elevation above ground level, then the website calculates the viewshed (cloak visibility) and the horizon limit as a kml polyline.

Having the DEM (.hgt file) for the desired region, how can I calculate the real horizon limit (as .kml polyline/polygon), not the viewshed, taking into account terrain in QGIS, just the limit of visible horizon (that it would be theoretically a circular polyline around the location without accounting for terrain).

Such a thing:

I would have such a horizon with QGIS, because the rendered horizon from the the cited website is not for each compass bearing, I am looking for a more accurate horizon line.

My goal is to compute at each bearing, the minimum angular elevation of sun or other objects in order to be seen from the location of observer.

Answer 1

WIthin QGIS you can use Grass tool r.horizon.height to calculate sun height for desired bearings. Result is raster map representing horizon height for each point. However, since you need this for each bearing (I presume 360?) this would mean you would have to generate 360 maps.

One way is by using Viewshed Analysis plugin avalible in QGIS which can be used to calculate horizon. If you vectorize result of that you could calculate height angles in combination with DEM.

Answer 2

The output of the 'website' you mention is generated using 90 meter srtm. You can use easily and quickly use output from the 'website' to generate 'Stellarium' horizons/landscapes.

There are other methods for generating landscapes using 30 meter srtm.

Have a look at the Stellarium astronomy program and one of the sample landscapes and assess if you had equivalent landscapes for your areas of interest whether you can pull the info you need back from Stellarium.

If this meets your requirements repost and I'll explain how to do it.

Answer 3

I wrote up some code that calculates the visible horizon as a polygon with 360 vertices

The output is a little rough but it works

from pathlib import Path

"""
##########################################################
User options
"""

#Variable assignment
inDEM        = 'C:/DEM.tif'
viewDist     = 3000
easting = '523926'
northing = '5258156'
crs = '[EPSG:28355]'

compressOptions         = 'COMPRESS=ZSTD|NUM_THREADS=ALL_CPUS|PREDICTOR=1|ZSTD_LEVEL=1|BIGTIFF=IF_SAFER|TILED=YES'

"""
##########################################################
Variable assignment for processing
"""

#Get the location of the initial image for storage of processing files
rootProcessDirectory = str(Path(inDEM).parent.absolute()).replace('\\','/') + '/'

#Set up the layer name for the raster calculations
inDEMName = inDEM.split("/")
inDEMName = inDEMName[-1]
inDEMName = inDEMName[:len(inDEMName)-4]
inDEMName = inDEMName[:8]
demRas = QgsRasterLayer(inDEM) 
pixelSize = (demRas.rasterUnitsPerPixelX() + demRas.rasterUnitsPerPixelY())/2

processDirectory = rootProcessDirectory + inDEMName + 'Process' + '/'

if not os.path.exists(processDirectory): os.mkdir(processDirectory)

"""
##########################################################
Processing
"""

processing.run("gdal:viewshed", {'INPUT':inDEM,'BAND':1,'OBSERVER':easting + ', ' + northing + ' ' + crs,'OBSERVER_HEIGHT':1.5,'TARGET_HEIGHT':1.5,'MAX_DISTANCE':viewDist,
'OPTIONS':compressOptions,'EXTRA':'','OUTPUT':processDirectory + 'Viewshed.tif'})

processing.run("native:pixelstopoints", {'INPUT_RASTER':processDirectory + 'Viewshed.tif','RASTER_BAND':1,'FIELD_NAME':'VALUE','OUTPUT':processDirectory + 'ViewshedPoints.gpkg'})

processing.run("native:extractbyexpression", {'INPUT':processDirectory + 'ViewshedPoints.gpkg','EXPRESSION':' \"VALUE\"  =  255 ',
'OUTPUT':processDirectory + 'ViewshedPointsExtract.gpkg'})

processing.run("native:rastersampling", {'INPUT':processDirectory + 'ViewshedPointsExtract.gpkg','RASTERCOPY':inDEM,'COLUMN_PREFIX':'DEM','OUTPUT':processDirectory + 'ViewshedPointsExtractDEM.gpkg'})

processing.run("native:extractbyexpression", {'INPUT':processDirectory + 'ViewshedPointsExtractDEM.gpkg','EXPRESSION':' \"DEM1\" IS NOT NULL ', 'OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtract.gpkg'})

processing.run("native:fieldcalculator", {'INPUT':processDirectory + 'ViewshedPointsExtractDEMExtract.gpkg','FIELD_NAME':'Angle','FIELD_TYPE':1,'FIELD_LENGTH':0,'FIELD_PRECISION':0,
'FORMULA':'angle_at_vertex( make_line(  make_point(' + easting + ',' + northing + '), $geometry ),0)','OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngle.gpkg'})

processing.run("native:fieldcalculator", {'INPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngle.gpkg','FIELD_NAME':'MaxDistance','FIELD_TYPE':1,'FIELD_LENGTH':0,'FIELD_PRECISION':0,
'FORMULA':'length( make_line( make_point(' + easting + ',' + northing + '), $geometry ))','OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDist.gpkg'})

processing.run("native:aggregate", {'INPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDist.gpkg','GROUP_BY':'\"Angle\"',
'AGGREGATES':[{'aggregate': 'maximum','delimiter': ',','input': '\"fid\"','length': 0,'name': 'fid','precision': 0,'type': 4},{'aggregate': 'maximum','delimiter': ',',
'input': '\"VALUE\"','length': 0,'name': 'VALUE','precision': 0,'type': 6},{'aggregate': 'maximum','delimiter': ',','input': '\"Angle\"','length': 0,'name': 'Angle','precision': 0,'type': 2},
{'aggregate': 'maximum','delimiter': ',','input': '\"MaxDistance\"','length': 0,'name': 'MaxDistance','precision': 0,'type': 2}],
'OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregate.gpkg'})

processing.run("native:multiparttosingleparts", {'INPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregate.gpkg','OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingle.gpkg'})

processing.run("native:fieldcalculator", {'INPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingle.gpkg','FIELD_NAME':'Distance','FIELD_TYPE':1,'FIELD_LENGTH':0,'FIELD_PRECISION':0,
'FORMULA':'length( make_line( make_point(' + easting + ',' + northing + '), $geometry ))','OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingleWDist.gpkg'})

processing.run("native:extractbyexpression", {'INPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingleWDist.gpkg','EXPRESSION':' \"Distance\"  =  \"MaxDistance\" ',
'OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingleWDistMaxed.gpkg'})

processing.run("qgis:pointstopath", {'INPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingleWDistMaxed.gpkg','CLOSE_PATH':True,'ORDER_FIELD':'Angle','GROUP_FIELD':'','DATE_FORMAT':'',
'OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingleWDistMaxedPath.gpkg'})

processing.runAndLoadResults("qgis:linestopolygons", {'INPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingleWDistMaxedPath.gpkg'
,'OUTPUT':processDirectory + 'ViewshedPointsExtractDEMExtractWAngleWDistAggregateSingleWDistMaxedPathPolygon.gpkg'})