-
Notifications
You must be signed in to change notification settings - Fork 68
/
xyToLine.py
368 lines (337 loc) · 15 KB
/
xyToLine.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
"""
/***************************************************************************
* *
* 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. *
* *
***************************************************************************/
"""
import os
import math
from geographiclib.geodesic import Geodesic
from qgis.core import QgsCoordinateTransform, QgsPointXY, QgsFeature, QgsGeometry, QgsProject, QgsWkbTypes
from qgis.PyQt.QtGui import QIcon
from qgis.PyQt.QtCore import QUrl
from qgis.core import (
QgsProcessing,
QgsProcessingException,
QgsProcessingAlgorithm,
QgsProcessingParameterBoolean,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterCrs,
QgsProcessingParameterEnum,
QgsProcessingParameterField,
QgsProcessingParameterFeatureSink)
from .settings import settings, epsg4326, geod
from .utils import checkIdlCrossings, tr, GCgetPointsOnLine
# import traceback
class XYToLineAlgorithm(QgsProcessingAlgorithm):
"""
Algorithm for creating lines from two coordinates within a record.
"""
LINE_TYPE = ['Geodesic', 'Great Circle', 'Simple Line']
PrmInputLayer = 'InputLayer'
PrmOutputPointLayer = 'OutputPointLayer'
PrmOutputLineLayer = 'OutputLineLayer'
PrmInputCRS = 'InputCRS'
PrmOutputCRS = 'OutputCRS'
PrmLineType = 'LineType'
PrmStartUseLayerGeom = 'StartUseLayerGeom'
PrmStartXField = 'StartXField'
PrmStartYField = 'StartYField'
PrmEndUseLayerGeom = 'EndUseLayerGeom'
PrmEndXField = 'EndXField'
PrmEndYField = 'EndYField'
PrmShowStartPoint = 'ShowStartPoint'
PrmShowEndPoint = 'ShowEndPoint'
PrmDateLineBreak = 'DateLineBreak'
def initAlgorithm(self, config):
self.addParameter(
QgsProcessingParameterFeatureSource(
self.PrmInputLayer,
tr('Input layer'),
[QgsProcessing.TypeFile | QgsProcessing.TypeVectorPoint])
)
self.addParameter(
QgsProcessingParameterCrs(
self.PrmInputCRS,
tr('Input CRS for coordinates within the vector fields'),
'EPSG:4326')
)
self.addParameter(
QgsProcessingParameterCrs(
self.PrmOutputCRS,
tr('Output layer CRS'),
'EPSG:4326')
)
self.addParameter(
QgsProcessingParameterEnum(
self.PrmLineType,
tr('Line type'),
options=self.LINE_TYPE,
defaultValue=0,
optional=False)
)
self.addParameter(
QgsProcessingParameterBoolean(
self.PrmStartUseLayerGeom,
tr('Use the point geometry for the line starting point'),
False,
optional=True)
)
self.addParameter(
QgsProcessingParameterField(
self.PrmStartXField,
tr('Starting X Field (lon)'),
parentLayerParameterName=self.PrmInputLayer,
type=QgsProcessingParameterField.Any,
optional=True
)
)
self.addParameter(
QgsProcessingParameterField(
self.PrmStartYField,
tr('Starting Y Field (lat)'),
parentLayerParameterName=self.PrmInputLayer,
type=QgsProcessingParameterField.Any,
optional=True
)
)
self.addParameter(
QgsProcessingParameterBoolean(
self.PrmEndUseLayerGeom,
tr('Use the point geometry for the line ending point'),
False,
optional=True)
)
self.addParameter(
QgsProcessingParameterField(
self.PrmEndXField,
tr('Ending X Field (lon)'),
parentLayerParameterName=self.PrmInputLayer,
type=QgsProcessingParameterField.Any,
optional=True
)
)
self.addParameter(
QgsProcessingParameterField(
self.PrmEndYField,
tr('Ending Y Field (lat)'),
parentLayerParameterName=self.PrmInputLayer,
type=QgsProcessingParameterField.Any,
optional=True
)
)
self.addParameter(
QgsProcessingParameterBoolean(
self.PrmShowStartPoint,
tr('Show starting point'),
False,
optional=True)
)
self.addParameter(
QgsProcessingParameterBoolean(
self.PrmShowEndPoint,
tr('Show ending point'),
True,
optional=True)
)
self.addParameter(
QgsProcessingParameterBoolean(
self.PrmDateLineBreak,
tr('Break lines at -180, 180 boundary for better rendering'),
False,
optional=True)
)
self.addParameter(
QgsProcessingParameterFeatureSink(
self.PrmOutputLineLayer,
tr('Output line layer'))
)
self.addParameter(
QgsProcessingParameterFeatureSink(
self.PrmOutputPointLayer,
tr('Output point layer'),
optional=True,
createByDefault=True)
)
def processAlgorithm(self, parameters, context, feedback):
source = self.parameterAsSource(parameters, self.PrmInputLayer, context)
sourceCrs = self.parameterAsCrs(parameters, self.PrmInputCRS, context)
sinkCrs = self.parameterAsCrs(parameters, self.PrmOutputCRS, context)
lineType = self.parameterAsInt(parameters, self.PrmLineType, context)
startUseGeom = self.parameterAsBool(parameters, self.PrmStartUseLayerGeom, context)
startXcol = self.parameterAsString(parameters, self.PrmStartXField, context)
startYcol = self.parameterAsString(parameters, self.PrmStartYField, context)
endUseGeom = self.parameterAsBool(parameters, self.PrmEndUseLayerGeom, context)
endXcol = self.parameterAsString(parameters, self.PrmEndXField, context)
endYcol = self.parameterAsString(parameters, self.PrmEndYField, context)
showStart = self.parameterAsBool(parameters, self.PrmShowStartPoint, context)
showEnd = self.parameterAsBool(parameters, self.PrmShowEndPoint, context)
dateLine = self.parameterAsBool(parameters, self.PrmDateLineBreak, context)
if startUseGeom and endUseGeom:
msg = tr('The layer geometry cannot be used for both the starting and ending points.')
raise QgsProcessingException(msg)
if (startUseGeom or endUseGeom) and (source.wkbType() != QgsWkbTypes.Point):
msg = tr('In order to use the layer geometry for the start or ending points, the input layer must be of type Point')
raise QgsProcessingException(msg)
if (not startUseGeom and (not startXcol or not startYcol)) or (not endUseGeom and (not endXcol or not endYcol)):
msg = tr('Please select valid starting and ending point columns')
raise QgsProcessingException(msg)
if dateLine and lineType <= 1:
isMultiPart = True
else:
isMultiPart = False
if isMultiPart:
(lineSink, lineDest_id) = self.parameterAsSink(
parameters, self.PrmOutputLineLayer, context, source.fields(),
QgsWkbTypes.MultiLineString, sinkCrs)
else:
(lineSink, lineDest_id) = self.parameterAsSink(
parameters, self.PrmOutputLineLayer, context, source.fields(),
QgsWkbTypes.LineString, sinkCrs)
skip_pt = True if self.PrmOutputPointLayer not in parameters or parameters[self.PrmOutputPointLayer] is None else False
if (showStart or showEnd) and not skip_pt:
(ptSink, ptDest_id) = self.parameterAsSink(
parameters, self.PrmOutputPointLayer, context, source.fields(),
QgsWkbTypes.Point, sinkCrs)
else:
if showStart or showEnd:
feedback.pushInfo(tr('Output point layer was set to [skip output]. No point layer will be generated.'))
showStart = False
showEnd = False
else:
feedback.pushInfo(tr('No beginning or ending points were selected so a point layer will not be generated.'))
# Set up CRS transformations
geomCrs = source.sourceCrs()
if (startUseGeom or endUseGeom) and (geomCrs != epsg4326):
geomTo4326 = QgsCoordinateTransform(geomCrs, epsg4326, QgsProject.instance())
if sourceCrs != epsg4326:
sourceTo4326 = QgsCoordinateTransform(sourceCrs, epsg4326, QgsProject.instance())
if sinkCrs != epsg4326:
toSinkCrs = QgsCoordinateTransform(epsg4326, sinkCrs, QgsProject.instance())
featureCount = source.featureCount()
total = 100.0 / featureCount if featureCount else 0
numBad = 0
maxseglen = settings.maxSegLength * 1000.0
maxSegments = settings.maxSegments
beginning_ending_same = False
iterator = source.getFeatures()
for cnt, feature in enumerate(iterator):
if (cnt % 100 == 0) and feedback.isCanceled():
break
try:
if startUseGeom:
ptStart = feature.geometry().asPoint()
if geomCrs != epsg4326:
ptStart = geomTo4326.transform(ptStart)
else:
ptStart = QgsPointXY(float(feature[startXcol]), float(feature[startYcol]))
if sourceCrs != epsg4326:
ptStart = sourceTo4326.transform(ptStart)
if endUseGeom:
ptEnd = feature.geometry().asPoint()
if geomCrs != epsg4326:
ptEnd = geomTo4326.transform(ptEnd)
else:
ptEnd = QgsPointXY(float(feature[endXcol]), float(feature[endYcol]))
if sourceCrs != epsg4326:
ptEnd = sourceTo4326.transform(ptEnd)
pts = [ptStart]
if ptStart == ptEnd: # We cannot have a line that begins and ends at the same point
numBad += 1
beginning_ending_same = True
continue
if lineType == 0: # Geodesic
gline = geod.InverseLine(ptStart.y(), ptStart.x(), ptEnd.y(), ptEnd.x())
if gline.s13 > maxseglen:
n = int(math.ceil(gline.s13 / maxseglen))
if n > maxSegments:
n = maxSegments
seglen = gline.s13 / n
for i in range(1, n + 1):
s = seglen * i
g = gline.Position(s, Geodesic.LATITUDE | Geodesic.LONGITUDE)
pts.append(QgsPointXY(g['lon2'], g['lat2']))
else: # The line segment is too short so it is from ptStart to ptEnd
pts.append(ptEnd)
elif lineType == 1: # Great circle
pts = GCgetPointsOnLine(
ptStart.y(), ptStart.x(),
ptEnd.y(), ptEnd.x(),
settings.maxSegLength * 1000.0, # Put it in meters
settings.maxSegments + 1)
else: # Simple line
pts.append(ptEnd)
f = QgsFeature()
if isMultiPart:
outseg = checkIdlCrossings(pts)
if sinkCrs != epsg4326: # Convert each point to the output CRS
for y in range(len(outseg)):
for x, pt in enumerate(outseg[y]):
outseg[y][x] = toSinkCrs.transform(pt)
f.setGeometry(QgsGeometry.fromMultiPolylineXY(outseg))
else:
if sinkCrs != epsg4326: # Convert each point to the output CRS
for x, pt in enumerate(pts):
pts[x] = toSinkCrs.transform(pt)
f.setGeometry(QgsGeometry.fromPolylineXY(pts))
f.setAttributes(feature.attributes())
lineSink.addFeature(f)
if showStart:
f = QgsFeature()
if sinkCrs != epsg4326:
f.setGeometry(QgsGeometry.fromPointXY(toSinkCrs.transform(ptStart)))
else:
f.setGeometry(QgsGeometry.fromPointXY(ptStart))
f.setAttributes(feature.attributes())
ptSink.addFeature(f)
if showEnd:
f = QgsFeature()
if sinkCrs != epsg4326:
f.setGeometry(QgsGeometry.fromPointXY(toSinkCrs.transform(ptEnd)))
else:
f.setGeometry(QgsGeometry.fromPointXY(ptEnd))
f.setAttributes(feature.attributes())
ptSink.addFeature(f)
except Exception:
numBad += 1
'''s = traceback.format_exc()
feedback.pushInfo(s)'''
if cnt % 100 == 0: # Set the progress after every 100 entries
feedback.setProgress(int(cnt * total))
if numBad > 0:
if beginning_ending_same:
feedback.pushInfo(tr("One of more features had the same beginning and ending coordinate and are invalid."))
feedback.pushInfo(tr("{} out of {} features from the input layer were invalid and were ignored.".format(numBad, featureCount)))
r = {}
r[self.PrmOutputLineLayer] = lineDest_id
if showStart or showEnd:
r[self.PrmOutputPointLayer] = ptDest_id
return (r)
def name(self):
return 'xy2line'
def icon(self):
return QIcon(os.path.dirname(__file__) + '/images/xyline.svg')
def displayName(self):
return tr('XY to line')
def group(self):
return tr('Vector geometry')
def groupId(self):
return 'vectorgeometry'
def helpUrl(self):
file = os.path.dirname(__file__) + '/index.html'
if not os.path.exists(file):
return ''
return QUrl.fromLocalFile(file).toString(QUrl.FullyEncoded)
def shortHelpString(self):
file = os.path.dirname(__file__) + '/doc/XYtoLineAlgorithm.help'
if not os.path.exists(file):
return ''
with open(file) as helpf:
help = helpf.read()
return help
def createInstance(self):
return XYToLineAlgorithm()