How to use the geographiclib.geodesic.Geodesic.LONGITUDE function in geographiclib

To help you get started, we’ve selected a few geographiclib examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github NationalSecurityAgency / qgis-shapetools-plugin / createRose.py View on Github external
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))
github NationalSecurityAgency / qgis-shapetools-plugin / createHypocycloid.py View on Github external
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)
github NationalSecurityAgency / qgis-shapetools-plugin / xyToLine.py View on Github external
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)):
github NationalSecurityAgency / qgis-shapetools-plugin / poly2Geodesic.py View on Github external
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:
github NREL / OpenStudio / openstudiocore / src / geographic_lib / python / geographiclib / geodesicline.py View on Github external
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
github NationalSecurityAgency / qgis-shapetools-plugin / geodesicDensify.py View on Github external
# 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)
github NationalSecurityAgency / qgis-shapetools-plugin / createEllipse.py View on Github external
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)
github NationalSecurityAgency / qgis-shapetools-plugin / createPolygon.py View on Github external
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)
github NationalSecurityAgency / qgis-shapetools-plugin / ext-libs / geographiclib / geodesicline.py View on Github external
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)
github NationalSecurityAgency / qgis-shapetools-plugin / createStar.py View on Github external
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: