Checking intersection between two layers with different projection using PyQGIS

Question

I have a point layer and polygon layer

polygon layer contains espg:6344

point layer contains espg:4326

project Crs is espg: 4326

polygon Layer contains rectangles and points are the centroids of rectangles

now on map canvas data show correctly if i used select by location intersected points selected

but when try with python its intersection false

poly=QgsProject().intance().mapLayersByName("polygons")[0]
pnt=QgsProject().intance().mapLayersByName("points")[0]

geom_poly=poly.selectedFeatures()[0].geometry()
geom_pnt=pnt.selectedFeatures()[0].geometry()

print(geom_pnt.intersects(geom_poly))

output: >>> False

above due to projection. then i converted the projection both espg 4326 and espg 6344 but not work.

How to get intersected points with polygon if both are different projection?

Answer 1

Firstly, your layers are shown in the correct location on the map canvas because QGIS re-projects layers on the fly to the project CRS for visualization.

Secondly, you have not said how you "converted" the CRS of your layers but, if you right clicked on the layer in the Table of Contents and used 'Set Layer CRS' option, that is not the way to change layer CRS. If you did that you should set them back to their correct CRS before doing anything else. To re-project a layer you can either save/export the layer and choose a different target CRS or use the Reproject Layer algorithm from the Processing Toolbox.

Solution:

Since your question is about PyQGIS, you can easily work with layers in different coordinate reference systems in the code without having to reproject layers. To do this, you can simply transform one of the QgsGeometry objects you are working with to match the CRS of the other, using the QgsCoordinateTransform class, passing a source CRS, destination CRS and QgsProject instance to the class constructor. Then call transform() on whichever geometry is in the source CRS.

Like this:

poly = QgsProject.instance().mapLayersByName('polygons')[0]
pnt = QgsProject.instance().mapLayersByName('points')[0]

geom_poly = poly.selectedFeatures()[0].geometry()
geom_pnt = pnt.selectedFeatures()[0].geometry()

# Create a QgsCoordinateTransform object, passing source CRS, dest CRS & QgsProject
xform = QgsCoordinateTransform(pnt.crs(), poly.crs(), QgsProject.instance())
# transform pnt_geom to poly_geom CRS (pass QgsCoordinateTransform object to the transform() call
geom_pnt.transform(xform)

# now should print True
print(geom_pnt.intersects(geom_poly))

See the relevant parts of documentation for QgsCoordinateTransform and QgsGeometry classes:

https://api.qgis.org/api/classQgsCoordinateTransform.html#aefa41b71beec7d5e46dfd4d04655fc4a

https://api.qgis.org/api/classQgsGeometry.html#a45a9753a80833dfb65f05b93ed019a57

And if you want to loop through all features in both layers to find the points which intersect with each polygon, you can use this approach. This should be fast and efficient since we are using a QgsSpatialIndex.

poly = QgsProject.instance().mapLayersByName('polygons')[0]
pnt = QgsProject.instance().mapLayersByName('points')[0]

idx = QgsSpatialIndex(pnt.getFeatures())

# Create a QgsCoordinateTransform object, passing source CRS, dest CRS & QgsProject
xform = QgsCoordinateTransform(poly.crs(), pnt.crs(), QgsProject.instance())

for f in poly.getFeatures():
    geom_poly = f.geometry()
    geom_poly.transform(xform)
    bb = geom_poly.boundingBox()
    for ft in pnt.getFeatures(idx.intersects(bb)):
        if ft.geometry().intersects(geom_poly):
            #do something...

Answer 2

if you want to loop into all the features from one layer and compare them to another layer you can set up a QgsFeatureRequest to handle an automatic reprojection of the features on the fly.

You have to use the setDestinationCrs method and pass the crs you want to reproject the features in and the QgsCoordinateTransformContext taht you can grab directly from QgsProject.instance().transformContext().

Here a small code snippet:

project = QgsProject.instance()
poly = project.mapLayersByName("polygons")[0]
pnt = project.mapLayersByName("points")[0]

request = QgsFeatureRequest()
request.setDestinationCrs(pnt.crs(), project.transformContext())

for poly_feature in poly.getFeatures(request):
    poly_geom = poly_feature.geometry()

    for pnt_feature in pnt.getFeatures():
        pnt_geom = pnt_feature.geometry()

        if poly_geom.intersects(pnt_geom):
            #do whatever you want