Python QgsDistanceArea.ellipsoid方法代码示例

# 需要导入模块: from qgis.core import QgsDistanceArea [as 别名]
# 或者: from qgis.core.QgsDistanceArea import ellipsoid [as 别名]
    def testMeasureLineProjectedWorldPoints(self):
        #   +-+
        #   | |
        # +-+ +
        # checking returned length_mapunits/projected_points of diffferent world points with results from SpatiaLite ST_Project
        da_3068 = QgsDistanceArea()
        da_3068.setSourceCrs(QgsCoordinateReferenceSystem.fromOgcWmsCrs('EPSG:3068'), QgsProject.instance().transformContext())
        if (da_3068.sourceCrs().isGeographic()):
            da_3068.setEllipsoid(da_3068.sourceCrs().ellipsoidAcronym())
        self.assertEqual(da_3068.sourceCrs().authid(), 'EPSG:3068')
        print(("setting [{}] srid [{}] description [{}] isGeographic[{}] lengthUnits[{}] projectionAcronym[{}] ellipsoidAcronym[{}]".format(u'EPSG:3068', da_3068.sourceCrs().authid(), da_3068.sourceCrs().description(), da_3068.sourceCrs().isGeographic(), QgsUnitTypes.toString(da_3068.lengthUnits()), da_3068.sourceCrs().projectionAcronym(), da_3068.sourceCrs().ellipsoidAcronym())))
        da_wsg84 = QgsDistanceArea()
        da_wsg84.setSourceCrs(QgsCoordinateReferenceSystem.fromOgcWmsCrs('EPSG:4326'), QgsProject.instance().transformContext())
        if (da_wsg84.sourceCrs().isGeographic()):
            da_wsg84.setEllipsoid(da_wsg84.sourceCrs().ellipsoidAcronym())
        self.assertEqual(da_wsg84.sourceCrs().authid(), 'EPSG:4326')
        print(("setting [{}] srid [{}] description [{}] isGeographic[{}] lengthUnits[{}] projectionAcronym[{}] ellipsoidAcronym[{}] ellipsoid[{}]".format(u'EPSG:4326', da_wsg84.sourceCrs().authid(), da_wsg84.sourceCrs().description(), da_wsg84.sourceCrs().isGeographic(), QgsUnitTypes.toString(da_wsg84.lengthUnits()), da_wsg84.sourceCrs().projectionAcronym(), da_wsg84.sourceCrs().ellipsoidAcronym(), da_wsg84.ellipsoid())))
        da_4314 = QgsDistanceArea()
        da_4314.setSourceCrs(QgsCoordinateReferenceSystem.fromOgcWmsCrs('EPSG:4314'), QgsProject.instance().transformContext())
        if (da_4314.sourceCrs().isGeographic()):
            da_4314.setEllipsoid(da_4314.sourceCrs().ellipsoidAcronym())
        self.assertEqual(da_4314.sourceCrs().authid(), 'EPSG:4314')
        print(("setting [{}] srid [{}] description [{}] isGeographic[{}] lengthUnits[{}] projectionAcronym[{}] ellipsoidAcronym[{}]".format(u'EPSG:4314', da_4314.sourceCrs().authid(), da_4314.sourceCrs().description(), da_4314.sourceCrs().isGeographic(), QgsUnitTypes.toString(da_4314.lengthUnits()), da_4314.sourceCrs().projectionAcronym(), da_4314.sourceCrs().ellipsoidAcronym())))
        da_4805 = QgsDistanceArea()
        da_4805.setSourceCrs(QgsCoordinateReferenceSystem.fromOgcWmsCrs('EPSG:4805'), QgsProject.instance().transformContext())
        if (da_4805.sourceCrs().isGeographic()):
            da_4805.setEllipsoid(da_4805.sourceCrs().ellipsoidAcronym())
        self.assertEqual(da_4805.sourceCrs().authid(), 'EPSG:4805')
        print(("setting [{}] srid [{}] description [{}] isGeographic[{}] lengthUnits[{}] projectionAcronym[{}] ellipsoidAcronym[{}]".format(u'EPSG:4805', da_4805.sourceCrs().authid(), da_4805.sourceCrs().description(), da_4805.sourceCrs().isGeographic(), QgsUnitTypes.toString(da_4805.lengthUnits()), da_4805.sourceCrs().projectionAcronym(), da_4805.sourceCrs().ellipsoidAcronym())))
        # EPSG:5665 unknown, why?
        da_5665 = QgsDistanceArea()
        da_5665.setSourceCrs(QgsCoordinateReferenceSystem.fromOgcWmsCrs('EPSG:5665'), QgsProject.instance().transformContext())
        if (da_5665.sourceCrs().isGeographic()):
            da_5665.setEllipsoid(da_5665.sourceCrs().ellipsoidAcronym())
        print(("setting [{}] srid [{}] description [{}] isGeographic[{}] lengthUnits[{}] projectionAcronym[{}] ellipsoidAcronym[{}]".format(u'EPSG:5665', da_5665.sourceCrs().authid(), da_5665.sourceCrs().description(), da_5665.sourceCrs().isGeographic(), QgsUnitTypes.toString(da_5665.lengthUnits()), da_5665.sourceCrs().projectionAcronym(), da_5665.sourceCrs().ellipsoidAcronym())))
        #self.assertEqual(da_5665.sourceCrs().authid(), 'EPSG:5665')
        da_25833 = QgsDistanceArea()
        da_25833.setSourceCrs(QgsCoordinateReferenceSystem.fromOgcWmsCrs('EPSG:25833'), QgsProject.instance().transformContext())
        if (da_25833.sourceCrs().isGeographic()):
            da_25833.setEllipsoid(da_25833.sourceCrs().ellipsoidAcronym())
        print(("setting [{}] srid [{}] description [{}] isGeographic[{}] lengthUnits[{}] projectionAcronym[{}] ellipsoidAcronym[{}]".format(u'EPSG:25833', da_25833.sourceCrs().authid(), da_25833.sourceCrs().description(), da_25833.sourceCrs().isGeographic(), QgsUnitTypes.toString(da_25833.lengthUnits()), da_25833.sourceCrs().projectionAcronym(), da_25833.sourceCrs().ellipsoidAcronym())))
        self.assertEqual(da_25833.sourceCrs().authid(), 'EPSG:25833')

        # Berlin - Brandenburg Gate - Quadriga
        point_berlin_3068 = QgsPointXY(23183.38449999984, 21047.3225000017)
        point_berlin_3068_project = point_berlin_3068.project(1, (math.pi / 2))
        point_meter_result = QgsPointXY(0, 0)
        length_meter_mapunits, point_meter_result = da_3068.measureLineProjected(point_berlin_3068, 1.0, (math.pi / 2))
        pprint(point_meter_result)
        print('-I-> Berlin 3068 length_meter_mapunits[{}] point_meter_result[{}]'.format(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, da_3068.lengthUnits(), True), point_meter_result.asWkt()))
        self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 1, da_3068.lengthUnits(), True), '1.0 m')
        self.assertEqual(point_meter_result.toString(7), point_berlin_3068_project.toString(7))
        point_berlin_wsg84 = QgsPointXY(13.37770458660236, 52.51627178856762)
        point_berlin_wsg84_project = QgsPointXY(13.37771931736259, 52.51627178856669)
        length_meter_mapunits, point_meter_result = da_wsg84.measureLineProjected(point_berlin_wsg84, 1.0, (math.pi / 2))
        print('-I-> Berlin Wsg84 length_meter_mapunits[{}] point_meter_result[{}] ellipsoid[{}]'.format(QgsDistanceArea.formatDistance(length_meter_mapunits, 20, da_wsg84.lengthUnits(), True), point_meter_result.asWkt(), da_wsg84.ellipsoid()))
        # for unknown reasons, this is returning '0.00001473026 m' instead of '0.00001473026 deg' when using da_wsg84.lengthUnits()
        # self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits,11,da_wsg84.lengthUnits(),True), '0.00001473026 deg')
        self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 11, QgsUnitTypes.DistanceDegrees, True), '0.00001473026 deg')
        self.assertEqual(point_meter_result.toString(7), point_berlin_wsg84_project.toString(7))
        point_berlin_4314 = QgsPointXY(13.37944343021465, 52.51767872437083)
        point_berlin_4314_project = QgsPointXY(13.37945816324759, 52.5176787243699)
        length_meter_mapunits, point_meter_result = da_4314.measureLineProjected(point_berlin_4314, 1.0, (math.pi / 2))
        print('-I-> Berlin 4314 length_meter_mapunits[{}] point_meter_result[{}]'.format(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, da_4314.lengthUnits(), True), point_meter_result.asWkt()))
        self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 9, QgsUnitTypes.DistanceDegrees, True), '0.000014733 deg')
        self.assertEqual(point_meter_result.toString(7), point_berlin_4314_project.toString(7))
        point_berlin_4805 = QgsPointXY(31.04960570069176, 52.5174657497405)
        point_berlin_4805_project = QgsPointXY(31.04962043365347, 52.51746574973957)
        length_meter_mapunits, point_meter_result = da_4805.measureLineProjected(point_berlin_4805, 1.0, (math.pi / 2))
        print('-I-> Berlin 4805 length_meter_mapunits[{}] point_meter_result[{}]'.format(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, da_4805.lengthUnits(), True), point_meter_result.asWkt()))
        self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 9, QgsUnitTypes.DistanceDegrees, True), '0.000014733 deg')
        self.assertEqual(point_meter_result.toString(7), point_berlin_4805_project.toString(7))
        point_berlin_25833 = QgsPointXY(389918.0748318382, 5819698.772194743)
        point_berlin_25833_project = point_berlin_25833.project(1, (math.pi / 2))
        length_meter_mapunits, point_meter_result = da_25833.measureLineProjected(point_berlin_25833, 1.0, (math.pi / 2))
        print('-I-> Berlin 25833 length_meter_mapunits[{}] point_meter_result[{}]'.format(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, da_25833.lengthUnits(), True), point_meter_result.asWkt()))
        self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, da_25833.lengthUnits(), True), '1.0000000 m')
        self.assertEqual(point_meter_result.toString(7), point_berlin_25833_project.toString(7))
        if da_5665.sourceCrs().authid() != "":
            point_berlin_5665 = QgsPointXY(3389996.871728864, 5822169.719727578)
            point_berlin_5665_project = point_berlin_5665.project(1, (math.pi / 2))
            length_meter_mapunits, point_meter_result = da_5665.measureLineProjected(point_berlin_5665, 1.0, (math.pi / 2))
            print('-I-> Berlin 5665 length_meter_mapunits[{}] point_meter_result[{}]'.format(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, da_5665.lengthUnits(), True), point_meter_result.asWkt()))
            self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 1.0, da_5665.lengthUnits(), True), '1.0 m')
            self.assertEqual(point_meter_result.toString(7), point_berlin_5665_project.toString(7))
        print('\n12 points ''above over'' and on the Equator')
        point_wsg84 = QgsPointXY(25.7844, 71.1725)
        point_wsg84_project = QgsPointXY(25.78442775215388, 71.17249999999795)
        length_meter_mapunits, point_meter_result = da_wsg84.measureLineProjected(point_wsg84, 1.0, (math.pi / 2))
        print('-I-> Nordkap, Norway - Wsg84 - length_meter_mapunits[{}] point_meter_result[{}] '.format(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, da_wsg84.lengthUnits(), True), point_meter_result.asWkt()))
        self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, QgsUnitTypes.DistanceDegrees, True), '0.0000278 deg')
        self.assertEqual(point_meter_result.toString(7), point_wsg84_project.toString(7))
        point_wsg84 = QgsPointXY(24.95995, 60.16841)
        point_wsg84_project = QgsPointXY(24.95996801277454, 60.16840999999877)
        length_meter_mapunits, point_meter_result = da_wsg84.measureLineProjected(point_wsg84, 1.0, (math.pi / 2))
        print('-I-> Helsinki, Finnland - Wsg84 - length_meter_mapunits[{}] point_meter_result[{}] '.format(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, da_wsg84.lengthUnits(), True), point_meter_result.asWkt()))
        self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 8, QgsUnitTypes.DistanceDegrees, True), '0.00001801 deg')
        self.assertEqual(point_meter_result.toString(7), point_wsg84_project.toString(7))
        point_wsg84 = QgsPointXY(12.599278, 55.692861)
        point_wsg84_project = QgsPointXY(12.59929390161872, 55.69286099999897)
