Skip to content

Commit

Permalink
Improve raster layers self-saving, improve DEUG output, rm old files,…
Browse files Browse the repository at this point in the history 
… modular parameters
  • Loading branch information
@emanuelegissi
emanuelegissi committed Aug 24, 2023
1 parent 18d6f0a commit 3e73494
Show file tree
Hide file tree
Showing 10 changed files with 2,481 additions and 273 deletions.
201 changes: 72 additions & 129 deletions qgis2fds_export_algo.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,6 @@
QgsProcessingParameterRasterDestination,
QgsProcessingParameterVectorDestination,
QgsCoordinateReferenceSystem,
QgsCoordinateTransform,
QgsPointXY,
QgsPoint,
QgsField,
NULL,
edit,
Expand All @@ -27,7 +24,7 @@
from qgis.PyQt.QtCore import QVariant

import processing
import math, time
import time
from .qgis2fds_params import *
from .types import (
utils,
Expand All @@ -39,6 +36,7 @@
Texture,
Wind,
)
from . import utils2

DEBUG = True

Expand Down Expand Up @@ -94,13 +92,6 @@ def initAlgorithm(self, config=None):
optional=True,
)
)
self.addParameter(
QgsProcessingParameterVectorDestination(
"ExtentDebug",
"Extent debug",
optional=True,
)
)

def processAlgorithm(self, parameters, context, model_feedback):
feedback = QgsProcessingMultiStepFeedback(6, model_feedback)
Expand All @@ -115,54 +106,45 @@ def processAlgorithm(self, parameters, context, model_feedback):
"feedback": feedback,
"project": project,
}

# FIXME move to required place

chid = ChidParam.get(**kwargs)
fds_path = FDSPathParam.get(**kwargs)
extent_layer = ExtentLayerParam.get(**kwargs) # in project crs
pixel_size = PixelSizeParam.get(**kwargs)
origin = OriginParam.get(**kwargs) # in project crs
dem_layer = DEMLayerParam.get(**kwargs)
landuse_layer = LanduseLayerParam.get(**kwargs)
landuse_type_filepath = LanduseTypeFilepathParam.get(**kwargs)

fire_layer = FireLayer.get(**kwargs)
tex_layer = TexLayerParam.get(**kwargs)
tex_pixel_size = TexPixelSizeParam.get(**kwargs)
nmesh = NMeshParam.get(**kwargs)
cell_size = CellSizeParam.get(**kwargs)

export_obst = ExportOBSTParam.get(**kwargs)
t_begin = StartTimeParam.get(**kwargs)
t_end = EndTimeParam.get(**kwargs)
wind_filepath = WindFilepathParam.get(**kwargs)
text_filepath = TextFilepathParam.get(**kwargs)

# Check other params for None values
# origin is checked later

if landuse_layer:
if not landuse_type_filepath:
raise QgsProcessingException(
f"Specify the landuse type for the {landuse_layer.name()} landuse layer.\n"
)
if not cell_size:
cell_size = pixel_size
# Get extent_layer and origin
# Calc wgs84_extent and wgs84_origin

# Get wgs84_extent from extent_layer in WGS84
extent_layer = ExtentLayerParam.get(**kwargs) # in project crs
origin = OriginParam.get(**kwargs) # in project crs

