[FEATURE][processing] native drape features to z/m algorithms

Sets vertex z/m values to values sampled from a raster band.

Values can optionally be scaled using a (data definable) scale value.

python/plugins/processing/tests/testdata/custom/dem_lines.cpg

@@ -0,0 +1 @@
+UTF-8

python/plugins/processing/tests/testdata/custom/dem_lines.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

python/plugins/processing/tests/testdata/custom/dem_lines.qpj

@@ -0,0 +1 @@
+PROJCS["WGS 84 / Pseudo-Mercator",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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

python/plugins/processing/tests/testdata/expected/drape_lines.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

python/plugins/processing/tests/testdata/expected/drape_lines.qpj

@@ -0,0 +1 @@
+PROJCS["WGS 84 / Pseudo-Mercator",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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

python/plugins/processing/tests/testdata/expected/drape_lines_m.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

python/plugins/processing/tests/testdata/expected/drape_lines_m.qpj

@@ -0,0 +1 @@
+PROJCS["WGS 84 / Pseudo-Mercator",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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

python/plugins/processing/tests/testdata/expected/drape_points.prj

@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

python/plugins/processing/tests/testdata/expected/drape_points.qpj

@@ -0,0 +1 @@
+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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

python/plugins/processing/tests/testdata/expected/drape_points_m.prj

@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

python/plugins/processing/tests/testdata/expected/drape_points_m.qpj

@@ -0,0 +1 @@
+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.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml

@@ -5924,4 +5924,74 @@ tests:
         name: expected/create_parallel_lines.gml
         type: vector
 
+  - algorithm: native:drapetoz
+    name: Drape points to z
+    params:
+      BAND: 1
+      INPUT:
+        name: custom/sampling_points.gml
+        type: vector
+      NODATA: 0.0
+      RASTER:
+        name: dem.tif
+        type: raster
+      SCALE: 1.2
+    results:
+      OUTPUT:
+        name: expected/drape_points.shp
+        type: vector
+
+  - algorithm: native:drapetom
+    name: Drape points to m
+    params:
+      BAND: 1
+      INPUT:
+        name: custom/sampling_points.gml
+        type: vector
+      NODATA: 0.0
+      RASTER:
+        name: dem.tif
+        type: raster
+      SCALE: 1.2
+    results:
+      OUTPUT:
+        name: expected/drape_points_m.shp
+        type: vector
+
+  - algorithm: native:drapetoz
+    name: Drape lines to z
+    params:
+      BAND: 1
+      INPUT:
+        name: custom/dem_lines.shp
+        type: vector
+      NODATA: -9999.0
+      RASTER:
+        name: dem.tif
+        type: raster
+      SCALE: 1.0
+    results:
+      OUTPUT:
+        name: expected/drape_lines.shp
+        type: vector
+
+  - algorithm: native:drapetom
+    name: Drape lines to m
+    params:
+      BAND: 1
+      INPUT:
+        name: custom/dem_lines.shp
+        type: vector
+      NODATA: -9999.0
+      RASTER:
+        name: dem.tif
+        type: raster
+      SCALE: 1.0
+    results:
+      OUTPUT:
+        name: expected/drape_lines_m.shp
+        type: vector
+
+
+
 # See ../README.md for a description of the file format

src/analysis/CMakeLists.txt

@@ -31,6 +31,7 @@ SET(QGIS_ANALYSIS_SRCS
   processing/qgsalgorithmdbscanclustering.cpp
   processing/qgsalgorithmdifference.cpp
   processing/qgsalgorithmdissolve.cpp
+  processing/qgsalgorithmdrape.cpp
   processing/qgsalgorithmdropgeometry.cpp
   processing/qgsalgorithmdropmzvalues.cpp
   processing/qgsalgorithmexplode.cpp

src/analysis/processing/qgsalgorithmdrape.cpp

@@ -0,0 +1,252 @@
+/***************************************************************************
+                         qgsalgorithmdrape.cpp
+                         ---------------------
+    begin                : November 2017
+    copyright            : (C) 2017 by Nyall Dawson
+    email                : nyall dot dawson at gmail dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   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.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "qgsalgorithmdrape.h"
+
+///@cond PRIVATE
+
+
+QString QgsDrapeAlgorithmBase::group() const
+{
+  return QObject::tr( "Vector geometry" );
+}
+
+QString QgsDrapeAlgorithmBase::groupId() const
+{
+  return QStringLiteral( "vectorgeometry" );
+}
+
+QString QgsDrapeAlgorithmBase::outputName() const
+{
+  return QObject::tr( "Draped" );
+}
+
+void QgsDrapeAlgorithmBase::initParameters( const QVariantMap & )
+{
+  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "RASTER" ),
+                QObject::tr( "Raster layer" ) ) );
+  addParameter( new QgsProcessingParameterBand( QStringLiteral( "BAND" ),
+                QObject::tr( "Band number" ), 1, QStringLiteral( "RASTER" ) ) );
+
+  // nodata value
+  std::unique_ptr< QgsProcessingParameterNumber > nodata = qgis::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "NODATA" ),
+      QObject::tr( "Value for nodata or non-intersecting vertices" ), QgsProcessingParameterNumber::Double,
+      0.0 );
+  nodata->setIsDynamic( true );
+  nodata->setDynamicPropertyDefinition( QgsPropertyDefinition( QStringLiteral( "NODATA" ), QObject::tr( "Value for nodata or non-intersecting vertices" ), QgsPropertyDefinition::Double ) );
+  nodata->setDynamicLayerParameterName( QStringLiteral( "INPUT" ) );
+  addParameter( nodata.release() );
+
+  auto scaleParam = qgis::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "SCALE" ), QObject::tr( "Scale factor" ), QgsProcessingParameterNumber::Double, 1.0, false, 0.0 );
+  scaleParam->setIsDynamic( true );
+  scaleParam->setDynamicPropertyDefinition( QgsPropertyDefinition( QStringLiteral( "SCALE" ), QObject::tr( "Scale factor" ), QgsPropertyDefinition::Double ) );
+  scaleParam->setDynamicLayerParameterName( QStringLiteral( "INPUT" ) );
+  addParameter( scaleParam.release() );
+}
+
+bool QgsDrapeAlgorithmBase::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
+{
+  mNoData = parameterAsDouble( parameters, QStringLiteral( "NODATA" ), context );
+  mDynamicNoData = QgsProcessingParameters::isDynamic( parameters, QStringLiteral( "NODATA" ) );
+  if ( mDynamicNoData )
+    mNoDataProperty = parameters.value( QStringLiteral( "NODATA" ) ).value< QgsProperty >();
+
+  mScale = parameterAsDouble( parameters, QStringLiteral( "SCALE" ), context );
+  mDynamicScale = QgsProcessingParameters::isDynamic( parameters, QStringLiteral( "SCALE" ) );
+  if ( mDynamicScale )
+    mScaleProperty = parameters.value( QStringLiteral( "SCALE" ) ).value< QgsProperty >();
+
+  QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "RASTER" ), context );
+
+  if ( !layer )
+    throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "RASTER" ) ) );
+
+  mBand = parameterAsInt( parameters, QStringLiteral( "BAND" ), context );
+  if ( mBand < 1 || mBand > layer->bandCount() )
+    throw QgsProcessingException( QObject::tr( "Invalid band number for BAND (%1): Valid values for input raster are 1 to %2" ).arg( mBand )
+                                  .arg( layer->bandCount() ) );
+
+  std::unique_ptr< QgsRasterInterface > provider( layer->dataProvider()->clone() );
+  QgsRasterDataProvider *dp = dynamic_cast< QgsRasterDataProvider * >( provider.get() );
+  if ( !dp )
+    throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "RASTER" ) ) );
+
+  mRasterProvider.reset( dp );
+  provider.release();
+
+  return true;
+}
+
+QgsFeatureList QgsDrapeAlgorithmBase::processFeature( const QgsFeature &feature, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
+  if ( !mCreatedTransform )
+  {
+    mCreatedTransform = true;
+    mTransform = QgsCoordinateTransform( sourceCrs(), mRasterProvider->crs(), context.transformContext() );
+  }
+
+  QgsFeature f = feature;
+  if ( f.hasGeometry() )
+  {
+    QgsGeometry geometry = f.geometry();
+
+    double nodata = mNoData;
+    if ( mDynamicNoData )
+      nodata = mNoDataProperty.valueAsDouble( context.expressionContext(), nodata );
+
+    double scale = mScale;
+    if ( mDynamicScale )
+      scale = mScaleProperty.valueAsDouble( context.expressionContext(), scale );
+
+    prepareGeometry( geometry, nodata );
+
+    geometry.transformVertices( [ = ]( const QgsPoint & p )->QgsPoint
+    {
+      QgsPointXY t;
+      double val = nodata;
+      try
+      {
+        t = mTransform.transform( p );
+        bool ok = false;
+        val = mRasterProvider->sample( t, mBand, &ok );
+        if ( !ok )
+          val = nodata;
+        else
+          val *= scale;
+      }
+      catch ( QgsCsException & )
+      {
+        feedback->reportError( QObject::tr( "Transform error while reprojecting feature {}" ).arg( f.id() ) );
+      }
+
+      return drapeVertex( p, val );
+    } );
+
+    f.setGeometry( geometry );
+  }
+  return QgsFeatureList() << f;
+}
+
+
+//
+// QgsDrapeToZAlgorithm
+//
+
+QString QgsDrapeToZAlgorithm::name() const
+{
+  return QStringLiteral( "drapetoz" );
+}
+
+QString QgsDrapeToZAlgorithm::displayName() const
+{
+  return QObject::tr( "Drape (Sample raster to z)" );
+}
+
+QStringList QgsDrapeToZAlgorithm::tags() const
+{
+  return QObject::tr( "3d,vertex,vertices,elevation,sample" ).split( ',' );
+}
+
+QString QgsDrapeToZAlgorithm::shortHelpString() const
+{
+  return QObject::tr( "This algorithm sets the z value of every vertex in the feature geometry to a value sampled from a band within a raster layer." )
+         + QStringLiteral( "\n\n" )
+         + QObject::tr( "The raster values can optionally be scaled by a preset amount." );
+}
+
+QString QgsDrapeToZAlgorithm::shortDescription() const
+{
+  return QObject::tr( "Sets the z value for vertices to values sampled from a raster layer." );
+}
+
+QgsDrapeToZAlgorithm *QgsDrapeToZAlgorithm::createInstance() const
+{
+  return new QgsDrapeToZAlgorithm();
+}
+
+QgsWkbTypes::Type QgsDrapeToZAlgorithm::outputWkbType( QgsWkbTypes::Type inputWkbType ) const
+{
+  QgsWkbTypes::Type wkb = inputWkbType;
+  return QgsWkbTypes::addZ( wkb );
+}
+
+void QgsDrapeToZAlgorithm::prepareGeometry( QgsGeometry &geometry, double defaultVal ) const
+{
+  geometry.get()->addZValue( defaultVal );
+}
+
+QgsPoint QgsDrapeToZAlgorithm::drapeVertex( const QgsPoint &p, double rasterVal ) const
+{
+  return QgsPoint( p.wkbType(), p.x(), p.y(), rasterVal, p.m() );
+}
+
+//
+// QgsDrapeToMAlgorithm
+//
+
+QString QgsDrapeToMAlgorithm::name() const
+{
+  return QStringLiteral( "drapetom" );
+}
+
+QString QgsDrapeToMAlgorithm::displayName() const
+{
+  return QObject::tr( "Sample raster to m-values" );
+}
+
+QStringList QgsDrapeToMAlgorithm::tags() const
+{
+  return QObject::tr( "vertex,vertices,sample,measure" ).split( ',' );
+}
+
+QString QgsDrapeToMAlgorithm::shortHelpString() const
+{
+  return QObject::tr( "This algorithm sets the m-value for every vertex in the feature geometry to a value sampled from a band within a raster layer." )
+         + QStringLiteral( "\n\n" )
+         + QObject::tr( "The raster values can optionally be scaled by a preset amount." );
+}
+
+QString QgsDrapeToMAlgorithm::shortDescription() const
+{
+  return QObject::tr( "Sets the m-value for vertices to values sampled from a raster layer." );
+}
+
+QgsDrapeToMAlgorithm *QgsDrapeToMAlgorithm::createInstance() const
+{
+  return new QgsDrapeToMAlgorithm();
+}
+
+QgsWkbTypes::Type QgsDrapeToMAlgorithm::outputWkbType( QgsWkbTypes::Type inputWkbType ) const
+{
+  QgsWkbTypes::Type wkb = inputWkbType;
+  return QgsWkbTypes::addM( wkb );
+}
+
+void QgsDrapeToMAlgorithm::prepareGeometry( QgsGeometry &geometry, double defaultVal ) const
+{
+  geometry.get()->addMValue( defaultVal );
+}
+
+QgsPoint QgsDrapeToMAlgorithm::drapeVertex( const QgsPoint &p, double rasterVal ) const
+{
+  return QgsPoint( p.wkbType(), p.x(), p.y(), p.z(), rasterVal );
+}
+
+
+///@endcond
+
+