def testSubstitutionMap(self):
        """Test that we can use degree symbols in substitutions.
        """
        # Create a point and convert it to text containing a degree symbol.
        myPoint = QgsPointXY(12.3, -33.33)
        myCoordinates = myPoint.toDegreesMinutesSeconds(2)
        myTokens = myCoordinates.split(',')
        myLongitude = myTokens[0]
        myLatitude = myTokens[1]
        myText = 'Latitude: %s, Longitude: %s' % (myLatitude, myLongitude)

        # Load the composition with the substitutions
        myComposition = QgsComposition(QgsProject.instance())
        mySubstitutionMap = {'replace-me': myText}
        myFile = os.path.join(TEST_DATA_DIR, 'template-for-substitution.qpt')
        with open(myFile) as f:
            myTemplateContent = f.read()
        myDocument = QDomDocument()
        myDocument.setContent(myTemplateContent)
        myComposition.loadFromTemplate(myDocument, mySubstitutionMap)

        # We should be able to get map0
        myMap = myComposition.getComposerMapById(0)
        myMessage = ('Map 0 could not be found in template %s', myFile)
        assert myMap is not None, myMessage

    def testRenderMetersInMapUnits(self):

        crs_wsg84 = QgsCoordinateReferenceSystem.fromOgcWmsCrs('EPSG:4326')
        rt_extent = QgsRectangle(13.37768985634235, 52.51625705830762, 13.37771931686235, 52.51628651882762)
        point_berlin_wsg84 = QgsPointXY(13.37770458660236, 52.51627178856762)
        length_wsg84_mapunits = 0.00001473026350140572
        meters_test = 2.40
        da_wsg84 = QgsDistanceArea()
        da_wsg84.setSourceCrs(crs_wsg84, QgsProject.instance().transformContext())
        if (da_wsg84.sourceCrs().isGeographic()):
            da_wsg84.setEllipsoid(da_wsg84.sourceCrs().ellipsoidAcronym())
        length_meter_mapunits = da_wsg84.measureLineProjected(point_berlin_wsg84, 1.0, (math.pi / 2))
        meters_test_mapunits = meters_test * length_wsg84_mapunits
        meters_test_pixel = meters_test * length_wsg84_mapunits
        ms = QgsMapSettings()
        ms.setDestinationCrs(crs_wsg84)
        ms.setExtent(rt_extent)
        r = QgsRenderContext.fromMapSettings(ms)
        r.setExtent(rt_extent)
        self.assertEqual(r.extent().center().toString(7), point_berlin_wsg84.toString(7))
        c = QgsMapUnitScale()
        r.setDistanceArea(da_wsg84)
        result_test_painterunits = r.convertToPainterUnits(meters_test, QgsUnitTypes.RenderMetersInMapUnits, c)
        self.assertEqual(QgsDistanceArea.formatDistance(result_test_painterunits, 7, QgsUnitTypes.DistanceUnknownUnit, True), QgsDistanceArea.formatDistance(meters_test_mapunits, 7, QgsUnitTypes.DistanceUnknownUnit, True))
        result_test_mapunits = r.convertToMapUnits(meters_test, QgsUnitTypes.RenderMetersInMapUnits, c)
        self.assertEqual(QgsDistanceArea.formatDistance(result_test_mapunits, 7, QgsUnitTypes.DistanceDegrees, True), QgsDistanceArea.formatDistance(meters_test_mapunits, 7, QgsUnitTypes.DistanceDegrees, True))
        result_test_meters = r.convertFromMapUnits(meters_test_mapunits, QgsUnitTypes.RenderMetersInMapUnits)
        self.assertEqual(QgsDistanceArea.formatDistance(result_test_meters, 1, QgsUnitTypes.DistanceMeters, True), QgsDistanceArea.formatDistance(meters_test, 1, QgsUnitTypes.DistanceMeters, True))

class TestQgsPointXY(unittest.TestCase):

    def __init__(self, methodName):
        """Run once on class initialization."""
        unittest.TestCase.__init__(self, methodName)

    def setUp(self):
        self.mPoint = QgsPointXY(10.0, 10.0)

    def test_Point(self):
        myExpectedValue = 10.0
        myActualValue = self.mPoint.x()
        myMessage = 'Expected: %s Got: %s' % (myExpectedValue, myActualValue)
        assert myExpectedValue == myActualValue, myMessage

    def test_pointToString(self):
        myExpectedValue = '10, 10'
        myActualValue = self.mPoint.toString()
        myMessage = 'Expected: %s Got: %s' % (myExpectedValue, myActualValue)
        assert myExpectedValue == myActualValue, myMessage

    def test_hash(self):
        a = QgsPointXY(2.0, 1.0)
        b = QgsPointXY(2.0, 2.0)
        c = QgsPointXY(1.0, 2.0)
        d = QgsPointXY(1.0, 1.0)
        e = QgsPointXY(2.0, 1.0)
        assert a.__hash__() != b.__hash__()
        assert e.__hash__() == a.__hash__()

        mySet = set([a, b, c, d, e])
        assert len(mySet) == 4

    def __init__(self, plugin, filepath, title, screenExtent):
        QgsPluginLayer.__init__(
            self, FreehandRasterGeoreferencerLayer.LAYER_TYPE, title)
        self.plugin = plugin
        self.iface = plugin.iface

        self.title = title
        self.filepath = filepath
        self.screenExtent = screenExtent
        self.history = []
        # set custom properties
        self.setCustomProperty("title", title)
        self.setCustomProperty("filepath", self.filepath)

        self.setValid(True)

        self.setTransparency(LayerDefaultSettings.TRANSPARENCY)
        self.setBlendModeByName(LayerDefaultSettings.BLEND_MODE)

        # dummy data: real init is done in intializeLayer
        self.center = QgsPointXY(0, 0)
        self.rotation = 0.0
        self.xScale = 1.0
        self.yScale = 1.0

        self.error = False
        self.initializing = False
        self.initialized = False
        self.initializeLayer(screenExtent)
        self._extent = None

        self.provider = FreehandRasterGeoreferencerLayerProvider(self)

    def testPoint(self):
        point = QgsReferencedPointXY(QgsPointXY(1.0, 2.0), QgsCoordinateReferenceSystem('epsg:3111'))
        self.assertEqual(point.x(), 1.0)
        self.assertEqual(point.y(), 2.0)
        self.assertEqual(point.crs().authid(), 'EPSG:3111')

        point.setCrs(QgsCoordinateReferenceSystem('epsg:28356'))
        self.assertEqual(point.crs().authid(), 'EPSG:28356')

        # in variant
        v = QVariant(QgsReferencedPointXY(QgsPointXY(3.0, 4.0), QgsCoordinateReferenceSystem('epsg:3111')))
        self.assertEqual(v.value().x(), 3.0)
        self.assertEqual(v.value().y(), 4.0)
        self.assertEqual(v.value().crs().authid(), 'EPSG:3111')

        # to QgsPointXY
        p = QgsPointXY(point)
        self.assertEqual(p.x(), 1.0)
        self.assertEqual(p.y(), 2.0)

    def test_hash(self):
        a = QgsPointXY(2.0, 1.0)
        b = QgsPointXY(2.0, 2.0)
        c = QgsPointXY(1.0, 2.0)
        d = QgsPointXY(1.0, 1.0)
        e = QgsPointXY(2.0, 1.0)
        assert a.__hash__() != b.__hash__()
        assert e.__hash__() == a.__hash__()

        mySet = set([a, b, c, d, e])
        assert len(mySet) == 4

    def checkPoint(self, point):
        state = 'yellow'

        if type(point) != QgsPointXY:
            try:
                if type(point) == list and len(point) == 2:
                    point = QgsPointXY(point[0], point[1])
                else:
                    point = QgsPointXY(point)
            except TypeError:
                return state, None
            
        if self.dhm != {}:
            extent = self.dhm['extent']
            [extLx, extHy, extHx, extLy] = extent
            
            if extLx <= float(point.x()) <= extHx \
                    and extLy <= float(point.y()) <= extHy:
                state = 'green'
            else:
                state = 'red'

        return state, point

 def readXml(self, node, context):
     self.readCustomProperties(node)
     self.title = self.customProperty("title", "")
     self.filepath = self.customProperty("filepath", "")
     self.xScale = float(self.customProperty("xScale", 1.0))
     self.yScale = float(self.customProperty("yScale", 1.0))
     self.rotation = float(self.customProperty("rotation", 0.0))
     xCenter = float(self.customProperty("xCenter", 0.0))
     yCenter = float(self.customProperty("yCenter", 0.0))
     self.center = QgsPointXY(xCenter, yCenter)
     self.setTransparency(int(self.customProperty(
         "transparency", LayerDefaultSettings.TRANSPARENCY)))
     self.setBlendModeByName(self.customProperty(
         "blendMode", LayerDefaultSettings.BLEND_MODE))
     return True

def run(sources_layer_path, receivers_layer_path, emission_pts_layer_path, research_ray):
    sources_layer = QgsVectorLayer(sources_layer_path, "input layer", "ogr")
    receivers_layer = QgsVectorLayer(receivers_layer_path, "output layer", "ogr")

    sources_feat_all = sources_layer.dataProvider().getFeatures()

    receivers_feat_all_dict = {}
    receivers_feat_all = receivers_layer.dataProvider().getFeatures()
    receivers_spIndex = QgsSpatialIndex()
    for receivers_feat in receivers_feat_all:
        receivers_spIndex.insertFeature(receivers_feat)
        receivers_feat_all_dict[receivers_feat.id()] = receivers_feat

    emission_pts_fields = QgsFields()
    emission_pts_fields.append(QgsField("id_emi", QVariant.Int))
    emission_pts_fields.append(QgsField("id_emi_source", QVariant.Int))
    emission_pts_fields.append(QgsField("id_source", QVariant.Int))
    emission_pts_fields.append(QgsField("d_rTOe", QVariant.Double, len=10, prec=2))
    # update for QGIS 3 converting VectorWriter to QgsVectorFileWriter
    # emission_pts_writer = VectorWriter(emission_pts_layer_path, None, emission_pts_fields, 0, sources_layer.crs())

    emission_pts_writer = QgsVectorFileWriter(emission_pts_layer_path, "System",
                                              emission_pts_fields, QgsWkbTypes.Point, sources_layer.crs(),
                                              "ESRI Shapefile")

    # initializes ray and emission point id
    emission_pt_id = 0

    for sources_feat in sources_feat_all:

        # researches the receiver points in a rectangle created by the research_ray
        # creates the search rectangle
        rect = QgsRectangle()
        rect.setXMinimum(sources_feat.geometry().boundingBox().xMinimum() - research_ray)
        rect.setXMaximum(sources_feat.geometry().boundingBox().xMaximum() + research_ray)
        rect.setYMinimum(sources_feat.geometry().boundingBox().yMinimum() - research_ray)
        rect.setYMaximum(sources_feat.geometry().boundingBox().yMaximum() + research_ray)

        receiver_pts_request = receivers_spIndex.intersects(rect)

        distance_min = []
        for receiver_pts_id in receiver_pts_request:
            receiver_pts_feat = receivers_feat_all_dict[receiver_pts_id]
            result = sources_feat.geometry().closestSegmentWithContext(receiver_pts_feat.geometry().asPoint())
            distance_min_tmp = sqrt(result[0])

            if distance_min_tmp <= research_ray:
                distance_min.append(distance_min_tmp)

        # defines segment max length
        if len(distance_min) >= 1:
            segment_max = min(distance_min) / 2
            if segment_max < 2:
                segment_max = 2
        else:
            continue

        # splits the sources line in emission points at a fix distance (minimum distance/2) and create the emission point layer
        # gets vertex
        sources_geom = sources_feat.geometry()
        if sources_geom.isMultipart():
            sources_geom.convertToSingleType()
        sources_feat_vertex_pt_all = sources_geom.asPolyline()

        emission_pt_id_road = 0

        for i in range(0, len(sources_feat_vertex_pt_all)):

            pt1 = QgsPointXY(sources_feat_vertex_pt_all[i])

            add_point_to_layer(emission_pts_writer, pt1,
                               [emission_pt_id, emission_pt_id_road, sources_feat.id(), segment_max])

            emission_pt_id = emission_pt_id + 1
            emission_pt_id_road = emission_pt_id_road + 1

            if i < len(sources_feat_vertex_pt_all) - 1:

                pt2 = QgsPoint(sources_feat_vertex_pt_all[i + 1])

                x1 = pt1.x()
                y1 = pt1.y()
                x2 = pt2.x()
                y2 = pt2.y()

                if y2 == y1:
                    dx = segment_max
                    dy = 0
                    m = 0
                elif x2 == x1:
                    dx = 0
                    dy = segment_max
                else:
                    m = (y2 - y1) / (x2 - x1)
                    dx = sqrt((segment_max ** 2) / (1 + m ** 2))
                    dy = sqrt(((segment_max ** 2) * (m ** 2)) / (1 + m ** 2))

                pt = pt1

                while compute_distance(pt, pt2) > segment_max:
#.........这里部分代码省略.........

 def setUp(self):
     self.mPoint = QgsPointXY(10.0, 10.0)

    def processAlgorithm(self, parameters, context, feedback):
        layer = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.INPUT_VECTOR), context)
        startPoints = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.START_POINTS), context)
        endPoint = self.getParameterValue(self.END_POINT)
        strategy = self.getParameterValue(self.STRATEGY)

        directionFieldName = self.getParameterValue(self.DIRECTION_FIELD)
        forwardValue = self.getParameterValue(self.VALUE_FORWARD)
        backwardValue = self.getParameterValue(self.VALUE_BACKWARD)
        bothValue = self.getParameterValue(self.VALUE_BOTH)
        defaultDirection = self.getParameterValue(self.DEFAULT_DIRECTION)
        bothValue = self.getParameterValue(self.VALUE_BOTH)
        defaultDirection = self.getParameterValue(self.DEFAULT_DIRECTION)
        speedFieldName = self.getParameterValue(self.SPEED_FIELD)
        defaultSpeed = self.getParameterValue(self.DEFAULT_SPEED)
        tolerance = self.getParameterValue(self.TOLERANCE)

        fields = QgsFields()
        fields.append(QgsField('start', QVariant.String, '', 254, 0))
        fields.append(QgsField('end', QVariant.String, '', 254, 0))
        fields.append(QgsField('cost', QVariant.Double, '', 20, 7))

        feat = QgsFeature()
        feat.setFields(fields)

        writer = self.getOutputFromName(
            self.OUTPUT_LAYER).getVectorWriter(fields, QgsWkbTypes.LineString, layer.crs(), context)

        tmp = endPoint.split(',')
        endPoint = QgsPointXY(float(tmp[0]), float(tmp[1]))

        directionField = -1
        if directionFieldName is not None:
            directionField = layer.fields().lookupField(directionFieldName)
        speedField = -1
        if speedFieldName is not None:
            speedField = layer.fields().lookupField(speedFieldName)

        director = QgsVectorLayerDirector(layer,
                                          directionField,
                                          forwardValue,
                                          backwardValue,
                                          bothValue,
                                          defaultDirection)

        distUnit = iface.mapCanvas().mapSettings().destinationCrs().mapUnits()
        multiplier = QgsUnitTypes.fromUnitToUnitFactor(distUnit, QgsUnitTypes.DistanceMeters)
        if strategy == 0:
            strategy = QgsNetworkDistanceStrategy()
        else:
            strategy = QgsNetworkSpeedStrategy(speedField,
                                               defaultSpeed,
                                               multiplier * 1000.0 / 3600.0)
            multiplier = 3600

        director.addStrategy(strategy)
        builder = QgsGraphBuilder(iface.mapCanvas().mapSettings().destinationCrs(),
                                  True,
                                  tolerance)

        feedback.pushInfo(self.tr('Loading start points...'))
        request = QgsFeatureRequest()
        request.setFlags(request.flags() ^ QgsFeatureRequest.SubsetOfAttributes)
        features = QgsProcessingUtils.getFeatures(startPoints, context, request)
        count = QgsProcessingUtils.featureCount(startPoints, context)

        points = [endPoint]
        for f in features:
            points.append(f.geometry().asPoint())

        feedback.pushInfo(self.tr('Building graph...'))
        snappedPoints = director.makeGraph(builder, points)

        feedback.pushInfo(self.tr('Calculating shortest paths...'))
        graph = builder.graph()

        idxEnd = graph.findVertex(snappedPoints[0])
        route = []

        total = 100.0 / count if count else 1
        for i in range(1, count + 1):
            idxStart = graph.findVertex(snappedPoints[i])
            tree, cost = QgsGraphAnalyzer.dijkstra(graph, idxStart, 0)

            if tree[idxEnd] == -1:
                msg = self.tr('There is no route from start point ({}) to end point ({}).'.format(points[i].toString(), endPoint.toString()))
                feedback.setProgressText(msg)
                QgsMessageLog.logMessage(msg, self.tr('Processing'), QgsMessageLog.WARNING)
                continue

            cost = 0.0
            current = idxEnd
            while current != idxStart:
                cost += graph.edge(tree[current]).cost(0)
                route.append(graph.vertex(graph.edge(tree[current]).inVertex()).point())
                current = graph.edge(tree[current]).outVertex()

            route.append(snappedPoints[i])
            route.reverse()

