Extract by location works improperly in PyQGIS standalone script

Question

I have been trying to extract polygons from one layer that intersect a polygon in another layer, and another extract that is disjointed, and looping this over several files organized by date.

When doing this in the QGIS GUI, it works as expected. The features from layer A that intersect layer B are extracted, as are the ones that are disjointed.

When I do the same procedure in my standalone script, the first date in the loop processes correctly, but for every following date the layer meant to contain intersecting features is empty, and all features are contained in the disjointed layer, including those that should be intersecting.

I have checked my code against the processing history in the GUI and see no obvious problems.

Here are some images to illustrate the problem: an album

Standalone script:

params = {'INPUT': 'C:/Users/ptbuffer{}{}20{}.shp'.format(dd, mm, yy),
          'PREDICATE': [0],
          'INTERSECT': 'C:/Users/dvdbuff{}{}20{}.shp'.format(dd, mm, yy),
          'OUTPUT': 'C:/Users/lochit{}{}20{}.shp'.format(dd, mm, yy)}

selectlochit = processing.run("qgis:extractbylocation", params)

params = {'INPUT': 'C:/Users/ptbuffer{}{}20{}.shp'.format(dd, mm, yy),
          'PREDICATE': [2],
          'INTERSECT': 'C:/Users/buff{}{}20{}.shp'.format(dd, mm, yy),
          'OUTPUT': 'C:/Users/locmiss{}{}20{}.shp'.format(dd, mm, yy)}

selectlocmiss = processing.run("native:extractbylocation", params)

Processing history:

processing.run("native:extractbylocation", {'INPUT':'C:/Users/ptbuffer30052011.shp',
       'PREDICATE':[0],
       'INTERSECT':'C:\\Users\\buff30052011.shp',
       'OUTPUT':'TEMPORARY_OUTPUT'})

processing.run("native:extractbylocation", {'INPUT':'C:/Users/ptbuffer30052011.shp',
       'PREDICATE':[2],
       'INTERSECT':'C:\\Users\\buff30052011.shp',
       'OUTPUT':'TEMPORARY_OUTPUT'})

I have the looped output written to a file like so, where condition is 0=disjoint, 1=intersect and 2 is the polygon in layer B. Area is in square meters. As you can see, the area of the intersecting polygons for every date but the first is -1.

DATE, CONDITION, AREA
20110529, 1, 13873909458.804
20110529, 0, 4970672583.033
20110529, 2, 908560692751.713
20110530, 1, -1.000
20110530, 0, 53305853234.103
20110530, 2, 1281098546690.083
20110531, 1, -1.000
20110531, 0, 4944271909.999
20110531, 2, 1459419315847.340

I have tested all possible parameter values (0-7) and none changed the result of the second "extract by location" operation.

I tested a sequence of 2 extract by locations outside of my loop to make sure my loop was not the source of error. The date shown in the example images represent the second "extract by location" shown in this code. The output of this test is the same as the original code, so the loop is not the problem. The test code and output follows:

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt5.QtCore import *
from PyQt5.QtGui import *
import calendar as c
import pandas as pd
import urllib.request
from timezonefinder import TimezoneFinder
import pendulum
tf = TimezoneFinder()
print("Import successful")


print("Output file initialized")

qgs = QgsApplication([], False)
QgsApplication.initQgis()

print("QGIS initialized")

print("Initializing processing...")
import processing
from processing.core.Processing import Processing
Processing.initialize()
from qgis.analysis import QgsNativeAlgorithms
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

print("Processing initialized")

blayer = QgsVectorLayer('C:/Users/ptbuffer29052011.shp', "ptbuffer", "ogr")
if not blayer.isValid():
    print("Layer fail!")
print("Point buffer layer added successfully")
QgsProject.instance().addMapLayer(blayer)

polyblayer = QgsVectorLayer('C:/Users/buff29052011.shp', "buffer", "ogr")
if not polyblayer.isValid():
    print("Layer fail!")
