Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
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())
arange = 360.0 / k
angle = -arange / 2.0
astep = arange / cnt
for i in range(k):
aoffset = arange * (k - 1)
index = 0
while index < cnt:
r = dist[index] * radius
g = geod.Direct(pt.y(), pt.x(), angle + aoffset + startAngle, r, Geodesic.LATITUDE | Geodesic.LONGITUDE)
pts.append(QgsPointXY(g['lon2'], g['lat2']))
angle += astep
index += 1
# repeat the very first point to close the polygon
pts.append(pts[0])
# 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))
continue
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)
x = r * (cusps2 - 1.0) * math.cos(a) + r * math.cos((cusps2 - 1.0) * a)
y = r * (cusps2 - 1.0) * math.sin(a) - r * math.sin((cusps2 - 1.0) * a)
a2 = math.degrees(math.atan2(y, x)) + sangle
dist = math.sqrt(x * x + y * y)
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)
ptEnd = sourceTo4326.transform(ptEnd)
pts = [ptStart]
if ptStart == ptEnd: # We cannot have a line that begins and ends at the same point
numBad += 1
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)):
if layercrs != epsg4326: # Convert to 4326
ptStart = transto4326.transform(ptStart)
pts = [ptStart]
for x in range(1,numpoints):
ptEnd = QgsPoint(points[x][0], points[x][1])
if layercrs != epsg4326: # Convert to 4326
ptEnd = transto4326.transform(ptEnd)
l = self.geod.InverseLine(ptStart.y(), ptStart.x(), ptEnd.y(), ptEnd.x())
n = int(math.ceil(l.s13 / maxseglen))
if n > maxSegments:
n = maxSegments
seglen = l.s13 / n
for i in range(1,n):
s = seglen * i
g = l.Position(s, Geodesic.LATITUDE | Geodesic.LONGITUDE | Geodesic.LONG_UNROLL)
pts.append( QgsPoint(g['lon2'], g['lat2']) )
pts.append(ptEnd)
ptStart = ptEnd
if layercrs != epsg4326: # Convert each point to the output CRS
for x, pt in enumerate(pts):
pts[x] = transfrom4326.transform(pt)
ptset.append(pts)
if len(ptset) > 0:
featureout = QgsFeature()
featureout.setGeometry(QgsGeometry.fromPolygon(ptset))
featureout.setAttributes(feature.attributes())
ppoly.addFeatures([featureout])
else:
returned. The :class:`~geographiclib.geodesicline.GeodesicLine`
object must have been constructed with the DISTANCE_IN capability.
"""
from geographiclib.geodesic import Geodesic
result = {'lat1': self.lat1,
'lon1': self.lon1 if outmask & Geodesic.LONG_UNROLL else
Math.AngNormalize(self.lon1),
'azi1': self.azi1, 's12': s12}
a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = self._GenPosition(
False, s12, outmask)
outmask &= Geodesic.OUT_MASK
result['a12'] = a12
if outmask & Geodesic.LATITUDE: result['lat2'] = lat2
if outmask & Geodesic.LONGITUDE: result['lon2'] = lon2
if outmask & Geodesic.AZIMUTH: result['azi2'] = azi2
if outmask & Geodesic.REDUCEDLENGTH: result['m12'] = m12
if outmask & Geodesic.GEODESICSCALE:
result['M12'] = M12; result['M21'] = M21
if outmask & Geodesic.AREA: result['S12'] = S12
return result
# If the input is not 4326 we need to convert it to that and then back to the output CRS
ptStart = QgsPointXY(points[0][0], points[0][1])
if layercrs != epsg4326: # Convert to 4326
ptStart = transto4326.transform(ptStart)
pts = [ptStart]
for x in range(1, numpoints):
ptEnd = QgsPointXY(points[x][0], points[x][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)
ptStart = ptEnd
if layercrs != epsg4326: # Convert each point to the output CRS
for x, pt in enumerate(pts):
pts[x] = transfrom4326.transform(pt)
ptset.append(pts)
multiset.append(ptset)
if len(multiset) > 0:
featureout = QgsFeature()
featureout.setGeometry(QgsGeometry.fromMultiPolygonXY(multiset))
featureout.setAttributes(feature.attributes())
sink.addFeature(featureout)
ab = sma * smi
step = 18.0 * smi / sma
if step < 1.0:
minimum = step
else:
minimum = 1.0
maxang = math.pi / 6 * minimum
delta = ab * math.pi / segments
pts = []
azi = 0
while azi < math.tau:
cos_azi = math.cos(azi)
sin_azi = math.sin(azi)
rad = ab / math.sqrt(sma * sma * sin_azi * sin_azi + smi * smi * cos_azi * cos_azi)
g = geod.Direct(lat, lon, math.degrees(azi) + orient, rad, Geodesic.LATITUDE | Geodesic.LONGITUDE)
pts.append(QgsPointXY(g['lon2'], g['lat2']))
delo = delta / (rad * rad)
if maxang < delo:
delo = maxang
azi += delo
# Append the starting point to close the shape
pts.append(pts[0])
return(pts)
else:
s = sides
if anglecol:
startangle = float(feature[anglecol])
else:
startangle = angle
if distcol:
d = float(feature[distcol]) * measureFactor
else:
d = defaultDist
pts = []
i = s
while i >= 0:
a = (i * 360.0 / s) + startangle
i -= 1
g = geod.Direct(pt.y(), pt.x(), a, d, Geodesic.LATITUDE | Geodesic.LONGITUDE)
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()
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)
attr.append(pt_orig_y)
B12 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C1a)
AB1 = (1 + self._A1m1) * (B12 - self._B11)
# sin(bet2) = cos(alp0) * sin(sig2)
sbet2 = self._calp0 * ssig2
# Alt: cbet2 = hypot(csig2, salp0 * ssig2)
cbet2 = math.hypot(self._salp0, self._calp0 * csig2)
if cbet2 == 0:
# I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
cbet2 = csig2 = Geodesic.tiny_
# tan(alp0) = cos(sig2)*tan(alp2)
salp2 = self._salp0; calp2 = self._calp0 * csig2 # No need to normalize
if outmask & Geodesic.DISTANCE:
s12 = self._b * ((1 + self._A1m1) * sig12 + AB1) if arcmode else s12_a12
if outmask & Geodesic.LONGITUDE:
# tan(omg2) = sin(alp0) * tan(sig2)
somg2 = self._salp0 * ssig2; comg2 = csig2 # No need to normalize
E = Math.copysign(1, self._salp0) # East or west going?
# omg12 = omg2 - omg1
omg12 = (E * (sig12
- (math.atan2( ssig2, csig2) -
math.atan2( self._ssig1, self._csig1))
+ (math.atan2(E * somg2, comg2) -
math.atan2(E * self._somg1, self._comg1)))
if outmask & Geodesic.LONG_UNROLL
else math.atan2(somg2 * self._comg1 - comg2 * self._somg1,
comg2 * self._comg1 + somg2 * self._somg1))
lam12 = omg12 + self._A3c * (
sig12 + (Geodesic._SinCosSeries(True, ssig2, csig2, self._C3a)
- self._B31))
lon12 = math.degrees(lam12)
if startAngleCol:
sangle = float(feature[startAngleCol])
else:
sangle = startAngle
if starPointsCol:
spoints = float(feature[starPointsCol])
shalf = (360.0 / spoints) / 2.0
else:
spoints = numPoints
shalf = half
i = spoints - 1
while i >= 0:
i -= 1
angle = (i * 360.0 / spoints) + sangle
g = geod.Direct(pt.y(), pt.x(), angle, oradius, Geodesic.LATITUDE | Geodesic.LONGITUDE)
pts.append(QgsPointXY(g['lon2'], g['lat2']))
g = geod.Direct(pt.y(), pt.x(), angle - shalf, iradius, Geodesic.LATITUDE | Geodesic.LONGITUDE)
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()
if shapetype == 0:
f.setGeometry(QgsGeometry.fromPolygonXY([pts]))
else:
f.setGeometry(QgsGeometry.fromPolylineXY(pts))
attr = feature.attributes()
if export_geom: