0

I am doing a viewshed analysis on a DEM .tiff file. I would like to convert the viewshed raster into vector geometry in the format of .geojson, which i have achieved, although, the result is not good. There are big "in sight" areas which are filled with singular "out of sight" points, and vise-versa. These singular alonestanding points disturbs and complicates the image.

This is my current result after converting the raster image to a vector layer: enter image description here

enter image description here

I have tried some processing tools that simplify the geometry in some way, but with no good result. This is a try with the "Vector Geometry / Simplify" tool:

enter image description here

Not a good result...

These are the steps i would imagine a tool suitable for this problem should do:

  1. Merge features smaller than a given size
  2. Simplify edges of the remaining geometries (specifically, remove points that are coincident between two neighboring points iteratively.)

The final data will hopefully be not that many features with relatively simple geometry with not so many points. My intent will eventually be to send the data over a HTTP server.

I have looked into the QGIS processing tools, but I cant seem to find any appropriate tool for solving my case. How can I simplify and "clean" the analysed shedview vector layer?

My final result will be code in Python, preferably using the osgeo library, maybe also other libraries, maybe geopandas.

Improve this question
2
  • 1
    Maybe clean up isolated pixels before rasterising gis.stackexchange.com/q/483056/2856 Commented Jul 25 at 0:01
  • 1
    I can confirm that the processing tool "Sieve" helps alot as a first "filter", thank you! I'll keep you updated on the rest of the process. Commented Jul 25 at 15:57

3 Answers 3

Reset to default
0

If you want to clean up the vector data after viewshed analysis:

It might be better to prepare the raster data beforehand to reduce the amount of noise.

Improve this answer
0

You can select small polygons by expression and eliminate them into larger adjacent ones.

enter image description here

I want to eliminate polygons with an area less than 2000: $area<2000

enter image description here

Improve this answer
0

I have now solved the problem as i initially intended to solve it. These are the processing steps:

  1. Perform viewshed analysis on raster image (GDAL/Raster miscellaneous/Viewshed)
  2. Use the "Sieve" tool to clean up image (GDAL/Raster analysis/Sieve)
  3. Convert the raster images to polygons (GDAL/Raster conversion/Polygonize (raster to vector)
  4. Simplify the feature geometry (Vector geometry/Simplify)

A notable caveat in this process is:

  • for Sieve to work, both the visible (typically valued 255) and the non-visible (typically valued 0) points must exist on the image. You cannot remove the non-visible points from the image.
  • Conversely, when converting the image from raster to vector, the non-visible poitns must be removed in order to only create polygons from the visible points.

Therefore, the non-visible points must be removed/hidden between step 2. and 3.

I did also manage to implement this in python code:

from osgeo import gdal, ogr, osr from uuid import uuid1 import json filepath = "./static/FO_DSM_2015_FOTM_10M.tif" async def getViewshed(lon: float, lat: float): id = uuid1() drv_geojson: gdal.Driver = ogr.GetDriverByName("GeoJSON") drv_gtiff: gdal.Driver = gdal.GetDriverByName("GTiff") # Open the TIFF file using GDAL dem_dataset: gdal.Dataset = gdal.Open(filepath) # cast(Dataset | None, gdal.Open(filepath)) if not dem_dataset: raise FileNotFoundError(f"Unable to open file: {filepath}") # Read the raster band (assuming the height/elevation is in the first band) dem_band: gdal.Band = dem_dataset.GetRasterBand(1) # cast(Band | None, dataset.GetRasterBand(1)) if not dem_band: raise ValueError("No raster band found in the file.") projection: str = dem_dataset.GetProjection() geotransform = dem_dataset.GetGeoTransform() # 1. Viewshed viewshed_dataset: gdal.Dataset = gdal.ViewshedGenerate( srcBand = dem_band, driverName = 'GTiff', targetRasterName = f'./static/{id}_viewshed.tif', creationOptions = ['COMPRESS=LZW'], observerX = lon, observerY = lat, observerHeight = 2.0, targetHeight = 0, visibleVal = 255.0, invisibleVal = 0.0, outOfRangeVal = 0.0, noDataVal = 1.0, dfCurvCoeff = 1.0, mode = 1, maxDistance = 120000.0 ) viewshed_band: gdal.Band = viewshed_dataset.GetRasterBand(1) # 2. Sieve sieve_dataset: gdal.Dataset = drv_gtiff.Create( utf8_path = f"./static/{id}_sieve.tif", xsize = dem_dataset.RasterXSize, ysize = dem_dataset.RasterYSize, bands = 1, eType = dem_band.DataType ) sieve_dataset.SetGeoTransform(geotransform) sieve_band: gdal.Band = sieve_dataset.GetRasterBand(1) gdal.SieveFilter( srcBand = viewshed_band, maskBand = None, dstBand = sieve_band, threshold = 25, connectedness = 4, ) sieve_band.SetNoDataValue(0.0) sieve_band.FlushCache() sieve_mask: gdal.Band = sieve_band.GetMaskBand() sp_ref = osr.SpatialReference() sp_ref.ImportFromWkt(projection) viewshed_geojson = drv_geojson.CreateDataSource(f"./static/{id}_vector.geojson") viewshed_layer: ogr.Layer = viewshed_geojson.CreateLayer("viewshed", srs = sp_ref) field_definition = ogr.FieldDefn("HA", ogr.OFTInteger) viewshed_layer.CreateField(field_definition) field_index = viewshed_layer.GetLayerDefn().GetFieldIndex("HA") gdal.Polygonize( srcBand = sieve_band, maskBand = sieve_mask, outLayer = viewshed_layer, iPixValField = field_index ) viewshed_geojson = None viewshed_geojson = ogr.Open(f"./static/{id}_vector.geojson", 1) viewshed_layer = viewshed_geojson.GetLayer() print("Feature count:", viewshed_layer.GetFeatureCount()) for feature in viewshed_layer: geometry: ogr.Geometry = feature.GetGeometryRef() if geometry is not None: simplified_geometry: ogr.Geometry = geometry.Simplify(35) feature.SetGeometry(simplified_geometry) viewshed_layer.SetFeature(feature) viewshed_geojson = None with open(f"./static/{id}_vector.geojson", "r") as file: json_data = json.load(file) return json_data

When implementing this code, there were some gotchas. One of them were that processed data was not immediately available for further processing. I had to close and reopen the file in order to have the processed data available. This code can surely be cleaned up, although it works!

Although this solved my initial problem statement, the result is not that pretty with very pointy and edgy polygons. Therefore i am considering trying to find a solution based on heatmaps, which might give a more visual-pleasing result.

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.