sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-10-04 00:04:03 -04:00
Code Issues Packages Projects Releases Wiki Activity

port Find projection algorithm to C++

Browse Source
This commit is contained in:
Alexander Bruy 2025-03-21 08:50:28 +00:00 committed by Nyall Dawson
parent 3e050475df
commit 43e8499071
8 changed files with 234 additions and 188 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
python/plugins/processing/algs/help/qgis.yaml
View File

@ -34,13 +34,6 @@ qgis:executesql: >
The result of the query will be added as a new layer.
qgis:findprojection: >
This algorithm allows creation of a shortlist of possible candidate coordinate reference systems for a layer with an unknown projection.
The expected area which the layer should reside in must be specified via the target area parameter.
The algorithm operates by testing the layer's extent in every known reference system and listing any in which the bounds would fall near the target area if the layer was in this projection.
qgis:generatepointspixelcentroidsalongline: >
This algorithm generates a point vector layer from an input raster and line layer.
The points correspond to the pixel centroids that intersect the line layer.

178
python/plugins/processing/algs/qgis/FindProjection.py
View File

@ -1,178 +0,0 @@
"""
***************************************************************************
FindProjection.py
-----------------
Date : February 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. *
* *
***************************************************************************
"""
__author__ = "Nyall Dawson"
__date__ = "February 2017"
__copyright__ = "(C) 2017, Nyall Dawson"
import os
from qgis.core import (
QgsGeometry,
QgsFeature,
QgsFeatureSink,
QgsField,
QgsFields,
QgsCoordinateReferenceSystem,
QgsCoordinateTransform,
QgsCoordinateTransformContext,
QgsWkbTypes,
QgsProcessingException,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterExtent,
QgsProcessingParameterCrs,
QgsProcessingParameterFeatureSink,
QgsProcessingParameterDefinition,
)
from qgis.PyQt.QtCore import QMetaType
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]
class FindProjection(QgisAlgorithm):
INPUT = "INPUT"
TARGET_AREA = "TARGET_AREA"
TARGET_AREA_CRS = "TARGET_AREA_CRS"
OUTPUT = "OUTPUT"
def tags(self):
return self.tr(
"crs,srs,coordinate,reference,system,guess,estimate,finder,determine"
).split(",")
def group(self):
return self.tr("Vector general")
def groupId(self):
return "vectorgeneral"
def __init__(self):
super().__init__()
def initAlgorithm(self, config=None):
self.addParameter(
QgsProcessingParameterFeatureSource(self.INPUT, self.tr("Input layer"))
)
extent_parameter = QgsProcessingParameterExtent(
self.TARGET_AREA, self.tr("Target area for layer")
)
self.addParameter(extent_parameter)
# deprecated
crs_param = QgsProcessingParameterCrs(
self.TARGET_AREA_CRS, "Target area CRS", optional=True
)
crs_param.setFlags(
crs_param.flags() | QgsProcessingParameterDefinition.Flag.FlagHidden
)
self.addParameter(crs_param)
self.addParameter(
QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr("CRS candidates"))
)
def name(self):
return "findprojection"
def displayName(self):
return self.tr("Find projection")
def processAlgorithm(self, parameters, context, feedback):
source = self.parameterAsSource(parameters, self.INPUT, context)
if source is None:
raise QgsProcessingException(
self.invalidSourceError(parameters, self.INPUT)
)
extent = self.parameterAsExtent(parameters, self.TARGET_AREA, context)
target_crs = self.parameterAsExtentCrs(parameters, self.TARGET_AREA, context)
if self.TARGET_AREA_CRS in parameters:
c = self.parameterAsCrs(parameters, self.TARGET_AREA_CRS, context)
if c.isValid():
target_crs = c
target_geom = QgsGeometry.fromRect(extent)
fields = QgsFields()
fields.append(QgsField("auth_id", QMetaType.Type.QString, "", 20))
(sink, dest_id) = self.parameterAsSink(
parameters,
self.OUTPUT,
context,
fields,
QgsWkbTypes.Type.NoGeometry,
QgsCoordinateReferenceSystem(),
)
if sink is None:
raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))
# make intersection tests nice and fast
engine = QgsGeometry.createGeometryEngine(target_geom.constGet())
engine.prepareGeometry()
layer_bounds = QgsGeometry.fromRect(source.sourceExtent())
crses_to_check = QgsCoordinateReferenceSystem.validSrsIds()
total = 100.0 / len(crses_to_check)
found_results = 0
transform_context = QgsCoordinateTransformContext()
for current, srs_id in enumerate(crses_to_check):
if feedback.isCanceled():
break
candidate_crs = QgsCoordinateReferenceSystem.fromSrsId(srs_id)
if not candidate_crs.isValid():
continue
transform_candidate = QgsCoordinateTransform(
candidate_crs, target_crs, transform_context
)
transform_candidate.setBallparkTransformsAreAppropriate(True)
transform_candidate.disableFallbackOperationHandler(True)
transformed_bounds = QgsGeometry(layer_bounds)
try:
if transformed_bounds.transform(transform_candidate) != 0:
continue
except:
continue
try:
if engine.intersects(transformed_bounds.constGet()):
feedback.pushInfo(
self.tr("Found candidate CRS: {}").format(
candidate_crs.authid()
)
)
f = QgsFeature(fields)
f.setAttributes([candidate_crs.authid()])
sink.addFeature(f, QgsFeatureSink.Flag.FastInsert)
found_results += 1
except:
continue
feedback.setProgress(int(current * total))
if found_results == 0:
feedback.reportError(self.tr("No matching projections found"))
sink.finalize()
return {self.OUTPUT: dest_id}

