sky_view_factor package added

urbantools/svf/code/create_patches.py

@@ -0,0 +1,52 @@
+import numpy as np
+
+def create_patches(patch_option):
+
+    deg2rad = np.pi/180
+
+    # patch_option = 1 = 145 patches (Robinson & Stone, 2004)
+    # patch_option = 2 = 153 patches (Wallenberg et al., 2022)
+    # patch_option = 3 = 306 patches -> test 
+    # patch_option = 4 = 612 patches -> test
+
+    skyvaultalt = np.atleast_2d([])
+    skyvaultazi = np.atleast_2d([])
+
+    # Creating skyvault of patches of constant radians (Tregeneza and Sharples, 1993)
+    # Patch option 1, 145 patches, Original Robinson & Stone (2004) after Tregenza (1987)/Tregenza & Sharples (1993)
+    if patch_option == 1:
+        annulino = np.array([0, 12, 24, 36, 48, 60, 72, 84, 90])
+        skyvaultaltint = np.array([6, 18, 30, 42, 54, 66, 78, 90]) # Robinson & Stone (2004)
+        azistart = np.array([0, 4, 2, 5, 8, 0, 10, 0]) # Fredrik/Nils
+        patches_in_band = np.array([30, 30, 24, 24, 18, 12, 6, 1]) # Robinson & Stone (2004)
+    # Patch option 2, 153 patches, Wallenberg et al. (2022)
+    elif patch_option == 2:
+        annulino = np.array([0, 12, 24, 36, 48, 60, 72, 84, 90])
+        skyvaultaltint = np.array([6, 18, 30, 42, 54, 66, 78, 90]) # Robinson & Stone (2004)
+        azistart = np.array([0, 4, 2, 5, 8, 0, 10, 0]) # Fredrik/Nils
+        patches_in_band = np.array([31, 30, 28, 24, 19, 13, 7, 1]) # Nils
+    # Patch option 3, 306 patches, test
+    elif patch_option == 3:
+        annulino = np.array([0, 12, 24, 36, 48, 60, 72, 84, 90])
+        skyvaultaltint = np.array([6, 18, 30, 42, 54, 66, 78, 90]) # Robinson & Stone (2004)
+        azistart = np.array([0, 4, 2, 5, 8, 0, 10, 0]) # Fredrik/Nils
+        patches_in_band = np.array([31*2, 30*2, 28*2, 24*2, 19*2, 13*2, 7*2, 1]) # Nils
+    # Patch option 4, 612 patches, test
+    elif patch_option == 4:
+        annulino = np.array([0, 4.5, 9, 15, 21, 27, 33, 39, 45, 51, 57, 63, 69, 75, 81, 90]) # Nils
+        skyvaultaltint = np.array([3, 9, 15, 21, 27, 33, 39, 45, 51, 57, 63, 69, 75, 81, 90]) # Nils
+        patches_in_band = np.array([31*2, 31*2, 30*2, 30*2, 28*2, 28*2, 24*2, 24*2, 19*2, 19*2, 13*2, 13*2, 7*2, 7*2, 1]) # Nils
+        azistart = np.array([0, 0, 4, 4, 2, 2, 5, 5, 8, 8, 0, 0, 10, 10, 0]) # Nils
+
+    skyvaultaziint = np.array([360/patches for patches in patches_in_band])
+
+    for j in range(0, skyvaultaltint.shape[0]):
+        for k in range(0, patches_in_band[j]):
+            skyvaultalt = np.append(skyvaultalt, skyvaultaltint[j])
+            skyvaultazi = np.append(skyvaultazi, k*skyvaultaziint[j] + azistart[j])
+
+    # skyvaultzen = (90 - skyvaultalt) * deg2rad
+    # skyvaultalt = skyvaultalt * deg2rad
+    # skyvaultazi = skyvaultazi * deg2rad
+
+    return skyvaultalt, skyvaultazi, annulino, skyvaultaltint, patches_in_band, skyvaultaziint, azistart

urbantools/svf/code/misc.py

@@ -0,0 +1,94 @@
+__author__ = 'xlinfr'
+
+import numpy as np
+from osgeo import gdal, osr
+from osgeo.gdalconst import GDT_Float32
+
+
+# Slope and aspect used in SEBE and Wall aspect
+def get_ders(dsm, scale):
+    # dem,_,_=read_dem_grid(dem_file)
+    dx = 1/scale
+    # dx=0.5
+    fy, fx = np.gradient(dsm, dx, dx)
+    asp, grad = cart2pol(fy, fx, 'rad')
+    slope = np.arctan(grad)
+    asp = asp * -1
+    asp = asp + (asp < 0) * (np.pi * 2)
+    return slope, asp
+
+
+def cart2pol(x, y, units='deg'):
+    radius = np.sqrt(x**2 + y**2)
+    theta = np.arctan2(y, x)
+    if units in ['deg', 'degs']:
+        theta = theta * 180 / np.pi
+    return theta, radius
+
+
+def saveraster(gdal_data, filename, raster):
+    rows = gdal_data.RasterYSize
+    cols = gdal_data.RasterXSize
+
+    outDs = gdal.GetDriverByName("GTiff").Create(filename, cols, rows, int(1), GDT_Float32)
+    outBand = outDs.GetRasterBand(1)
+
+    # write the data
+    outBand.WriteArray(raster, 0, 0)
+    # flush data to disk, set the NoData value and calculate stats
+    outBand.FlushCache()
+    outBand.SetNoDataValue(-9999)
+
+    # georeference the image and set the projection
+    outDs.SetGeoTransform(gdal_data.GetGeoTransform())
+    outDs.SetProjection(gdal_data.GetProjection())
+
+def saverasternd(gdal_data, filename, raster):
+    rows = gdal_data.RasterYSize
+    cols = gdal_data.RasterXSize
+
+    outDs = gdal.GetDriverByName("GTiff").Create(filename, cols, rows, int(1), GDT_Float32)
+    outBand = outDs.GetRasterBand(1)
+
+    # write the data
+    outBand.WriteArray(raster, 0, 0)
+    # flush data to disk, set the NoData value and calculate stats
+    outBand.FlushCache()
+    # outBand.SetNoDataValue(-9999)
+
+    # georeference the image and set the projection
+    outDs.SetGeoTransform(gdal_data.GetGeoTransform())
+    outDs.SetProjection(gdal_data.GetProjection())
+
+def xy2latlon(crsWtkIn, x, y):
+    old_cs = osr.SpatialReference()
+    old_cs.ImportFromWkt(crsWtkIn)
+
+    wgs84_wkt = """
+    GEOGCS["WGS 84",
+        DATUM["WGS_1984",
+            SPHEROID["WGS 84",6378137,298.257223563,
+                AUTHORITY["EPSG","7030"]],
+            AUTHORITY["EPSG","6326"]],
+        PRIMEM["Greenwich",0,
+            AUTHORITY["EPSG","8901"]],
+        UNIT["degree",0.01745329251994328,
+            AUTHORITY["EPSG","9122"]],
+        AUTHORITY["EPSG","4326"]]"""
+
+    new_cs = osr.SpatialReference()
+    new_cs.ImportFromWkt(wgs84_wkt)
+
+    transform = osr.CoordinateTransformation(old_cs, new_cs)
+
+    lonlat = transform.TransformPoint(x, y)
+
+    gdalver = float(gdal.__version__[0])
+    if gdalver >= 3.:
+        lon = lonlat[1] #changed to gdal 3
+        lat = lonlat[0] #changed to gdal 3
+    else:
+        lon = lonlat[0] #changed to gdal 2
+        lat = lonlat[1] #changed to gdal 2
+
+    return lat, lon

urbantools/svf/code/run_svf_standalone.py