#.........这里部分代码省略.........

# 需要导入模块: from qgis.core import QgsDistanceArea [as 别名]
# 或者: from qgis.core.QgsDistanceArea import ellipsoid [as 别名]
    def by_seconds(self, seconds):
        """move by variable distance"""

        def iterations(distance,h):
            FIe1=0
            LAMe1=0
            FI=100
            LAM=100

            # iterations
            while fabs(FI-FIe1)>0.0000000001 and fabs(LAM -LAMe1)>0.0000000001:
                FI=FIe1
                LAM=LAMe1
                for i in range(0,int(ceil(distance/h))):
                    kfi=[]
                    klam=[]
                    kazi=[]
                    kfi.append(cos(azi[i])/(a*(1-(e2))/(pow((sqrt(1-(e2)*pow(sin(fi[i]),2))),3))))
                    klam.append(sin(azi[i])/((a/(sqrt(1-(e2)*pow(sin(fi[i]),2))))*cos(fi[i])))
                    kazi.append(sin(azi[i])*tan(fi[i])/(a/(sqrt(1-(e2)*pow(sin(fi[i]),2)))))
                    for j in range(1,3):
                        kfi.append(cos(azi[i]+kazi[j-1]*h/2)/(a*(1-(e2))/(pow(sqrt(1-(e2)*pow(sin(fi[i]+kfi[j-1]*h/2),2)),3))))
                        klam.append(sin(azi[i]+kazi[j-1]*h/2)/((a/(sqrt(1-(e2)*pow(sin(fi[i]+kfi[j-1]*h/2),2))))*cos(fi[i]+kfi[j-1]*h/2)))
                        kazi.append(sin(azi[i]+kazi[j-1]*h/2)*tan(fi[i]+kfi[j-1]*h/2)/(a/(sqrt(1-(e2)*pow(sin(fi[i]+kfi[j-1]*h/2),2)))))

                    kfi.append(cos(azi[i]+kazi[2]*h)/(a*(1-(e2))/(pow(sqrt(1-(e2)*pow(sin(fi[i]+kfi[2]*h),2)),3))))
                    klam.append(sin(azi[i]+kazi[2]*h)/((a/(sqrt(1-(e2)*pow(sin(fi[i]+kfi[2]*h),2))))*cos(fi[i]+kfi[2]*h)))
                    kazi.append(sin(azi[i]+kazi[2]*h)*tan(fi[i]+kfi[2]*h)/(a/(sqrt(1-(e2)*pow(sin(fi[i]+kfi[2]*h),2)))))

                    fi.append(fi[i]+(h/6.0)*(kfi[0]+2*kfi[1]+2*kfi[2]+kfi[3]))
                    lam.append(lam[i]+(h/6.0)*(klam[0]+2*klam[1]+2*klam[2]+klam[3]))
                    azi.append(azi[i]+(h/6.0)*(kazi[0]+2*kazi[1]+2*kazi[2]+kazi[3]))

                FIe1=fi[i+1]
                LAMe1=lam[i+1]
                h=h/2
                fi[1:]=[]
                lam[1:]=[]
                azi[1:]=[]

            return FIe1,LAMe1

        self._check()
        header=self.inputfile.readline()
        beforeLat=header.split('Lat_deg')
        numberOfLatColumn=beforeLat[0].split(',')
        beforeLong=header.split('Lon_deg')
        numberOfLonColumn=beforeLong[0].split(',')
        beforeSec=header.split('Gtm_sec')
        numberOfSecColumn=beforeSec[0].split(',')
        self.outputfile.write(header)

        d = QgsDistanceArea()
        d.setEllipsoid('WGS84')
        d.setEllipsoidalMode(True)
        d.ellipsoid()
        a = 6378137.0 # WGS84 ellipsoid parametres
        e2 = 0.081819190842622
        line1=self.inputfile.readline()

        if seconds>0:
            linePos=self.inputfile.tell()
            line1=line1.split(',')
            while line1:
                self.inputfile.seek(linePos)
                line2=self.inputfile.readline()
                linePos=self.inputfile.tell()
                moveTime=1*seconds
                outline=1*line1
                inline=line2.split(',')

                while moveTime>0: # for case of more than 1 second
                    if line2:
                        line2=line2.split(',')
                        p1=QgsPoint(float(line1[len(numberOfLonColumn)-1]),float(line1[len(numberOfLatColumn)-1]))
                        p2=QgsPoint(float(line2[len(numberOfLonColumn)-1]),float(line2[len(numberOfLatColumn)-1]))

                        if p1!=p2:
                            aziA = d.bearing(p1,p2)
                            l = d.computeDistanceBearing(p1,p2)[0]

                            if moveTime>(float(line2[len(numberOfSecColumn)-1])-float(line1[len(numberOfSecColumn)-1])):
                                moveTime=moveTime-(float(line2[len(numberOfSecColumn)-1])-float(line1[len(numberOfSecColumn)-1]))
                            elif moveTime!=0 and moveTime!=(float(line2[len(numberOfSecColumn)-1])-float(line1[len(numberOfSecColumn)-1])): #first geodetic problem
                                distance=l/(float(line2[len(numberOfSecColumn)-1])-float(line1[len(numberOfSecColumn)-1]))*moveTime
                                h=distance/2.0
                                fi=[float(line1[len(numberOfLatColumn)-1])*pi/180]
                                lam=[float(line1[len(numberOfLonColumn)-1])*pi/180]
                                azi=[aziA]

                                FIe1,LAMe1 = iterations(distance,h)
                                moveTime=0

                            else:
                                FIe1=float(line2[len(numberOfLatColumn)-1])*pi/180
                                LAMe1=float(line2[len(numberOfLonColumn)-1])*pi/180
                                moveTime=0
                                break
                        else:
                            if moveTime>(float(line2[len(numberOfSecColumn)-1])-float(line1[len(numberOfSecColumn)-1])):
#.........这里部分代码省略.........