2
python/plugins/processing/algs/qgis/QgisAlgorithmProvider.py
View File

@ -30,7 +30,6 @@ from .BoxPlot import BoxPlot
from .EliminateSelection import EliminateSelection
from .ExecuteSQL import ExecuteSQL
from .FieldPyculator import FieldsPyculator
from .FindProjection import FindProjection
from .Heatmap import Heatmap
from .HubDistanceLines import HubDistanceLines
from .HubDistancePoints import HubDistancePoints
@ -90,7 +89,6 @@ class QgisAlgorithmProvider(QgsProcessingProvider):
EliminateSelection(),
ExecuteSQL(),
FieldsPyculator(),
FindProjection(),
Heatmap(),
HubDistanceLines(),
HubDistancePoints(),

2
python/plugins/processing/tests/testdata/qgis_algorithm_tests3.yaml vendored
View File

@ -396,7 +396,7 @@ tests:
# This is VERY slow on proj 6+ builds due to the excessive amount of message log noise it causes to be emitted.
# We need a way to disable the default handlers on individual QgsCoordinateTransform objects in order to avoid
# this and speed the algorithm back up
# - algorithm: qgis:findprojection
# - algorithm: native:findprojection
# name: Find projection
# params:
# INPUT:

1
src/analysis/CMakeLists.txt
View File

