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?
'PREDICATE': [1]
? Whta if you change'PREDICATE': [1]
into'PREDICATE': "['intersects']"
or'PREDICATE': "[u'intersects']"
? Do that problematic layers give any error on their geometry consistency?[1]
as predicate? I do have[0]
in QGIS 3.4