# -*- coding: utf-8 -*-
import cupy as cp
import gc
"""
/***************************************************************************
ProcessingUMEP
A QGIS plugin
UMEP for processing toolbox
Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
-------------------
begin : 2020-04-02
copyright : (C) 2020 by Fredrik Lindberg
email : fredrikl@gvc.gu.se
Modifications by:
Jessica Monahan 2025
Modified so code can run without QGIS, and can have 3D input
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
"""
__author__ = 'Fredrik Lindberg'
__date__ = '2020-04-02'
__copyright__ = '(C) 2020 by Fredrik Lindberg'
# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'
from osgeo import gdal
import os
import numpy as np
from ..util import misc
from ..functions import svf_functions as svf
[docs]
class ProcessingSkyViewFactorAlgorithm():
"""
A class to process the SkyViewFactor (SVF) algorithm inputs.
Parameters:
INPUT_DSM (str): Digital Surface Model (DSM) input file.
INPUT_CDSM (str): Canopy DSM input file.
INPUT_TDSM (str): Trunk zone DSM input file.
OUTPUT_DIR (str): Output directory.
OUTPUT_FILE (str): Output file.
TRANS_VEG (float or str): Vegetation transmissivity value or input file.
INPUT_THEIGHT (float): Trunk height (percentage, if no trunk zone DSM).
use_veg (bool): Create SVF with trees.
tdsm_exist (bool): Does the trunk DSM exist.
ANISO (bool): Create output for anisotropic sky.
INPUT_EXTRAHEIGHT (float): Highest possible height of the z-component of the solar vector,
dependent on max solar height for location.
INPUT_DTM (str): Digital Terrain Model (DTM) input file.
MULT_DSMS (str): For 3D SOLWEIG: input of layered DSM file.
Outputs:
The class produces and stores the following files:
1. SVF files stored as `.tif` for buildings, and, if chosen, vegetation separately.
2. Directional SVF files stored as `.tif` for north, east, south, and west, for buildings and optionally vegetation separately.
3. Shadow matrices saved as `.npz` files for all shadow computations involved.
"""
def __init__(self, INPUT_DSM, INPUT_CDSM, OUTPUT_DIR, OUTPUT_FILE, INPUT_DTM=None, INPUT_EXTRAHEIGHT=6, INPUT_MULT_DSMS = None, INPUT_TDSM=None, USE_VEG=True, TRANS_VEG=15,
TSDM_EXIST=False, INPUT_THEIGHT=25.0, ANISO=True):
self.INPUT_DSM = INPUT_DSM
self.INPUT_CDSM = INPUT_CDSM
self.INPUT_TDSM = INPUT_TDSM
self.INPUT_DTM = INPUT_DTM
self.USE_VEG = USE_VEG
self.TRANS_VEG = TRANS_VEG
self.TSDM_EXIST = TSDM_EXIST
self.INPUT_THEIGHT = INPUT_THEIGHT
self.INPUT_EXTRAHEIGHT = INPUT_EXTRAHEIGHT
self.MULT_DSMS = INPUT_MULT_DSMS
self.ANISO = ANISO
self.OUTPUT_DIR = OUTPUT_DIR
self.OUTPUT_FILE = OUTPUT_FILE
[docs]
def processAlgorithm(self):
"""
This algorithm is a processing version of SkyViewFactor
"""
# InputParameters
outputDir = self.OUTPUT_DIR
outputFile = self.OUTPUT_FILE
dsm_path = self.INPUT_DSM
# usevegdem = self.USE_VEG
transVeg = float(self.TRANS_VEG)
vegdsm_path = self.INPUT_CDSM
vegdsm2_path = self.INPUT_TDSM
# tdsmExists = self.parameterAsBool(parameters, self.TSDM_EXIST, context.dxf)
trunkr = float(self.INPUT_THEIGHT)
aniso = bool(self.ANISO)
dtm_path = self.INPUT_DTM
print('Initiating algorithm')
if dtm_path is not None:
gdal_dtm = gdal.Open(dtm_path)
dtm = cp.array(gdal_dtm.ReadAsArray().astype(float), dtype=cp.float32)
else:
dtm = None
gdal_dsm = gdal.Open(dsm_path)
dsm = cp.array(gdal_dsm.ReadAsArray(), dtype=cp.float32)
# response to issue #85
nd = gdal_dsm.GetRasterBand(1).GetNoDataValue()
dsm[dsm == nd] = 0.
dsm_min = dsm.min()
print(dsm_min)
if 0 <= dsm_min < self.INPUT_EXTRAHEIGHT:
dsmraise = self.INPUT_EXTRAHEIGHT - dsm_min
elif dsm_min < 0:
dsmraise = cp.abs(dsm.min()) + self.INPUT_EXTRAHEIGHT
else:
dsmraise = 0
dsm += dsmraise
print('DSM raised with ' + str(dsmraise) + 'm.')
if dtm_path is not None:
dtm += dsmraise
sizex = dsm.shape[0]
sizey = dsm.shape[1]
geotransform = gdal_dsm.GetGeoTransform()
scale = 1 / geotransform[1]
trans = transVeg / 100.0
if vegdsm_path:
usevegdem = 1
print('Vegetation scheme activated')
# vegdsm = self.parameterAsRasterLayer(parameters, self.INPUT_CDSM, context.dxf)
# if vegdsm is None:
# raise QgsProcessingException("Error: No valid vegetation DSM selected")
# load raster
gdal_vegdsm = gdal.Open(vegdsm_path)
vegdsm = cp.array(gdal_vegdsm.ReadAsArray(), dtype=cp.float32)
vegdsm = cp.where(vegdsm <= 0, np.nan, vegdsm)
vegsizex = vegdsm.shape[0]
vegsizey = vegdsm.shape[1]
if not (vegsizex == sizex) & (vegsizey == sizey):
raise Exception("Error in Vegetation Canopy DSM: All rasters must be of same extent and resolution")
if vegdsm2_path:
# vegdsm2 = self.parameterAsRasterLayer(parameters, self.INPUT_TDSM, context.dxf)
# if vegdsm2 is None:
# raise QgsProcessingException("Error: No valid Trunk zone DSM selected")
# load raster
gdal_vegdsm2 = gdal.Open(vegdsm2_path)
vegdsm2 = cp.array(gdal_vegdsm2.ReadAsArray(), dtype=cp.float32)
vegdsm2 = cp.where(vegdsm2 <= 0, np.nan, vegdsm2)
else:
trunkratio = trunkr / 100.0
vegdsm2 = vegdsm * trunkratio
vegsizex = vegdsm2.shape[0]
vegsizey = vegdsm2.shape[1]
if not (vegsizex == sizex) & (vegsizey == sizey):
raise Exception("Error in Trunk Zone DSM: All rasters must be of same extent and resolution")
else:
rows = dsm.shape[0]
cols = dsm.shape[1]
vegdsm = np.zeros([rows, cols])
vegdsm2 = 0.
usevegdem = 0
if aniso == 1:
print('Calculating SVF using 153 iterations')
# ret = svf.svfForProcessing153(dsm, vegdsm, vegdsm2, scale, usevegdem)
ret = svf.svfForProcessing153(dsm, dtm, vegdsm, vegdsm2, scale, usevegdem)
else:
print('Calculating SVF using 655 iterations')
ret = svf.svfForProcessing655(dsm, dtm, vegdsm, vegdsm2, scale, usevegdem)
filename = outputFile
if not os.path.exists(outputDir + '/svfs'):
os.makedirs(outputDir + '/svfs')
if ret is not None:
svfbu = ret["svf"].get()
svfbuE = ret["svfE"].get()
svfbuS = ret["svfS"].get()
svfbuW = ret["svfW"].get()
svfbuN = ret["svfN"].get()
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svf.tif', svfbu)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfE.tif', svfbuE)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfS.tif', svfbuS)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfW.tif', svfbuW)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfN.tif', svfbuN)
if usevegdem == 0:
svftotal = svfbu
else:
# report the result
svfveg = ret["svfveg"].get()
svfEveg = ret["svfEveg"].get()
svfSveg = ret["svfSveg"].get()
svfWveg = ret["svfWveg"].get()
svfNveg = ret["svfNveg"].get()
svfaveg = ret["svfaveg"].get()
svfEaveg = ret["svfEaveg"].get()
svfSaveg = ret["svfSaveg"].get()
svfWaveg = ret["svfWaveg"].get()
svfNaveg = ret["svfNaveg"].get()
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfveg.tif', svfveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfEveg.tif', svfEveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfSveg.tif', svfSveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfWveg.tif', svfWveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfNveg.tif', svfNveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfaveg.tif', svfaveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfEaveg.tif', svfEaveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfSaveg.tif', svfSaveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfWaveg.tif', svfWaveg)
misc.saveraster(gdal_dsm, outputDir + '/svfs/' + 'svfNaveg.tif', svfNaveg)
trans = transVeg / 100.0
svftotal = (svfbu - (1 - svfveg) * (1 - trans))
misc.saveraster(gdal_dsm, filename, svftotal)
# Save shadow images for SOLWEIG 2019a
if aniso == 1:
shmat = ret["shmat"].get()
vegshmat = ret["vegshmat"].get()
vbshvegshmat = ret["vbshvegshmat"].get()
# wallshmat = ret["wallshmat"]
# wallsunmat = ret["wallsunmat"]
# wallshvemat = ret["wallshvemat"]
# facesunmat = ret["facesunmat"]
np.savez_compressed(outputDir + '/' + "shadowmats.npz", shadowmat=shmat, vegshadowmat=vegshmat,
vbshmat=vbshvegshmat)
# np.savez_compressed(outputDir + '/' + "shadowmats.npz", shadowmat=shmat, vegshadowmat=vegshmat, vbshmat=vbshvegshmat)
# vbshvegshmat=vbshvegshmat, wallshmat=wallshmat, wallsunmat=wallsunmat,
# facesunmat=facesunmat, wallshvemat=wallshvemat)
print("Sky View Factor: SVF grid(s) successfully generated")
gc.collect()
# Clear memory pools
cp._default_memory_pool.free_all_blocks()
cp._default_pinned_memory_pool.free_all_blocks()
return {self.OUTPUT_DIR: outputDir, self.OUTPUT_FILE: outputFile}
[docs]
def processAlgorithm_3d(self):
"""
This algorithm is a 3D processing version of SkyViewFactor
"""
# InputParameters
outputDir = self.OUTPUT_DIR
outputFile = self.OUTPUT_FILE
dsms_path = self.MULT_DSMS
# usevegdem = self.USE_VEG
transVeg = float(self.TRANS_VEG)
vegdsm_path = self.INPUT_CDSM
vegdsm2_path = self.INPUT_TDSM
# tdsmExists = self.parameterAsBool(parameters, self.TSDM_EXIST, context.dxf)
trunkr = float(self.INPUT_THEIGHT)
aniso = bool(self.ANISO)
dtm_path = self.INPUT_DTM
if dtm_path is not None:
gdal_dtm = gdal.Open(dtm_path)
dtm = cp.array(gdal_dtm.ReadAsArray().astype(float), dtype=cp.float32)
else:
dtm = None
print('Initiating algorithm')
gdal_dsms = gdal.Open(dsms_path)
layers = gdal_dsms.RasterCount
dsms = cp.stack([cp.array(gdal_dsms.GetRasterBand(i).ReadAsArray(), dtype=cp.float32) for i in range(1, layers + 1)],
axis=0)
dsm_min = dsms[0].min()
if 0 <= dsm_min < self.INPUT_EXTRAHEIGHT:
dsmraise = self.INPUT_EXTRAHEIGHT - dsm_min
elif dsm_min < 0:
dsmraise = cp.abs(dsm_min) + self.INPUT_EXTRAHEIGHT
else:
dsmraise = 0
dsms += dsmraise
if dtm_path is not None:
dtm += dsmraise
sizex = dsms[0].shape[0]
sizey = dsms[0].shape[1]
geotransform = gdal_dsms.GetGeoTransform()
scale = 1 / geotransform[1]
trans = transVeg / 100.0
if vegdsm_path is not None:
usevegdem = 1
print('Vegetation scheme activated')
# load raster
gdal_vegdsm = gdal.Open(vegdsm_path)
vegdsm = cp.array(gdal_vegdsm.ReadAsArray(), dtype=cp.float32)
vegdsm = cp.where(vegdsm <= 0, np.nan, vegdsm)
vegsizex = vegdsm.shape[0]
vegsizey = vegdsm.shape[1]
if not (vegsizex == sizex) & (vegsizey == sizey):
raise Exception("Error in Vegetation Canopy DSM: All rasters must be of same extent and resolution")
if vegdsm2_path:
# load raster
gdal_vegdsm2 = gdal.Open(vegdsm2_path)
vegdsm2 = cp.array(gdal_vegdsm2.ReadAsArray(), dtype=cp.float32)
vegdsm2 = cp.where(vegdsm2 <= 0, np.nan, vegdsm2)
else:
trunkratio = trunkr / 100.0
vegdsm2 = vegdsm * trunkratio
vegsizex = vegdsm2.shape[0]
vegsizey = vegdsm2.shape[1]
if not (vegsizex == sizex) & (vegsizey == sizey):
raise Exception("Error in Trunk Zone DSM: All rasters must be of same extent and resolution")
else:
rows = dsms[0].shape[0]
cols = dsms[0].shape[1]
vegdsm = cp.zeros([rows, cols])
vegdsm2 = 0.
usevegdem = 0
if aniso == 1:
print('Calculating SVF using 153 iterations')
ret = svf.svfForProcessing153_3d(dsms, dtm, vegdsm, vegdsm2, scale, usevegdem)
else:
print('Calculating SVF using 655 iterations')
ret = svf.svfForProcessing655_3d(dsms, dtm, vegdsm, vegdsm2, scale, usevegdem)
filename = outputFile
if not os.path.exists(outputDir + '/svfs'):
os.makedirs(outputDir + '/svfs')
if ret is not None:
svfbu = ret["svf"].get()
svfbuE = ret["svfE"].get()
svfbuS = ret["svfS"].get()
svfbuW = ret["svfW"].get()
svfbuN = ret["svfN"].get()
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svf.tif', svfbu)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfE.tif', svfbuE)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfS.tif', svfbuS)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfW.tif', svfbuW)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfN.tif', svfbuN)
if os.path.isfile(outputDir + '/' + 'svfs.zip'):
os.remove(outputDir + '/' + 'svfs.zip')
if usevegdem == 0:
svftotal = svfbu
else:
# report the result
svfveg = ret["svfveg"].get()
svfEveg = ret["svfEveg"].get()
svfSveg = ret["svfSveg"].get()
svfWveg = ret["svfWveg"].get()
svfNveg = ret["svfNveg"].get()
svfaveg = ret["svfaveg"].get()
svfEaveg = ret["svfEaveg"].get()
svfSaveg = ret["svfSaveg"].get()
svfWaveg = ret["svfWaveg"].get()
svfNaveg = ret["svfNaveg"].get()
misc.saveraster(gdal_vegdsm, outputDir + '/svfs/' + 'svfveg.tif', svfveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfEveg.tif', svfEveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfSveg.tif', svfSveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfWveg.tif', svfWveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfNveg.tif', svfNveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfaveg.tif', svfaveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfEaveg.tif', svfEaveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfSaveg.tif', svfSaveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfWaveg.tif', svfWaveg)
misc.saveraster(gdal_dsms, outputDir + '/svfs/' + 'svfNaveg.tif', svfNaveg)
trans = transVeg / 100.0
svftotal = (svfbu - (1 - svfveg) * (1 - trans))
misc.saveraster(gdal_dsms, filename, svftotal)
# Save shadow images for SOLWEIG 2019a
if aniso == 1:
shmat = ret["shmat"].get()
vegshmat = ret["vegshmat"].get()
vbshvegshmat = ret["vbshvegshmat"].get()
np.savez_compressed(outputDir + '/' + "shadowmats.npz", shadowmat=shmat, vegshadowmat=vegshmat,
vbshmat=vbshvegshmat) # ,
# vbshvegshmat=vbshvegshmat, wallshmat=wallshmat, wallsunmat=wallsunmat,
# facesunmat=facesunmat, wallshvemat=wallshvemat)
print("Sky View Factor: SVF grid(s) successfully generated")
gc.collect()
# Clear memory pools
cp._default_memory_pool.free_all_blocks()
cp._default_pinned_memory_pool.free_all_blocks()
return {self.OUTPUT_DIR: outputDir, self.OUTPUT_FILE: outputFile}