@ -151,6 +151,7 @@ set(QGIS_ANALYSIS_SRCS
processing/qgsalgorithmfilter.cpp
processing/qgsalgorithmfilterbygeometry.cpp
processing/qgsalgorithmfiltervertices.cpp
processing/qgsalgorithmfindprojection.cpp
processing/qgsalgorithmfixgeometries.cpp
processing/qgsalgorithmflattenrelationships.cpp
processing/qgsalgorithmforcerhr.cpp

180
src/analysis/processing/qgsalgorithmfindprojection.cpp Normal file
View File

@ -0,0 +1,180 @@
/***************************************************************************
qgsalgorithmfindprojection.cpp
---------------------
begin : March 2025
copyright : (C) 2025 by Alexander Bruy
email : alexander dot bruy 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 "qgsalgorithmfindprojection.h"
#include "qgsgeometryengine.h"
///@cond PRIVATE
QString QgsFindProjectionAlgorithm::name() const
{
return QStringLiteral( "findprojection" );
}
QString QgsFindProjectionAlgorithm::displayName() const
{
return QObject::tr( "Find projection" );
}
QStringList QgsFindProjectionAlgorithm::tags() const
{
return QObject::tr( "crs,srs,coordinate,reference,system,guess,estimate,finder,determine" ).split( ',' );
}
QString QgsFindProjectionAlgorithm::group() const
{
return QObject::tr( "Vector general" );
}
QString QgsFindProjectionAlgorithm::groupId() const
{
return QStringLiteral( "vectorgeneral" );
}
QString QgsFindProjectionAlgorithm::shortHelpString() const
{
return QObject::tr( "Creates a list of possible candidate coordinate reference systems for a layer "
"with an unknown projection.\n\n"
" The expected area which the layer should reside in must be specified via the "
"target area parameter.\n\n"
"The algorithm operates by testing the layer's extent in every known reference "
"system and listing any in which the bounds would fall near the target area if "
"the layer was in this projection." );
}
QgsFindProjectionAlgorithm *QgsFindProjectionAlgorithm::createInstance() const
{
return new QgsFindProjectionAlgorithm();
}
void QgsFindProjectionAlgorithm::initAlgorithm( const QVariantMap & )
{
addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ) ) );
addParameter( new QgsProcessingParameterExtent( QStringLiteral( "TARGET_AREA" ), QObject::tr( "Target area for layer" ) ) );
// deprecated
auto crsParam = std::make_unique<QgsProcessingParameterCrs>( QStringLiteral( "TARGET_AREA_CRS" ), QObject::tr( "Target area CRS" ), QVariant(), true );
crsParam->setFlags( crsParam->flags() | Qgis::ProcessingParameterFlag::Hidden );
addParameter( crsParam.release() );
addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "CRS candidates" ), Qgis::ProcessingSourceType::Vector ) );
}
QVariantMap QgsFindProjectionAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
std::unique_ptr<QgsProcessingFeatureSource> source( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
if ( !source )
throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
const QgsRectangle extent = parameterAsExtent( parameters, QStringLiteral( "TARGET_AREA" ), context );
QgsCoordinateReferenceSystem targetCrs = parameterAsExtentCrs( parameters, QStringLiteral( "TARGET_AREA" ), context );
if ( parameters.contains( QStringLiteral( "TARGET_AREA_CRS" ) ) )
{
QgsCoordinateReferenceSystem crs = parameterAsCrs( parameters, QStringLiteral( "TARGET_AREA_CRS" ), context );
if ( crs.isValid() )
{
targetCrs = crs;
}
}
QgsGeometry targetGeometry = QgsGeometry::fromRect( extent );
QgsFields fields;
fields.append( QgsField( "auth_id", QMetaType::Type::QString, "", 20 ) );
QString dest;
std::unique_ptr<QgsFeatureSink> sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, Qgis::WkbType::NoGeometry, QgsCoordinateReferenceSystem() ) );
if ( !sink )
{
throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
}
std::unique_ptr<QgsGeometryEngine> engine( QgsGeometry::createGeometryEngine( targetGeometry.constGet() ) );
engine->prepareGeometry();
QgsGeometry layerBounds = QgsGeometry::fromRect( source->sourceExtent() );
const QList< long > crsesToCheck = QgsCoordinateReferenceSystem::validSrsIds();
double step = crsesToCheck.count() > 0 ? 100.0 / crsesToCheck.count() : 0;
long foundResults = 0;
long i = 0;
QgsCoordinateTransformContext transformContext;
for ( long srsId : crsesToCheck )
{
if ( feedback->isCanceled() )
{
break;
}
QgsCoordinateReferenceSystem candidateCrs = QgsCoordinateReferenceSystem::fromSrsId( srsId );
if ( !candidateCrs.isValid() )
{
continue;
}
QgsCoordinateTransform transformCandidate = QgsCoordinateTransform( candidateCrs, targetCrs, transformContext );
transformCandidate.setBallparkTransformsAreAppropriate( true );
transformCandidate.disableFallbackOperationHandler( true );
QgsGeometry transformedBounds = QgsGeometry( layerBounds );
try
{
if ( transformedBounds.transform( transformCandidate ) != Qgis::GeometryOperationResult::Success )
{
continue;
}
}
catch ( QgsCsException & )
{
continue;
}
try
{
if ( engine->intersects( transformedBounds.constGet() ) )
{
feedback->pushInfo( QObject::tr( "Found candidate CRS: %1." ).arg( candidateCrs.authid() ) );
QgsFeature f = QgsFeature( fields );
f.setAttributes( QgsAttributes() << candidateCrs.authid() );
sink->addFeature( f, QgsFeatureSink::Flag::FastInsert );
foundResults++;
}
}
catch ( QgsCsException & )
{
continue;
}
feedback->setProgress( i * step );
i++;
}
if ( foundResults == 0 )
{
feedback->reportError( QObject::tr( "No matching projections found." ) );
}
sink->finalize();
QVariantMap outputs;
outputs.insert( QStringLiteral( "OUTPUT" ), dest );
return outputs;
}
///@endcond

50
src/analysis/processing/qgsalgorithmfindprojection.h Normal file
View File

@ -0,0 +1,50 @@
/***************************************************************************
qgsalgorithmfindprojection.h
---------------------
begin : March 2025
copyright : (C) 2025 by Alexander Bruy
email : alexander dot bruy 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. *
* *
***************************************************************************/
#ifndef QGSALGORITHMFINDPROJECTION_H
#define QGSALGORITHMFINDPROJECTION_H
#define SIP_NO_FILE
#include "qgis_sip.h"
#include "qgsprocessingalgorithm.h"
///@cond PRIVATE
/**
* Native find projection algorithm.
*/
class QgsFindProjectionAlgorithm : public QgsProcessingAlgorithm
{
public:
QgsFindProjectionAlgorithm() = default;
void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override;
QString name() const override;
QString displayName() const override;
QStringList tags() const override;
QString group() const override;
QString groupId() const override;
QString shortHelpString() const override;
QgsFindProjectionAlgorithm *createInstance() const override SIP_FACTORY;
protected:
QVariantMap processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
};
///@endcond PRIVATE
#endif // QGSALGORITHMFINDPROJECTION_H

2
src/analysis/processing/qgsnativealgorithms.cpp
View File

@ -132,6 +132,7 @@
#include "qgsalgorithmfilter.h"
#include "qgsalgorithmfilterbygeometry.h"
#include "qgsalgorithmfiltervertices.h"
#include "qgsalgorithmfindprojection.h"
#include "qgsalgorithmfixgeometries.h"
#include "qgsalgorithmflattenrelationships.h"
#include "qgsalgorithmforcerhr.h"
@ -449,6 +450,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
addAlgorithm( new QgsFilterByLayerTypeAlgorithm() );
addAlgorithm( new QgsFilterVerticesByM() );
addAlgorithm( new QgsFilterVerticesByZ() );
addAlgorithm( new QgsFindProjectionAlgorithm() );
addAlgorithm( new QgsFixGeometriesAlgorithm() );
addAlgorithm( new QgsFlattenRelationshipsAlgorithm() );
addAlgorithm( new QgsForceRHRAlgorithm() );