Skip to content

Commit c87bb20

Browse files
dwolfschdwolfsch
committed
spatial Voellmy friction model: add VariableVoellmyShapeToRaster to generate rasters from shapes
Adds functionality for generating raster files for mu and xsi values from a DEM and shapefiles including: A main script to rasterize mu and xsi based on polygon shapefiles and attribute fields. A corresponding run script and configuration file for user customization.
1 parent bff5409 commit c87bb20

File tree

5 files changed

+224
-9
lines changed

5 files changed

+224
-9
lines changed

‎avaframe/com6RockAvalanche/readMeVariableVoellmyShapeToRaster.txt‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
The VariableVoellmyShapeToRaster.py script allows the user to define spatially different values for the voellmy Parameters mu and xsi, with the use of polygon shapefiles. For the extent of a DEM raster, all the areas that are not covered by a polygon get assigned a default mu or xsi value. The script then converts this Information into a raster mu and a raster xsi file, which can then be used in Avaframe Simulation runs, using the "spatialVoellmy" friction model.
2+
3+
First, set up the Config File and provide inputs:
4+
oInputs:
5+
dem: Path to the DEM Raster that is later on used for the avaframe simulation. This is needed, because the mu and xsi output rasters need to be the exact same size. Usually this file lies in \AvaFrame\avaframe\data\ *yourAvalancheDir*\Inputs
6+
mu_shapefile: Path to the mu shapefile, that is then converted to a raster file. Be aware, that the attribute has to be named “xsi”.
7+
xsi_shapefile: Path to the xsi shapefile, that is then converted to a raster file. Be aware, that the attribute has to be named “mu”.
8+
oDefaults:
9+
default_mu: this is the default mu value, that gets assigned to all areas in the raster, that are not covered by shapefile-polygons
10+
default_xsi: this is the default xsi value, that gets assigned to all areas in the raster, that are not covered by shapefile-polygons
11+
12+
RunScript:
13+
oOnce everything is set up, run the script “runVariableVoellmyShapeToRaster.py”
14+
oIf libraries are missing use: pip install *name of missing library

‎avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.ini‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
[INPUT]
2+
# Input DEM raster file path
3+
# Example: /path/to/dem.asc
4+
dem =
5+
6+
# Input shapefile for mu
7+
# Example: /path/to/mu_shapefile.shp
8+
# be aware, that the attribute name has to be "mu"
9+
mu_shapefile =
10+
11+
# Input shapefile for xsi
12+
# Example: /path/to/xsi_shapefile.shp
13+
# be aware, that the attribute name has to be "xsi"
14+
xsi_shapefile =
15+
16+
[DEFAULTS]
17+
# Default value for areas not covered by shapefiles
18+
default_mu = 0.2
19+
20+
# Default value for areas not covered by shapefiles
21+
default_xsi = 2500
22+
23+
[OUTPUT]
24+
# Important: For the variable Voellmy calculations in the com1DFA algorithm to work, it is mandatory, that the files are stored in:
25+
# avaframe\data\*yourAvalancheDir*\Inputs\RASTERS\
26+
# also, the file names need to be of format *_mu.asc and *_xi.asc and both are required and same extent as DEM
27+
# Output raster mu
28+
# Example: /path/to/output_mu.asc
29+
mu_raster =
30+
31+
# Output raster xsi
32+
# Example: /path/to/output_xi.asc
33+
xsi_raster =

‎avaframe/com6RockAvalanche/variableVoellmyShapeToRaster.py‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
# -*- coding: utf-8 -*-
2+
"""
3+
Created on Tue Jan 21 11:43:18 2025
4+
5+
@author: Domi
6+
"""
7+
8+
import rasterio
9+
import numpy as np
10+
import configparser
11+
from rasterio.features import rasterize
12+
from shapely.geometry import shape, mapping
13+
from in2Trans.shpConversion import SHP2Array
14+
import logging
15+
16+
# Configure logging
17+
logging.basicConfig(level=logging.DEBUG)
18+
log = logging.getLogger(__name__)
19+
20+
def generateMuXsiRasters(configPath):
21+
"""
22+
Generate raster files for μ and ξ based on input DEM and shapefiles.
23+
24+
Parameters
25+
----------
26+
configPath : str
27+
Path to the configuration file.
28+
29+
Returns
30+
-------
31+
None
32+
"""
33+
# Load configuration
34+
config = configparser.ConfigParser()
35+
config.read(configPath)
36+
37+
demPath = config['INPUT']['dem']
38+
muShapefile = config['INPUT']['mu_shapefile']
39+
xsiShapefile = config['INPUT']['xsi_shapefile']
40+
defaultMu = float(config['DEFAULTS']['default_mu'])
41+
defaultXsi = float(config['DEFAULTS']['default_xsi'])
42+
muOutputPath = config['OUTPUT']['mu_raster']
43+
xsiOutputPath = config['OUTPUT']['xsi_raster']
44+
45+
# Read DEM
46+
with rasterio.open(demPath) as demSrc:
47+
demData = demSrc.read(1)
48+
demTransform = demSrc.transform
49+
demCrs = demSrc.crs
50+
demShape = demData.shape
51+
52+
# Helper function to rasterize shapefiles
53+
def rasterizeShapefile(shapefilePath, defaultValue, attributeName):
54+
if not shapefilePath:
55+
return np.full(demShape, defaultValue, dtype=np.float32)
56+
57+
shpData = SHP2Array(shapefilePath)
58+
shapes = []
59+
for i in range(shpData['nFeatures']):
60+
start = int(shpData['Start'][i])
61+
length = int(shpData['Length'][i])
62+
coords = [(shpData['x'][j], shpData['y'][j]) for j in range(start, start + length)]
63+
poly = shape({'type': 'Polygon', 'coordinates': [coords]})
64+
value = shpData['attributes'][i][attributeName] # Extract the attribute value
65+
shapes.append((mapping(poly), value))
66+
67+
rasterData = rasterize(
68+
shapes,
69+
out_shape=demShape,
70+
transform=demTransform,
71+
fill=defaultValue,
72+
all_touched=True,
73+
dtype=np.float32
74+
)
75+
return rasterData
76+
77+
# Generate μ and ξ rasters
78+
log.info("Rasterizing μ shapefile.")
79+
muRaster = rasterizeShapefile(muShapefile, defaultMu, "mu")
80+
81+
log.info("Rasterizing ξ shapefile.")
82+
xsiRaster = rasterizeShapefile(xsiShapefile, defaultXsi, "xsi")
83+
84+
# Save output rasters
85+
def saveRaster(outputPath, data):
86+
with rasterio.open(
87+
outputPath,
88+
'w',
89+
driver='GTiff',
90+
height=data.shape[0],
91+
width=data.shape[1],
92+
count=1,
93+
dtype=data.dtype,
94+
crs=demCrs,
95+
transform=demTransform,
96+
) as dst:
97+
dst.write(data, 1)
98+
99+
log.info("Saving μ raster to %s", muOutputPath)
100+
saveRaster(muOutputPath, muRaster)
101+
102+
log.info("Saving ξ raster to %s", xsiOutputPath)
103+
saveRaster(xsiOutputPath, xsiRaster)
104+
105+
log.info("Raster generation completed.")

‎avaframe/in2Trans/shpConversion.py‎

