Skip to content

Commit

Permalink
[FEATURE][processing] Algorithm to find an unknown layer's projection
Browse files Browse the repository at this point in the history 
If you have a layer with an unknown CRS, this algorithm gives a list
of possible candidate CRSes which the layer could be in.

It allows users to set the area (and corresponding CRS) which they know
the layer should be located near. The algorithm then tests every CRS
in the database to see what candidate CRSes would cause the layer
to be located at that preset area.

It's much faster than it sounds!! (just a couple of seconds)

Sponsored by SMEC/Surbana Jurong
  • Loading branch information
@nyalldawson
nyalldawson committed Feb 13, 2017
1 parent 3ecafb3 commit 91eea30
Show file tree
Hide file tree
Showing 7 changed files with 200 additions and 1 deletion.
7 changes: 7 additions & 0 deletions python/plugins/processing/algs/help/qgis.yaml
View file
Expand Up @@ -200,6 +200,13 @@ qgis:extractspecificnodes: >
qgis:fieldcalculator: >
This algorithm computes a new vector layer with the same features of the input layer, but with an additional attribute. The values of this new attribute are computed from each feature using a mathematical formula, based on te properties and attributes of the feature.

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. Additionally, the coordinate reference system for this target area must also be set.

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:fixeddistancebuffer: >
This algorithm computes a buffer area for all the features in an input layer, using a fixed distance.

Expand Down
116 changes: 116 additions & 0 deletions python/plugins/processing/algs/qgis/FindProjection.py
View file
@@ -0,0 +1,116 @@
# -*- coding: utf-8 -*-

"""
***************************************************************************
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'

# This will get replaced with a git SHA1 when you do a git archive

__revision__ = '$Format:%H$'

import os
import codecs

from qgis.core import (QgsGeometry,
QgsRectangle,
QgsCoordinateReferenceSystem,
QgsCoordinateTransform)

from processing.core.GeoAlgorithm import GeoAlgorithm
from processing.core.parameters import ParameterVector
from processing.core.parameters import ParameterCrs
from processing.core.parameters import ParameterExtent
from processing.core.outputs import OutputHTML
from processing.tools import dataobjects


pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]


class FindProjection(GeoAlgorithm):

INPUT_LAYER = 'INPUT_LAYER'
TARGET_AREA = 'TARGET_AREA'
TARGET_AREA_CRS = 'TARGET_AREA_CRS'
OUTPUT_HTML_FILE = 'OUTPUT_HTML_FILE'

def defineCharacteristics(self):
self.name, self.i18n_name = self.trAlgorithm('Find projection')
self.group, self.i18n_group = self.trAlgorithm('Vector general tools')
self.tags = self.tr('crs,srs,coordinate,reference,system,guess,estimate,finder,determine')

self.addParameter(ParameterVector(self.INPUT_LAYER,
self.tr('Input layer')))
extent_parameter = ParameterExtent(self.TARGET_AREA,
self.tr('Target area for layer'),
self.INPUT_LAYER)
extent_parameter.skip_crs_check = True
self.addParameter(extent_parameter)
self.addParameter(ParameterCrs(self.TARGET_AREA_CRS, 'Target area CRS'))

self.addOutput(OutputHTML(self.OUTPUT_HTML_FILE,
self.tr('Candidates')))

def processAlgorithm(self, feedback):
layer = dataobjects.getObjectFromUri(
self.getParameterValue(self.INPUT_LAYER))

extent = self.getParameterValue(self.TARGET_AREA).split(',')
target_crs = QgsCoordinateReferenceSystem(self.getParameterValue(self.TARGET_AREA_CRS))

target_geom = QgsGeometry.fromRect(QgsRectangle(float(extent[0]), float(extent[2]),
float(extent[1]), float(extent[3])))

output_file = self.getOutputValue(self.OUTPUT_HTML_FILE)

# make intersection tests nice and fast
engine = QgsGeometry.createGeometryEngine(target_geom.geometry())
engine.prepareGeometry()

layer_bounds = QgsGeometry.fromRect(layer.extent())

results = []

for srs_id in QgsCoordinateReferenceSystem.validSrsIds():
candidate_crs = QgsCoordinateReferenceSystem.fromSrsId(srs_id)
if not candidate_crs.isValid():
continue

transform_candidate = QgsCoordinateTransform(candidate_crs, target_crs)
transformed_bounds = QgsGeometry(layer_bounds)
try:
if not transformed_bounds.transform(transform_candidate) == 0:
continue
except:
continue

if engine.intersects(transformed_bounds.geometry()):
results.append(candidate_crs.authid())

self.createHTML(output_file, results)

def createHTML(self, outputFile, candidates):
with codecs.open(outputFile, 'w', encoding='utf-8') as f:
f.write('<html><head>\n')
f.write('<meta http-equiv="Content-Type" content="text/html; \
charset=utf-8" /></head><body>\n')
for c in candidates:
f.write('<p>' + c + '</p>\n')
f.write('</body></html>\n')
3 changes: 2 additions & 1 deletion python/plugins/processing/algs/qgis/QGISAlgorithmProvider.py
View file
Expand Up @@ -186,6 +186,7 @@
from .Polygonize import Polygonize
from .FixGeometry import FixGeometry
from .ExecuteSQL import ExecuteSQL
from .FindProjection import FindProjection

pluginPath = os.path.normpath(os.path.join(
os.path.split(os.path.dirname(__file__))[0], os.pardir))
Expand Down Expand Up @@ -255,7 +256,7 @@ def __init__(self):
ShortestPathPointToPoint(), ShortestPathPointToLayer(),
ShortestPathLayerToPoint(), ServiceAreaFromPoint(),
ServiceAreaFromLayer(), TruncateTable(), Polygonize(),
FixGeometry(), ExecuteSQL()
FixGeometry(), ExecuteSQL(), FindProjection()
]

if hasMatplotlib:
Expand Down
16 changes: 16 additions & 0 deletions python/plugins/processing/tests/testdata/custom/find_projection.gfs
View file
@@ -0,0 +1,16 @@
<GMLFeatureClassList>
<GMLFeatureClass>
<Name>find_projection</Name>
<ElementPath>find_projection</ElementPath>
<!--POINT-->
<GeometryType>1</GeometryType>
<SRSName>EPSG:28356</SRSName>
<DatasetSpecificInfo>
<FeatureCount>4</FeatureCount>
<ExtentXMin>326064.87345</ExtentXMin>
<ExtentXMax>328710.85361</ExtentXMax>
<ExtentYMin>6245538.80829</ExtentYMin>
<ExtentYMax>6247808.84283</ExtentYMax>
</DatasetSpecificInfo>
</GMLFeatureClass>
</GMLFeatureClassList>
34 changes: 34 additions & 0 deletions python/plugins/processing/tests/testdata/custom/find_projection.gml
View file
@@ -0,0 +1,34 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation=""
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>326064.8734538725</gml:X><gml:Y>6245538.808291084</gml:Y></gml:coord>
<gml:coord><gml:X>328710.8536064086</gml:X><gml:Y>6247808.842832778</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>

<gml:featureMember>
<ogr:find_projection fid="find_projection.0">
<ogr:geometryProperty><gml:Point srsName="EPSG:28356"><gml:coordinates>326597.10210384,6247528.48897066</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:find_projection>
</gml:featureMember>
<gml:featureMember>
<ogr:find_projection fid="find_projection.1">
<ogr:geometryProperty><gml:Point srsName="EPSG:28356"><gml:coordinates>327449.712845641,6246290.065952</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:find_projection>
</gml:featureMember>
<gml:featureMember>
<ogr:find_projection fid="find_projection.2">
<ogr:geometryProperty><gml:Point srsName="EPSG:28356"><gml:coordinates>328710.853606409,6247808.84283278</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:find_projection>
</gml:featureMember>
<gml:featureMember>
<ogr:find_projection fid="find_projection.3">
<ogr:geometryProperty><gml:Point srsName="EPSG:28356"><gml:coordinates>326064.873453873,6245538.80829108</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:find_projection>
</gml:featureMember>
</ogr:FeatureCollection>
11 changes: 11 additions & 0 deletions python/plugins/processing/tests/testdata/expected/find_projection.html
View file
@@ -0,0 +1,11 @@
<html><head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" /></head><body>
<p>EPSG:20256</p>
<p>EPSG:20356</p>
<p>EPSG:28356</p>
<p>EPSG:32356</p>
<p>EPSG:32556</p>
<p>EPSG:32756</p>
<p>IGNF:UTM56SW84</p>
<p>EPSG:5552</p>
</body></html>
14 changes: 14 additions & 0 deletions python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml
View file
Expand Up @@ -2347,3 +2347,17 @@ tests:
OUTPUT:
name: expected/voronoi_buffer.gml
type: vector

- algorithm: qgis:findprojection
name: Find projection
params:
INPUT_LAYER:
name: custom/find_projection.gml
type: vector
TARGET_AREA: 151.1198,151.1368,-33.9118,-33.9003
TARGET_AREA_CRS: EPSG:4326
results:
OUTPUT_HTML_FILE:
name: expected/find_projection.html
type: file

0 comments on commit 91eea30

Please sign in to comment.