@@ -0,0 +1,156 @@
+"""
+Standalone script for calculating the Sky View Factor (SVF) from a DSM in ASCII (.ASC) format.
+The code corrects the header of the .ASC file (converting XLLCENTER/YLLCENTER to XLLCORNER/YLLCORNER), 
+loads the raster, performs the SVF calculation with the UMEP algorithm (svfForProcessing153 or svfForProcessing655), 
+and finally saves the result in GeoTIFF and CSV format.
+
+Optionally, it generates a heatmap of the result.
+
+Author: Aldo, UMEP
+Dependencies: numpy, osgeo.gdal, pandas, matplotlib, UMEP functions (svf_functions.py)
+"""
+
+import numpy as np
+from osgeo import gdal
+import os
+import sys
+import matplotlib.pyplot as plt
+import pandas as pd
+# === Aggiungi percorso funzioni UMEP ===
+sys.path.append(os.path.join(os.getcwd(), "functions"))
+sys.path.append(os.path.join(os.getcwd(), "util"))
+import svf_functions as svf  # usa svfForProcessing153 o svfForProcessing655
+
+tile = "dsm_demo"
+
+# === Parametri ===
+input_asc_path = f"data/{tile}.ASC"       # <--- Inserisci il tuo DSM ASCII qui
+corrected_asc_path = f"data/corrected_dsm_{tile}.asc"
+output_svf_path = f"output/svf_{tile}.tif"
+
+ANISOTROPIC = True  # True = 153 direzioni, False = 655 direzioni
+TRANS_VEG = 0       # percentuale trasmissività vegetazione (0 se non presente)
+
+# === Funzione per correggere intestazione .ASC ===
+def fix_asc_header(input_path, output_path):
+    with open(input_path, 'r') as f:
+        lines = f.readlines()
+
+    header = {}
+    data_start_idx = 0
+    for i, line in enumerate(lines):
+        if line.strip().split()[0].upper() in ['NCOLS', 'NROWS', 'CELLSIZE', 'XLLCENTER', 'YLLCENTER', 'NODATA_VALUE']:
+            key, val = line.strip().split()
+            header[key.upper()] = float(val)
+        else:
+            data_start_idx = i
+            break
+
+    if 'XLLCENTER' in header and 'YLLCENTER' in header:
+        header['XLLCORNER'] = header['XLLCENTER'] - header['CELLSIZE'] / 2
+        header['YLLCORNER'] = header['YLLCENTER'] - header['CELLSIZE'] / 2
+
+    # Ricostruisci intestazione
+    new_header = [
+        f"NCOLS         {int(header['NCOLS'])}",
+        f"NROWS         {int(header['NROWS'])}",
+        f"XLLCORNER     {header['XLLCORNER']:.4f}",
+        f"YLLCORNER     {header['YLLCORNER']:.4f}",
+        f"CELLSIZE      {header['CELLSIZE']:.4f}",
+        f"NODATA_VALUE  {int(header['NODATA_VALUE'])}"
+    ]
+
+    with open(output_path, 'w') as f:
+        f.write('\n'.join(new_header) + '\n')
+        f.writelines(lines[data_start_idx:])
+
+    print(f"[INFO] Intestazione .ASC corretta: {output_path}")
+
+# === 1. Correggi intestazione .ASC ===
+fix_asc_header(input_asc_path, corrected_asc_path)
+
+# === 2. Carica DSM ===
+ds = gdal.Open(corrected_asc_path)
+dsm = ds.ReadAsArray().astype(float)
+
+nd = ds.GetRasterBand(1).GetNoDataValue()
+if nd is not None:
+    dsm[dsm == nd] = 0
+
+if dsm.min() < 0:
+    dsm += abs(dsm.min())
+
+gt = ds.GetGeoTransform()
+pixel_res = gt[1]
+scale = 1 / pixel_res
+rows, cols = dsm.shape
+
+# === 3. Dummy array vegetazione ===
+veg = np.zeros_like(dsm)
+veg2 = 0.0
+useveg = 0
+
+
+# === 4. Calcolo SVF ===
+class FakeFeedback:
+    def isCanceled(self):
+        return False
+    def setProgressText(self, text):
+        print(f"[FEEDBACK] {text}")
+    def setProgress(self, value):
+        print(f"[PROGRESS] {value}%")
+
+print("[INFO] Calcolo SVF...")
+if ANISOTROPIC:
+    feedback = FakeFeedback()
+    result = svf.svfForProcessing153(dsm, veg, veg2, scale, useveg, pixel_res, False, None, feedback)
+else:
+    result = svf.svfForProcessing655(dsm, veg, veg2, scale, useveg, feedback=None)
+
+svf_total = result["svf"]
+
+# === 5. Salva output ===
+driver = gdal.GetDriverByName('GTiff')
+out_ds = driver.Create(output_svf_path, cols, rows, 1, gdal.GDT_Float32)
+out_ds.SetGeoTransform(gt)
+out_ds.SetProjection(ds.GetProjection())
+out_ds.GetRasterBand(1).WriteArray(svf_total)
+out_ds.GetRasterBand(1).SetNoDataValue(0)
+out_ds.FlushCache()
+
+print(f"[SUCCESS] SVF salvato in: {output_svf_path}")
+
+
+# === 6. Crea CSV con coordinate e SVF ===
+print("[INFO] Generazione CSV...")
+
+x0, dx, _, y0, _, dy = gt
+x_coords = x0 + dx * np.arange(cols)
+y_coords = y0 + dy * np.arange(rows)
+xx, yy = np.meshgrid(x_coords, y_coords)
+
+flat_coords = np.column_stack((xx.flatten(), yy.flatten()))
+flat_svf = svf_total.flatten()
+
+# Rimuovi i punti nodata (dove SVF = 0 o NaN)
+valid_mask = ~np.isnan(flat_svf) & (flat_svf > 0)
+coords_strings = [f"{x:.2f},{y:.2f}" for x, y in flat_coords[valid_mask]]
+
+df = pd.DataFrame({
+    "coordinate": coords_strings,
+    "svf": flat_svf[valid_mask]
+})
+
+csv_output_path = f"output/svf_{tile}.csv"
+df.to_csv(csv_output_path, index=False)
+print(f"[SUCCESS] CSV salvato in: {csv_output_path}")
+
+
+print("[INFO] Visualizzazione Heatmap...")
+plt.figure(figsize=(10, 8))
+plt.imshow(svf_total, cmap='coolwarm', origin='upper')
+plt.colorbar(label="Sky View Factor")
+plt.title("SVF Heatmap")
+plt.axis("off")
+plt.tight_layout()
+plt.show()