Lines changed: 12 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -58,7 +58,7 @@ def SHP2Array(infile, defname=None):
5858
number of features per line (parts)
5959
6060
"""
61-
# Input shapefile
61+
# Input shapefile
6262
sf = shapefile.Reader(str(infile))
6363
infile = pathlib.Path(infile)
6464
# set defaults for variables
@@ -94,6 +94,9 @@ def SHP2Array(infile, defname=None):
9494
start = 0
9595
nParts = []
9696

97+
# New: Create an empty list to store attributes
98+
attributes = []
99+
97100
for n, (item, rec) in enumerate(zip(shps, sf.records())):
98101
pts = item.points
99102
# if feature has no points - ignore
@@ -112,18 +115,20 @@ def SHP2Array(infile, defname=None):
112115

113116
# check if records are available and extract
114117
if records:
115-
# loop through fields
118+
# Extract attributes for the feature
119+
attr_dict = {}
116120
for (name, typ, size, deci), value in zip(sf.fields[1:], records[n].record):
117-
# get entity name
118121
name = name.lower()
122+
attr_dict[name] = value # Store attributes in dictionary
123+
124+
# Specific field handling (existing code)
119125
if name == "name":
120126
layername = str(value)
121127
if (name == "thickness") or (name == "d0"):
122128
thickness = value
123129
if name == "ci95":
124130
ci95 = value
125131
if name == "slope":
126-
# for dams
127132
slope = value
128133
if name == "rho":
129134
rho = value
@@ -133,11 +138,8 @@ def SHP2Array(infile, defname=None):
133138
iso = value
134139
if name == "layer":
135140
layerN = value
136-
# if name is still empty go through file again and take Layer instead
137-
if (type(layername) is bytes) or (layername is None):
138-
for (name, typ, size, deci), value in zip(sf.fields[1:], records[n].record):
139-
if name == "Layer":
140-
layername = value
141+
142+
attributes.append(attr_dict) # Add the attribute dictionary to the list
141143

142144
# if layer still not defined, use generic
143145
if layername is None:
@@ -173,6 +175,7 @@ def SHP2Array(infile, defname=None):
173175
SHPdata["layerName"] = layerNameList
174176
SHPdata["nParts"] = nParts
175177
SHPdata["nFeatures"] = len(Start)
178+
SHPdata["attributes"] = attributes # Add attributes to SHPdata
176179

177180
sf.close()
178181

‎avaframe/runVariableVoellmyShapeToRaster.py‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
# -*- coding: utf-8 -*-
2+
"""
3+
Created on Tue Jan 21 11:49:44 2025
4+
5+
@author: Domi
6+
"""
7+
8+
import argparse
9+
import pathlib
10+
from com6RockAvalanche.variableVoellmyShapeToRaster import generateMuXsiRasters
11+
from avaframe.in3Utils import cfgUtils
12+
from avaframe.in3Utils import logUtils
13+
from avaframe.in3Utils import initializeProject as initProj
14+
15+
def runMuXsiWorkflow(configPath=''):
16+
"""
17+
Run the workflow to generate \u03bc and \u03be rasters.
18+
19+
Parameters
20+
----------
21+
configPath : str
22+
Path to the configuration file.
23+
24+
Returns
25+
-------
26+
None
27+
"""
28+
logName = 'runMuXsi'
29+
30+
# Load general configuration file
31+
cfgMain = cfgUtils.getGeneralConfig()
32+
33+
# Load configuration file path from general config if not provided
34+
if configPath:
35+
cfgMain['MAIN']['configFile'] = configPath
36+
else:
37+
configPath = cfgMain['MAIN']['configFile']
38+
39+
configPath = pathlib.Path(configPath)
40+
41+
# Start logging
42+
log = logUtils.initiateLogger(configPath.parent, logName)
43+
log.info('MAIN SCRIPT')
44+
log.info('Using configuration file: %s', configPath)
45+
46+
# Clean input directory(ies) of old work files
47+
initProj.cleanSingleAvaDir(configPath.parent, deleteOutput=False)
48+
49+
# Run the raster generation process
50+
generateMuXsiRasters(str(configPath))
51+
52+
log.info('Workflow completed successfully.')
53+
54+
if __name__ == '__main__':
55+
parser = argparse.ArgumentParser(description='Run \u03bc and \u03be raster generation workflow')
56+
parser.add_argument('configPath', metavar='c', type=str, nargs='?', default='',
57+
help='Path to the configuration file')
58+
59+
args = parser.parse_args()
60+
runMuXsiWorkflow(str(args.configPath))

0 commit comments

Comments
 (0)