Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"""Construct a PolygonArea object
:param earth: a :class:`~geographiclib.geodesic.Geodesic` object
:param polyline: if true, treat object as a polyline instead of a polygon
Initially the polygon has no vertices.
"""
from geographiclib.geodesic import Geodesic
self.earth = earth
"""The geodesic object (readonly)"""
self.polyline = polyline
"""Is this a polyline? (readonly)"""
self.area0 = 4 * math.pi * earth._c2
"""The total area of the ellipsoid in meter^2 (readonly)"""
self._mask = (Geodesic.LATITUDE | Geodesic.LONGITUDE |
Geodesic.DISTANCE |
(Geodesic.EMPTY if self.polyline else
Geodesic.AREA | Geodesic.LONG_UNROLL))
if not self.polyline: self._areasum = Accumulator()
self._perimetersum = Accumulator()
self.num = 0
"""The current number of points in the polygon (readonly)"""
self.lat1 = Math.nan
"""The current latitude in degrees (readonly)"""
self.lon1 = Math.nan
"""The current longitude in degrees (readonly)"""
self.Clear()
for x in range(numpairs):
azimuth = values[x << 1] + declination
distance = values[(x << 1) + 1] * measureFactor
g = geod.Direct(pt.y(), pt.x(), azimuth, distance, Geodesic.LATITUDE | Geodesic.LONGITUDE)
pt = QgsPoint(g['lon2'], g['lat2']) # Keep this in EPSG:4326
pt_trans = transform.transform(g['lon2'], g['lat2']) # Transformed version
feat = QgsFeature(layer.fields())
feat.setGeometry(QgsGeometry.fromPointXY(pt_trans))
layer.addFeature(feat)
else: # It will either be a Line or Polygon
ptStart = transform.transform(self.pt.x(), self.pt.y())
pts = [ptStart]
for x in range(numpairs):
azimuth = values[x << 1] + declination
distance = values[(x << 1) + 1] * measureFactor
g = geod.Direct(pt.y(), pt.x(), azimuth, distance, Geodesic.LATITUDE | Geodesic.LONGITUDE)
pt = QgsPoint(g['lon2'], g['lat2']) # Keep this in EPSG:4326
pt_trans = transform.transform(g['lon2'], g['lat2']) # Transformed version
pts.append(pt_trans)
if layer.geometryType() == QgsWkbTypes.PolygonGeometry or closeline:
pts.append(ptStart)
feat = QgsFeature(layer.fields())
if layer.geometryType() == QgsWkbTypes.LineGeometry:
feat.setGeometry(QgsGeometry.fromPolylineXY(pts))
else:
feat.setGeometry(QgsGeometry.fromPolygonXY([pts]))
layer.addFeatures([feat])
layer.updateExtents()
self.iface.mapCanvas().refresh()
self.close()
pt = transform.transform(g['lon2'], g['lat2'])
feat = QgsFeature(layer.fields())
feat.setGeometry(QgsGeometry.fromPointXY(pt))
layer.addFeature(feat)
else: # It will either be a LineString or MultiLineString
maxseglen = settings.maxSegLength * 1000.0 # Needs to be in meters
maxSegments = settings.maxSegments
gline = geod.Line(pt.y(), pt.x(), azimuth)
n = int(math.ceil(distance / maxseglen))
if n > maxSegments:
n = maxSegments
seglen = distance / n
pts = []
for i in range(0, n + 1):
s = seglen * i
g = gline.Position(s, Geodesic.LATITUDE | Geodesic.LONGITUDE | Geodesic.LONG_UNROLL)
ptc = transform.transform(g['lon2'], g['lat2'])
pts.append(ptc)
feat = QgsFeature(layer.fields())
if layer.wkbType() == QgsWkbTypes.LineString:
feat.setGeometry(QgsGeometry.fromPolylineXY(pts))
else:
feat.setGeometry(QgsGeometry.fromMultiPolylineXY([pts]))
layer.addFeatures([feat])
layer.updateExtents()
self.iface.mapCanvas().refresh()
self.close()
# make sure the coordinates are in EPSG:4326
pt = self.transform.transform(pt.x(), pt.y())
lat = pt.y()
lon = pt.x()
angle = 0
while angle < 360:
if innerCol == -1:
iRadius = defInnerRadius
else:
iRadius = float(feature[innerCol]) * measureFactor
if outerCol == -1:
oRadius = defOuterRadius
else:
oRadius = float(feature[outerCol]) * measureFactor
if iRadius != 0:
g = geod.Direct(lat, lon, angle, iRadius, Geodesic.LATITUDE | Geodesic.LONGITUDE)
ptsi.append(QgsPointXY(g['lon2'], g['lat2']))
g = geod.Direct(lat, lon, angle, oRadius, Geodesic.LATITUDE | Geodesic.LONGITUDE)
ptso.append(QgsPointXY(g['lon2'], g['lat2']))
angle += ptSpacing
if iRadius != 0:
ptsi.append(ptsi[0])
ptso.append(ptso[0])
# If the Output crs is not 4326 transform the points to the proper crs
if self.outputCRS != epsg4326:
if iRadius != 0:
for x, ptout in enumerate(ptsi):
ptsi[x] = self.transformOut.transform(ptout)
for x, ptout in enumerate(ptso):
ptso[x] = self.transformOut.transform(ptout)
pts = []
pt = feature.geometry().asPoint()
pt_orig_x = pt.x()
pt_orig_y = pt.y()
# make sure the coordinates are in EPSG:4326
if srcCRS != epsg4326:
pt = geomTo4326.transform(pt.x(), pt.y())
angle = 0.0
while angle <= 360.0:
a = math.radians(angle)
sina = math.sin(a)
x = 16 * sina * sina * sina
y = 13 * math.cos(a) - 5 * math.cos(2 * a) - 2 * math.cos(3 * a) - math.cos(4 * a)
dist = math.sqrt(x * x + y * y) * radius2 / 17.0
a2 = math.degrees(math.atan2(y, x)) + sangle
g = geod.Direct(pt.y(), pt.x(), a2, dist, Geodesic.LATITUDE | Geodesic.LONGITUDE)
pts.append(QgsPointXY(g['lon2'], g['lat2']))
angle += step
# If the Output crs is not 4326 transform the points to the proper crs
if srcCRS != epsg4326:
for x, ptout in enumerate(pts):
pts[x] = toSinkCrs.transform(ptout)
f = QgsFeature()
if shapetype == 0:
f.setGeometry(QgsGeometry.fromPolygonXY([pts]))
else:
f.setGeometry(QgsGeometry.fromPolylineXY(pts))
attr = feature.attributes()
if export_geom:
attr.append(pt_orig_x)
distance = defaultDist
pt = feature.geometry().asPoint()
pt_orig_x = pt.x()
pt_orig_y = pt.y()
# make sure the coordinates are in EPSG:4326
if srcCRS != epsg4326:
pt = geomTo4326.transform(pt.x(), pt.y())
pts = [pt]
gline = geod.Line(pt.y(), pt.x(), bearing)
n = int(math.ceil(distance / maxseglen))
if n > maxSegments:
n = maxSegments
seglen = distance / n
for i in range(1, n + 1):
s = seglen * i
g = gline.Position(s, Geodesic.LATITUDE | Geodesic.LONGITUDE | Geodesic.LONG_UNROLL)
pts.append(QgsPointXY(g['lon2'], g['lat2']))
# If the Output crs is not 4326 transform the points to the proper crs
if srcCRS != epsg4326:
for x, ptout in enumerate(pts):
pts[x] = toSinkCrs.transform(ptout)
f = QgsFeature()
f.setGeometry(QgsGeometry.fromPolylineXY(pts))
attr = feature.attributes()
if export_geom:
attr.append(pt_orig_x)
attr.append(pt_orig_y)
f.setAttributes(attr)
sink.addFeature(f)
except Exception:
def getLinePts(self, distance, pt1, pt2):
canvasCrs = self.canvas.mapSettings().destinationCrs()
transform = QgsCoordinateTransform(epsg4326, canvasCrs, QgsProject.instance())
pt1c = transform.transform(pt1.x(), pt1.y())
pt2c = transform.transform(pt2.x(), pt2.y())
if distance < 10000:
return [pt1c, pt2c]
gline = geod.InverseLine(pt1.y(), pt1.x(), pt2.y(), pt2.x())
n = int(math.ceil(distance / 10000.0))
if n > 20:
n = 20
seglen = distance / n
pts = [pt1c]
for i in range(1, n):
s = seglen * i
g = gline.Position(s, Geodesic.LATITUDE | Geodesic.LONGITUDE | Geodesic.LONG_UNROLL)
ptc = transform.transform(g['lon2'], g['lat2'])
pts.append(ptc)
pts.append(pt2c)
return pts
if discardVertices:
ptStart = QgsPointXY(seg[0][0][0], seg[0][0][1])
if layercrs != epsg4326: # Convert to 4326
ptStart = transto4326.transform(ptStart)
pts = [ptStart]
numpoints = len(seg[numseg - 1])
ptEnd = QgsPointXY(seg[numseg - 1][numpoints - 1][0], seg[numseg - 1][numpoints - 1][1])
if layercrs != epsg4326: # Convert to 4326
ptEnd = transto4326.transform(ptEnd)
gline = geod.InverseLine(ptStart.y(), ptStart.x(), ptEnd.y(), ptEnd.x())
if gline.s13 > maxseglen:
n = int(math.ceil(gline.s13 / maxseglen))
seglen = gline.s13 / n
for i in range(1, n):
s = seglen * i
g = gline.Position(s, Geodesic.LATITUDE | Geodesic.LONGITUDE | Geodesic.LONG_UNROLL)
pts.append(QgsPointXY(g['lon2'], g['lat2']))
pts.append(ptEnd)
if layercrs != epsg4326: # Convert each point back to the output CRS
for x, pt in enumerate(pts):
pts[x] = transfrom4326.transform(pt)
fline.setGeometry(QgsGeometry.fromPolylineXY(pts))
else:
if not feature.geometry().isMultipart():
line = seg[0]
numpoints = len(line)
ptStart = QgsPointXY(line[0][0], line[0][1])
if layercrs != epsg4326: # Convert to 4326
ptStart = transto4326.transform(ptStart)
pts = [ptStart]
for x in range(1, numpoints):