#.........这里部分代码省略.........

    def balanced(features, graph, feedback, balance=0, min_colors=4):
        feature_colors = {}
        # start with minimum number of colors in pool
        color_pool = set(range(1, min_colors + 1))

        # calculate count of neighbours
        neighbour_count = defaultdict(int)
        for feature_id, neighbours in graph.node_edge.items():
            neighbour_count[feature_id] += len(neighbours)

        # sort features by neighbour count - we want to handle those with more neighbours first
        sorted_by_count = [feature_id for feature_id in sorted(neighbour_count.items(),
                                                               key=operator.itemgetter(1),
                                                               reverse=True)]
        # counts for each color already assigned
        color_counts = defaultdict(int)
        color_areas = defaultdict(float)
        for c in color_pool:
            color_counts[c] = 0
            color_areas[c] = 0

        total = 10.0 / len(sorted_by_count)
        i = 0

        for (feature_id, n) in sorted_by_count:
            if feedback.isCanceled():
                break

            # first work out which already assigned colors are adjacent to this feature
            adjacent_colors = set()
            for neighbour in graph.node_edge[feature_id]:
                if neighbour in feature_colors:
                    adjacent_colors.add(feature_colors[neighbour])

            # from the existing colors, work out which are available (ie non-adjacent)
            available_colors = color_pool.difference(adjacent_colors)

            feature_color = -1
            if len(available_colors) == 0:
                # no existing colors available for this feature, so add new color to pool and repeat
                min_colors += 1
                return ColoringAlgorithm.balanced(features, graph, feedback, balance, min_colors)
            else:
                if balance == 0:
                    # choose least used available color
                    counts = [(c, v) for c, v in color_counts.items() if c in available_colors]
                    feature_color = sorted(counts, key=operator.itemgetter(1))[0][0]
                    color_counts[feature_color] += 1
                elif balance == 1:
                    areas = [(c, v) for c, v in color_areas.items() if c in available_colors]
                    feature_color = sorted(areas, key=operator.itemgetter(1))[0][0]
                    color_areas[feature_color] += features[feature_id].geometry().area()
                elif balance == 2:
                    min_distances = {c: sys.float_info.max for c in available_colors}
                    this_feature_centroid = QgsPointXY(features[feature_id].geometry().centroid().geometry())

                    # find features for all available colors
                    other_features = {f_id: c for (f_id, c) in feature_colors.items() if c in available_colors}

                    # loop through these, and calculate the minimum distance from this feature to the nearest
                    # feature with each assigned color
                    for other_feature_id, c in other_features.items():
                        if feedback.isCanceled():
                            break

                        other_geometry = features[other_feature_id].geometry()
                        other_centroid = QgsPointXY(other_geometry.centroid().geometry())

                        distance = this_feature_centroid.distanceSquared(other_centroid)
                        if distance < min_distances[c]:
                            min_distances[c] = distance

                    # choose color such that minimum distance is maximised! ie we want MAXIMAL separation between
                    # features with the same color
                    feature_color = sorted(min_distances, key=min_distances.__getitem__, reverse=True)[0]

            feature_colors[feature_id] = feature_color

            i += 1
            feedback.setProgress(70 + int(i * total))

        return feature_colors

    def processAlgorithm(self, parameters, context, feedback):
        layer = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.INPUT_VECTOR), context)
        startPoint = self.getParameterValue(self.START_POINT)
        endPoint = self.getParameterValue(self.END_POINT)
        strategy = self.getParameterValue(self.STRATEGY)

        directionFieldName = self.getParameterValue(self.DIRECTION_FIELD)
        forwardValue = self.getParameterValue(self.VALUE_FORWARD)
        backwardValue = self.getParameterValue(self.VALUE_BACKWARD)
        bothValue = self.getParameterValue(self.VALUE_BOTH)
        defaultDirection = self.getParameterValue(self.DEFAULT_DIRECTION)
        bothValue = self.getParameterValue(self.VALUE_BOTH)
        defaultDirection = self.getParameterValue(self.DEFAULT_DIRECTION)
        speedFieldName = self.getParameterValue(self.SPEED_FIELD)
        defaultSpeed = self.getParameterValue(self.DEFAULT_SPEED)
        tolerance = self.getParameterValue(self.TOLERANCE)

        fields = QgsFields()
        fields.append(QgsField('start', QVariant.String, '', 254, 0))
        fields.append(QgsField('end', QVariant.String, '', 254, 0))
        fields.append(QgsField('cost', QVariant.Double, '', 20, 7))

        writer = self.getOutputFromName(
            self.OUTPUT_LAYER).getVectorWriter(fields, QgsWkbTypes.LineString, layer.crs(), context)

        tmp = startPoint.split(',')
        startPoint = QgsPointXY(float(tmp[0]), float(tmp[1]))
        tmp = endPoint.split(',')
        endPoint = QgsPointXY(float(tmp[0]), float(tmp[1]))

        directionField = -1
        if directionFieldName is not None:
            directionField = layer.fields().lookupField(directionFieldName)
        speedField = -1
        if speedFieldName is not None:
            speedField = layer.fields().lookupField(speedFieldName)

        director = QgsVectorLayerDirector(layer,
                                          directionField,
                                          forwardValue,
                                          backwardValue,
                                          bothValue,
                                          defaultDirection)

        distUnit = iface.mapCanvas().mapSettings().destinationCrs().mapUnits()
        multiplier = QgsUnitTypes.fromUnitToUnitFactor(distUnit, QgsUnitTypes.DistanceMeters)
        if strategy == 0:
            strategy = QgsNetworkDistanceStrategy()
        else:
            strategy = QgsNetworkSpeedStrategy(speedField,
                                               defaultSpeed,
                                               multiplier * 1000.0 / 3600.0)
            multiplier = 3600

        director.addStrategy(strategy)
        builder = QgsGraphBuilder(iface.mapCanvas().mapSettings().destinationCrs(),
                                  True,
                                  tolerance)
        feedback.pushInfo(self.tr('Building graph...'))
        snappedPoints = director.makeGraph(builder, [startPoint, endPoint])

        feedback.pushInfo(self.tr('Calculating shortest path...'))
        graph = builder.graph()
        idxStart = graph.findVertex(snappedPoints[0])
        idxEnd = graph.findVertex(snappedPoints[1])

        tree, cost = QgsGraphAnalyzer.dijkstra(graph, idxStart, 0)
        if tree[idxEnd] == -1:
            raise GeoAlgorithmExecutionException(
                self.tr('There is no route from start point to end point.'))

        route = []
        cost = 0.0
        current = idxEnd
        while current != idxStart:
            cost += graph.edge(tree[current]).cost(0)
            route.append(graph.vertex(graph.edge(tree[current]).inVertex()).point())
            current = graph.edge(tree[current]).outVertex()

        route.append(snappedPoints[0])
        route.reverse()

        self.setOutputValue(self.TRAVEL_COST, cost / multiplier)

        feedback.pushInfo(self.tr('Writing results...'))
        geom = QgsGeometry.fromPolyline(route)
        feat = QgsFeature()
        feat.setFields(fields)
        feat['start'] = startPoint.toString()
        feat['end'] = endPoint.toString()
        feat['cost'] = cost / multiplier
        feat.setGeometry(geom)
        writer.addFeature(feat)
        del writer

    def testMeasureLineProjected(self):
        #   +-+
        #   | |
        # +-+ +
        # test setting/getting the source CRS
        da_3068 = QgsDistanceArea()
        da_wsg84 = QgsDistanceArea()

        da_3068.setSourceCrs(QgsCoordinateReferenceSystem.fromOgcWmsCrs('EPSG:3068'), QgsProject.instance().transformContext())
        if (da_3068.sourceCrs().isGeographic()):
            da_3068.setEllipsoid(da_3068.sourceCrs().ellipsoidAcronym())
        print(("setting [{}] srid [{}] description [{}]".format(u'Soldner Berlin', da_3068.sourceCrs().authid(), da_3068.sourceCrs().description())))
        self.assertEqual(da_3068.sourceCrs().authid(), 'EPSG:3068')
        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[{}]".format(u'Wsg84', da_wsg84.sourceCrs().authid(), da_wsg84.sourceCrs().description(), da_wsg84.sourceCrs().isGeographic())))
        # print(("-- projectionAcronym[{}] ellipsoidAcronym[{}] toWkt[{}] mapUnits[{}] toProj4[{}]".format(da_wsg84.sourceCrs().projectionAcronym(),da_wsg84.sourceCrs().ellipsoidAcronym(), da_wsg84.sourceCrs().toWkt(),da_wsg84.sourceCrs().mapUnits(),da_wsg84.sourceCrs().toProj4())))
        print(("Testing Position change for[{}] years[{}]".format(u'Ampelanlage - Potsdamer Platz, Verkehrsinsel', u'1924 and 1998')))

        # 1924-10-24 SRID=3068;POINT(23099.49 20296.69)
        # 1924-10-24 SRID=4326;POINT(13.37650707988041 52.50952361017194)
        # 1998-10-02 SRID=3068;POINT(23082.30 20267.80)
        # 1998-10-02 SRID=4326;POINT(13.37625537334001 52.50926345498337)
        # values returned by SpatiaLite
        point_soldner_1924 = QgsPointXY(23099.49, 20296.69)
        point_soldner_1998 = QgsPointXY(23082.30, 20267.80)
        distance_soldner_meters = 33.617379
        azimuth_soldner_1924 = 3.678339
        # ST_Transform(point_soldner_1924,point_soldner_1998,4326)
        point_wsg84_1924 = QgsPointXY(13.37650707988041, 52.50952361017194)
        point_wsg84_1998 = QgsPointXY(13.37625537334001, 52.50926345498337)
        # ST_Distance(point_wsg84_1924,point_wsg84_1998,1)
        distance_wsg84_meters = 33.617302
        # ST_Distance(point_wsg84_1924,point_wsg84_1998)
        # distance_wsg84_mapunits=0.000362
        distance_wsg84_mapunits_format = QgsDistanceArea.formatDistance(0.000362, 7, QgsUnitTypes.DistanceDegrees, True)
        # ST_Azimuth(point_wsg84_1924,point_wsg84_1998)
        azimuth_wsg84_1924 = 3.674878
        # ST_Azimuth(point_wsg84_1998,point_wsg84_1998)
        azimuth_wsg84_1998 = 0.533282
        # ST_Project(point_wsg84_1924,33.617302,3.674878)
        # SRID=4326;POINT(13.37625537318728 52.50926345503591)
        point_soldner_1998_project = QgsPointXY(13.37625537318728, 52.50926345503591)
        # ST_Project(point_wsg84_1998,33.617302,0.533282)
        # SRID=4326;POINT(13.37650708009255 52.50952361009799)
        point_soldner_1924_project = QgsPointXY(13.37650708009255, 52.50952361009799)

        distance_qpoint = point_soldner_1924.distance(point_soldner_1998)
        azimuth_qpoint = point_soldner_1924.azimuth(point_soldner_1998)
        point_soldner_1998_result = point_soldner_1924.project(distance_qpoint, azimuth_qpoint)

        point_soldner_1924_result = QgsPointXY(0, 0)
        point_soldner_1998_result = QgsPointXY(0, 0)
        # Test meter based projected point from point_1924 to point_1998
        length_1998_mapunits, point_soldner_1998_result = da_3068.measureLineProjected(point_soldner_1924, distance_soldner_meters, azimuth_qpoint)
        self.assertEqual(point_soldner_1998_result.toString(6), point_soldner_1998.toString(6))
        # Test degree based projected point from point_1924 1 meter due East
        point_wsg84_meter_result = QgsPointXY(0, 0)
        point_wsg84_1927_meter = QgsPointXY(13.37652180838435, 52.50952361017102)
        length_meter_mapunits, point_wsg84_meter_result = da_wsg84.measureLineProjected(point_wsg84_1924, 1.0, (math.pi / 2))
        self.assertEqual(QgsDistanceArea.formatDistance(length_meter_mapunits, 7, QgsUnitTypes.DistanceDegrees, True), '0.0000147 deg')
        self.assertEqual(point_wsg84_meter_result.toString(7), point_wsg84_1927_meter.toString(7))

        point_wsg84_1998_result = QgsPointXY(0, 0)
        length_1928_mapunits, point_wsg84_1998_result = da_wsg84.measureLineProjected(point_wsg84_1924, distance_wsg84_meters, azimuth_wsg84_1924)
        self.assertEqual(QgsDistanceArea.formatDistance(length_1928_mapunits, 7, QgsUnitTypes.DistanceDegrees, True), distance_wsg84_mapunits_format)
        self.assertEqual(point_wsg84_1998_result.toString(7), point_wsg84_1998.toString(7))

    def processAlgorithm(self, parameters, context, feedback):
        layer = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.INPUT_VECTOR), context)

        rasterPath = str(self.getParameterValue(self.INPUT_RASTER))

        rasterDS = gdal.Open(rasterPath, gdal.GA_ReadOnly)
        geoTransform = rasterDS.GetGeoTransform()
        rasterDS = None

        fields = QgsFields()
        fields.append(QgsField('id', QVariant.Int, '', 10, 0))
        fields.append(QgsField('poly_id', QVariant.Int, '', 10, 0))
        fields.append(QgsField('point_id', QVariant.Int, '', 10, 0))

        writer = self.getOutputFromName(self.OUTPUT_LAYER).getVectorWriter(fields, QgsWkbTypes.Point,
                                                                           layer.crs(), context)

        outFeature = QgsFeature()
        outFeature.setFields(fields)

        fid = 0
        polyId = 0
        pointId = 0

        features = QgsProcessingUtils.getFeatures(layer, context)
        total = 100.0 / layer.featureCount() if layer.featureCount() else 0
        for current, f in enumerate(features):
            geom = f.geometry()
            bbox = geom.boundingBox()

            xMin = bbox.xMinimum()
            xMax = bbox.xMaximum()
            yMin = bbox.yMinimum()
            yMax = bbox.yMaximum()

            (startRow, startColumn) = raster.mapToPixel(xMin, yMax, geoTransform)
            (endRow, endColumn) = raster.mapToPixel(xMax, yMin, geoTransform)

            # use prepared geometries for faster intersection tests
            engine = QgsGeometry.createGeometryEngine(geom.geometry())
            engine.prepareGeometry()

            for row in range(startRow, endRow + 1):
                for col in range(startColumn, endColumn + 1):
                    (x, y) = raster.pixelToMap(row, col, geoTransform)
                    point = QgsPointXY()
                    point.setX(x)
                    point.setY(y)

                    if engine.contains(point):
                        outFeature.setGeometry(QgsGeometry(point))
                        outFeature['id'] = fid
                        outFeature['poly_id'] = polyId
                        outFeature['point_id'] = pointId

                        fid += 1
                        pointId += 1

                        writer.addFeature(outFeature, QgsFeatureSink.FastInsert)

            pointId = 0
            polyId += 1

            feedback.setProgress(int(current * total))

        del writer

 def testQgsPointXYRepr(self):
     p = QgsPointXY(123.456, 987.654)
     self.assertTrue(p.__repr__().startswith('<QgsPointXY: POINT(123.456'))

    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(QgsDistanceAr