QgsProject.instance().addMapLayer(polyblayer)

# SELECT BY LOCATION
params1 = {'INPUT': 'C:/Users/ptbuffer29052011.shp',
          'PREDICATE': [0],
          'INTERSECT': 'C:/Users/buff29052011.shp',
          'OUTPUT': 'C:/Users/lochit29052011.shp'}
selectlochit = processing.run("qgis:extractbylocation", params1)
lochit = QgsVectorLayer(selectlochit['OUTPUT'], "selhit", "ogr")
if not lochit.isValid():
    print("Layer fail!")
QgsProject.instance().addMapLayer(lochit)
params2 = {'INPUT': 'C:/Users/ptbuffer29052011.shp',
          'PREDICATE': [2],
          'INTERSECT': 'C:/Users/buff29052011.shp',
          'OUTPUT': 'C:/Users/locmiss29052011.shp'}
selectlocmiss = processing.run("native:extractbylocation", params2)
locmiss = QgsVectorLayer(selectlocmiss['OUTPUT'], "selmiss", "ogr")
if not locmiss.isValid():
    print("Layer fail!")
QgsProject.instance().addMapLayer(locmiss)

print("Hit Feature count: ", lochit.featureCount()) # OUTPUT: 3
print("Miss Feature count: ", locmiss.featureCount()) # OUTPUT: 2


blayer = QgsVectorLayer('C:/Users/ptbuffer30052011.shp', "ptbuffer", "ogr")
if not blayer.isValid():
    print("Layer fail!")
print("Point buffer layer added successfully")
QgsProject.instance().addMapLayer(blayer)

polyblayer = QgsVectorLayer('C:/Users/buff30052011.shp', "buffer", "ogr")
if not polyblayer.isValid():
    print("Layer fail!")
QgsProject.instance().addMapLayer(polyblayer)

params1 = {'INPUT': 'C:/Users/ptbuffer30052011.shp',
       'PREDICATE': [0],
       'INTERSECT': 'C:/Users/buff30052011.shp',
       'OUTPUT': 'C:/Users/lochit30052011.shp'}
selectlochit = processing.run("qgis:extractbylocation", params1)
lochit = QgsVectorLayer(selectlochit['OUTPUT'], "selhit", "ogr")
if not lochit.isValid():
    print("Layer fail!")
QgsProject.instance().addMapLayer(lochit)
params2 = {'INPUT': 'C:/Users/ptbuffer30052011.shp',
       'PREDICATE': [2],
       'INTERSECT': 'C:/Users/buff30052011.shp',
       'OUTPUT': 'C:/Users/locmiss30052011.shp'}
selectlocmiss = processing.run("native:extractbylocation", params2)
locmiss = QgsVectorLayer(selectlocmiss['OUTPUT'], "selmiss", "ogr")
if not locmiss.isValid():
    print("Layer fail!")
QgsProject.instance().addMapLayer(locmiss)

print("Hit Feature count: ", lochit.featureCount()) # OUTPUT: 0
print("Miss Feature count: ", locmiss.featureCount()) # OUTPUT: 17

Any ideas on what I'm doing wrong here?

Answer 1

I solved the problem, though I'm still unsure what caused it.

I was testing to see if there were any projection issues, and found that there were with the following code:

    print(lochit.crs().isValid()) # True
    print(lochit.crs().authid()) # -blank-

For every date beyond the first in the loop, the crs is found to be valid, but authid() has no output. I then used the following code after the first processing step:

    if dd != ddi:   # checks if first iteration of loop
        crs = rplayer.crs()
        crs.createFromId(4326)  
        rplayer.setCrs(crs)

Without the first iteration check, the same problem occurs but only with the first date in the loop. This workaround, somehow, solves the problem. I don't understand what it means for a layer to have a valid CRS but for authid() to have no output, or why problems with the projections would only arise after the first date. I'm satisfied with this workaround, but any insight is welcome and appreciated.