Fixed a problem with CERN CRS, moved the DEM clipping to vector layer…

…, added timing

qgis2fds_export_algo.py

@@ -21,7 +21,7 @@
  QgsCoordinateTransform,
  QgsPoint,
 )
-import processing, math
+import processing, math, time
 from .qgis2fds_params import *
 from .types import (
  utils,
@@ -78,13 +78,6 @@ def initAlgorithm(self, config=None):
  )
 
  if DEBUG:
- self.addParameter(
- QgsProcessingParameterRasterDestination(
- "ClippedDemLayer",
- "Clipped DEM layer",
- optional=True,
- )
- )
  self.addParameter(
  QgsProcessingParameterVectorDestination(
  "UtmDemPoints",
@@ -111,6 +104,7 @@ def initAlgorithm(self, config=None):
  def processAlgorithm(self, parameters, context, model_feedback):
  feedback = QgsProcessingMultiStepFeedback(9, model_feedback) # FIXME
  results, outputs, project = {}, {}, QgsProject.instance()
+ time0 = time.time()
 
  # Load parameter values
 
@@ -238,6 +232,7 @@ def processAlgorithm(self, parameters, context, model_feedback):
 
  # Create UTM sampling grid
  feedback.setProgressText("\nCreate UTM sampling grid...")
+ t0 = time.time()
  spacing = pixel_size
  alg_params = {
  "CRS": utm_crs,
@@ -256,6 +251,7 @@ def processAlgorithm(self, parameters, context, model_feedback):
  feedback=feedback,
  is_child_algorithm=True,
  )
+ feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
 
  # Check UTM sampling grid
  feedback.setProgressText("\nCheck UTM sampling grid...")
@@ -269,54 +265,56 @@ def processAlgorithm(self, parameters, context, model_feedback):
  if feedback.isCanceled():
  return {}
 
- # Clip DEM by needed extent
- feedback.setProgressText("\nClip DEM by needed extent...")
+ # Transform DEM pixels to points
+ feedback.setProgressText("\nTransform DEM pixels to points...")
+ t0 = time.time()
  alg_params = {
- "DATA_TYPE": 0, # use input layer data type
- "EXTRA": "",
- "INPUT": dem_layer, # in PROJWIN crs
- "NODATA": -999.0,
- "OPTIONS": "",
- "OVERCRS": True,
- "PROJWIN": clipped_dem_extent,
- "OUTPUT": DEBUG
- and parameters["ClippedDemLayer"]
- or QgsProcessing.TEMPORARY_OUTPUT,
+ "FIELD_NAME": "DEM",
+ "INPUT_RASTER": dem_layer,
+ "RASTER_BAND": 1,
+ "OUTPUT": QgsProcessing.TEMPORARY_OUTPUT,
  }
- outputs["ClippedDemLayer"] = processing.run(
- "gdal:cliprasterbyextent",
+ outputs["DemPixelsToPoints"] = processing.run(
+ "native:pixelstopoints",
  alg_params,
  context=context,
  feedback=feedback,
  is_child_algorithm=True,
  )
+ feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
 
- feedback.setCurrentStep(4)
+ feedback.setCurrentStep(2)
  if feedback.isCanceled():
  return {}
 
- # Transform clipped DEM pixels to points
- feedback.setProgressText("\nTransform clipped DEM pixels to points...")
+ # Clip DEM points by requested extent
+ # Spatial index does not improve the performance here
+ feedback.setProgressText("\nClip DEM points by requested extent...")
+ t0 = time.time()
  alg_params = {
- "FIELD_NAME": "DEM",
- "INPUT_RASTER": outputs["ClippedDemLayer"]["OUTPUT"],
- "RASTER_BAND": 1,
- "OUTPUT": QgsProcessing.TEMPORARY_OUTPUT,
+ "INPUT": outputs["DemPixelsToPoints"]["OUTPUT"],
+ "EXTENT": clipped_dem_extent,
+ "CLIP": False,
+ "OUTPUT": DEBUG
+ and parameters["UtmDemPoints"]
+ or QgsProcessing.TEMPORARY_OUTPUT,
  }
  outputs["DemPixelsToPoints"] = processing.run(
- "native:pixelstopoints",
+ "native:extractbyextent",
  alg_params,
  context=context,
  feedback=feedback,
  is_child_algorithm=True,
  )
+ feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
 
- feedback.setCurrentStep(5)
+ feedback.setCurrentStep(5) # FIXME
  if feedback.isCanceled():
  return {}
 
  # Reproject DEM points to UTM
  feedback.setProgressText("\nReproject DEM points to UTM...")
+ t0 = time.time()
  alg_params = {
  "CONVERT_CURVED_GEOMETRIES": False,
  "INPUT": outputs["DemPixelsToPoints"]["OUTPUT"],
@@ -333,16 +331,35 @@ def processAlgorithm(self, parameters, context, model_feedback):
  feedback=feedback,
  is_child_algorithm=True,
  )
+ feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
 
  feedback.setCurrentStep(6)
  if feedback.isCanceled():
  return {}
 
+ # # Create spatial index to speed up the next process
+ # feedback.setProgressText("\nCreate spatial index...")
+ # t0 = time.time()
+ # input = outputs["UtmDemPoints"]["OUTPUT"]
+ # alg_params = {
+ # "INPUT": input,
+ # }
+ # processing.run(
+ # "native:createspatialindex",
+ # alg_params,
+ # context=context,
+ # feedback=feedback,
+ # is_child_algorithm=True,
+ # )
+ # feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
+
  # Interpolate UTM DEM points to DEM raster layer (Grid, IDW with nearest neighbor searching)
  # aligned to UTM sampling grid
+ # Spatial index does not improve the performance here
  feedback.setProgressText(
  "\nInterpolate UTM DEM points to DEM raster layer (IDW)..."
  )
+ t0 = time.time()
  radius = max(dem_layer_rx, dem_layer_ry)
  e = utm_extent
  extra = f"-txe {e.xMinimum()} {e.xMaximum()} -tye {e.yMinimum()} {e.yMaximum()} -tr {pixel_size} {pixel_size}"
@@ -369,12 +386,14 @@ def processAlgorithm(self, parameters, context, model_feedback):
  feedback=feedback,
  is_child_algorithm=True,
  )
+ feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
 
  # # Interpolate UTM DEM points to DEM raster layer (TIN interpolation)
  # # aligned to UTM sampling grid
  # feedback.setProgressText(
  # "\nInterpolate UTM DEM points to DEM raster layer (TIN)..."
  # )
+ # t0 = time.time()
  # utm_dem_layer = outputs["UtmDemPoints"]["OUTPUT"]
  # layer_source = context.getMapLayer(utm_dem_layer).source()
  # interpolation_source = 0
@@ -397,6 +416,7 @@ def processAlgorithm(self, parameters, context, model_feedback):
  # feedback=feedback,
  # is_child_algorithm=True,
  # )
+ # feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
 
  feedback.setCurrentStep(7)
  if feedback.isCanceled():
@@ -406,6 +426,7 @@ def processAlgorithm(self, parameters, context, model_feedback):
  feedback.setProgressText(
  "\nSample Z values from interpolated DEM raster layer..."
  )
+ t0 = time.time()
  alg_params = {
  "BAND": 1,
  "INPUT": outputs["UtmGrid"]["OUTPUT"],
@@ -422,13 +443,15 @@ def processAlgorithm(self, parameters, context, model_feedback):
  feedback=feedback,
  is_child_algorithm=True,
  )
+ feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
 
  feedback.setCurrentStep(8)
  if feedback.isCanceled():
  return {}
 
  # Sample landuse values from landuse layer
  feedback.setProgressText("\nSample landuse values from landuse layer...")
+ t0 = time.time()
  if landuse_layer and landuse_type_filepath:
  alg_params = {
  "COLUMN_PREFIX": "landuse",
@@ -451,6 +474,7 @@ def processAlgorithm(self, parameters, context, model_feedback):
  )
  else:
  fds_grid_layer = context.getMapLayer(outputs["SetZFromDem"]["OUTPUT"])
+ feedback.setProgressText(f"time: {time.time()-t0:.1f}s")
 
  feedback.setCurrentStep(9)
  if feedback.isCanceled():