How to use the geographiclib.geodesic.Geodesic.LATITUDE 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 / ext-libs / geographiclib / polygonarea.py View on Github external
"""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()
github NationalSecurityAgency / qgis-shapetools-plugin / lineDigitizer.py View on Github external
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()
github NationalSecurityAgency / qgis-shapetools-plugin / azDigitizer.py View on Github external
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()
github NationalSecurityAgency / qgis-shapetools-plugin / vector2Shape.py View on Github external
# 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)
github NationalSecurityAgency / qgis-shapetools-plugin / createHeart.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())
            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)
github NationalSecurityAgency / qgis-shapetools-plugin / createLob.py View on Github external
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:
github NationalSecurityAgency / qgis-shapetools-plugin / geodesicMeasureTool.py View on Github external
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
github NationalSecurityAgency / qgis-shapetools-plugin / geodesicDensify.py View on Github external
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):