wgs84_crs = QgsCoordinateReferenceSystem("EPSG:4326")
wgs84_extent = _transform_extent(
wgs84_extent = utils2.transform_extent(
extent=extent_layer.extent(),
source_crs=extent_layer.crs(),
dest_crs=wgs84_crs,
)
feedback.setProgressText(f"\nWGS84 extent: {wgs84_extent}")

# Check origin and get wgs84_origin

if origin:
wgs84_origin = _transform_point(
point=origin, source_crs=project.crs(), dest_crs=wgs84_crs
wgs84_origin = utils2.transform_point(
point=origin,
source_crs=project.crs(),
dest_crs=wgs84_crs,
)
else:
wgs84_origin = wgs84_extent.center() # in project crs
wgs84_origin = wgs84_extent.center()

feedback.setProgressText(f"WGS84 origin: {wgs84_origin}")

# Calc applicable UTM crs at the origin
Expand All @@ -173,42 +155,40 @@ def processAlgorithm(self, parameters, context, model_feedback):

# Calc utm_origin

utm_origin = _transform_point(
utm_origin = utils2.transform_point(
point=wgs84_origin,
source_crs=wgs84_crs,
dest_crs=utm_crs,
)
feedback.setProgressText(f"UTM origin: {utm_origin}")

# Calc utm_extent
# Align it to pixel_size, better for sampling grid and DEM interpolation
# Get pixel_size
# Calc and adapt utm_extent:
# align it to the pixel_size, better for sampling grid and DEM interpolation

pixel_size = PixelSizeParam.get(**kwargs)

utm_extent = _transform_extent(
utm_extent = utils2.transform_extent(
extent=wgs84_extent,
source_crs=wgs84_crs,
dest_crs=utm_crs,
)
e = 1e-6 # epsilon used to nudge the native:creategrid algo
w = math.ceil(utm_extent.width() / pixel_size) * pixel_size + e
h = math.ceil(utm_extent.height() / pixel_size) * pixel_size + e
x_min, x_max, y_min, y_max = (
utm_extent.xMinimum(),
utm_extent.xMaximum(),
utm_extent.yMinimum(),
utm_extent.yMaximum(),
)
utm_extent.setXMaximum(x_min + w)
utm_extent.setYMinimum(y_max - h)
feedback.setProgressText(f"UTM extent: {utm_extent}")
feedback.setProgressText(f"Extent size: {x_max-x_min:.1f}x{y_max-y_min:.1f}m")

utm_extent = utils2.get_extent_multiple_of_pixels(
extent=utm_extent, pixel_size=pixel_size, epsilon=1e-6
) # epsilon used to nudge the native:creategrid algo

msg = f"UTM extent size: {utm_extent.width():.1f}x{utm_extent.height():.1f}m"
feedback.setProgressText(msg)

if DEBUG:
_run_extent_to_layer(
parameters=parameters,
utils2.show_extent(
context=context,
feedback=feedback,
extent=utm_extent,
extent_crs=utm_crs,
name="DEBUG UTM extent",
style="Extent layer style.qml",
)

feedback.setCurrentStep(1)
Expand All @@ -227,7 +207,7 @@ def processAlgorithm(self, parameters, context, model_feedback):
"VSPACING": pixel_size,
"HOVERLAY": 0.0,
"VOVERLAY": 0.0,
"OUTPUT": DEBUG and parameters["UtmGrid"] or QgsProcessing.TEMPORARY_OUTPUT,
"OUTPUT": QgsProcessing.TEMPORARY_OUTPUT,
}
outputs["UtmGrid"] = processing.run(
"native:creategrid",
Expand Down Expand Up @@ -260,31 +240,19 @@ def processAlgorithm(self, parameters, context, model_feedback):
utm_idem_extent = utm_extent.buffered(pixel_size / 2.0 - e)

if DEBUG:
_run_extent_to_layer(
parameters=parameters,
utils2.show_extent(
context=context,
feedback=feedback,
extent=utm_idem_extent,
extent_crs=utm_crs,
name="DEBUG UTM DEM extent",
style="Extent layer style.qml",
)

# Check DEM layer is local, otherwise save it

dem_layer = DEMLayerParam.get_local(
algo=self,
parameters=parameters,
context=context,
feedback=feedback,
project=project,
extent=_transform_extent(
extent=utm_idem_extent, source_crs=utm_crs, dest_crs=dem_layer.crs()
), # FIXME not useful here
)

return {} # FIXME

# Clip, reproject, and interpolate the DEM

dem_layer = DEMLayerParam.get(**kwargs, extent=utm_extent, extent_crs=utm_crs)

feedback.setProgressText("\nClip, reproject, and interpolate the DEM...")
t0 = time.time()
alg_params = {
Expand Down Expand Up @@ -330,7 +298,7 @@ def processAlgorithm(self, parameters, context, model_feedback):
"OFFSET": 0,
"RASTER": outputs["Interpolation"]["OUTPUT"],
"SCALE": 1,
"OUTPUT": DEBUG and parameters["UtmGrid"] or QgsProcessing.TEMPORARY_OUTPUT,
"OUTPUT": QgsProcessing.TEMPORARY_OUTPUT,
}
outputs["SetZFromDem"] = processing.run(
"native:setzfromraster",
Expand All @@ -348,14 +316,27 @@ def processAlgorithm(self, parameters, context, model_feedback):
# Sample landuse values from landuse layer
# The new data column name is bc1

if landuse_layer and landuse_type_filepath:
landuse_layer = LanduseLayerParam.get(
**kwargs, extent=utm_extent, extent_crs=utm_crs
)
landuse_type_filepath = LanduseTypeFilepathParam.get(**kwargs)

if landuse_layer:
feedback.setProgressText("\nSample landuse values from landuse layer...")
t0 = time.time()

# Check landuse type exists

if not landuse_type_filepath:
msg = f"Landuse type not available, cannot proceed."
raise QgsProcessingException(msg)

# Sample
alg_params = {
"COLUMN_PREFIX": "bc", # creates the bc1 field from landuse
"INPUT": outputs["SetZFromDem"]["OUTPUT"],
"RASTERCOPY": landuse_layer,
"OUTPUT": parameters["UtmGrid"],
"OUTPUT": QgsProcessing.TEMPORARY_OUTPUT,
}
outputs["SampleRasterValues"] = processing.run(
"native:rastersampling",
Expand All @@ -364,13 +345,12 @@ def processAlgorithm(self, parameters, context, model_feedback):
feedback=feedback,
is_child_algorithm=True,
)
# results["UtmGrid"] = outputs["SampleRasterValues"]["OUTPUT"]
sampling_layer = context.getMapLayer(
outputs["SampleRasterValues"]["OUTPUT"]
)
feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
else:
feedback.setProgressText("\nNo landuse layer or type.")
feedback.setProgressText("\nNo landuse layer.")
sampling_layer = context.getMapLayer(outputs["SetZFromDem"]["OUTPUT"])
with edit(sampling_layer): # create an empty bc1 field
attributes = list((QgsField("bc1", QVariant.Int),))
Expand All @@ -381,6 +361,14 @@ def processAlgorithm(self, parameters, context, model_feedback):
if feedback.isCanceled():
return {}

if DEBUG:
utils2.show_layer(
context=context,
feedback=feedback,
layer=sampling_layer,
name="DEBUG Sampling layer",
)

# Get landuse type
# (we need landuse_type.bc_out_default and .bc_in_default)

Expand Down Expand Up @@ -496,6 +484,12 @@ def processAlgorithm(self, parameters, context, model_feedback):
if feedback.isCanceled():
return {}

# Get cell_size

cell_size = CellSizeParam.get(**kwargs)
if not cell_size:
cell_size = pixel_size

domain = Domain(
feedback=feedback,
utm_crs=utm_crs,
Expand Down Expand Up @@ -547,57 +541,6 @@ def createInstance(self):
# FIXME move elsewhere


def _transform_extent(extent, source_crs, dest_crs):
_tr = QgsCoordinateTransform(source_crs, dest_crs, QgsProject.instance())
return _tr.transformBoundingBox(extent)


def _transform_point(point, source_crs, dest_crs):
_tr = QgsCoordinateTransform(source_crs, dest_crs, QgsProject.instance())
return _tr.transform(QgsPointXY(point))


def _run_create_spatial_index(parameters, context, feedback, vector_layer):
"""Create spatial index of vector layer to speed up the next process."""
feedback.setProgressText("\nCreate spatial index...")
t0 = time.time()
alg_params = {
"INPUT": vector_layer,
}
output = processing.run(
"native:createspatialindex",
alg_params,
context=context,
feedback=feedback,
is_child_algorithm=True,
)
feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
return output


def _run_extent_to_layer(parameters, context, feedback, extent, extent_crs):
"""Show extent as vector layer."""
x_min, x_max, y_min, y_max = (
extent.xMinimum(),
extent.xMaximum(),
extent.yMinimum(),
extent.yMaximum(),
)
extent_str = f"{x_min}, {x_max}, {y_min}, {y_max} [{extent_crs.authid()}]"
feedback.setProgressText(f"Extent to layer: {extent_str} ...")
alg_params = {
"INPUT": extent_str,
"OUTPUT": parameters["ExtentDebug"], # FIXME how to make it unlinked?
}
return processing.run(
"native:extenttolayer",
alg_params,
context=context,
feedback=feedback,
is_child_algorithm=True,
)


def _load_fire_layer_bc(
parameters, context, feedback, sampling_layer, fire_layer, bc_field, bc_default
):
Expand Down
Loading

0 comments on commit 3e73494

Please sign in to comment.