geoprocessing across multiple vector layers QGIS2

Question

I have 6 layers of point data (each layer is one year), where data in each layer are distributed throughout a grid, and the grid is identical in each layer. I'd like to extract data from each layer, excluding those data where there is not a full six years of data in the grid square. For example, if a grid square does not have 6 years of data, then I'm not interested in it. If a grid does have data for all 6 years, then I want to keep it. The result will be a series of layers where the table for each layer has the same number of grid square references.

I've explore database options with this - my data is such that for each grid square per year, there are multiple (ie daily) occurrences within that year, so the Access query becomes very weird quickly, so I thought there might be a gis solution. I considered the geoprocessing tools like clipping or intersection, but it seems you can only do 2 layers at a time. I'm using QGIS2 but I'm open to a GRASS solutions as well.

Any pointers in the right direction would be great. Many thanks, Paul.

Answer 1

It is possible in Python, pure Python without QGIS with modules such as Fiona, Pyshp, shapely with rtree (see rtree python polygon index) and with PyQGIS.

The problem: 6 points layers and an unique grid:

A possible solution with PyQGIS:

1) Using a spatial index for the grid (see Using Spatial Index in Using Vector Layers, Using a QGIS spatial index to speed up your code or How to do a spatial search without select() using PyQGIS?):

grid = iface.activeLayer()
index = QgsSpatialIndex() 
for elem in grid.getFeatures():
     index.insertFeature(elem)

or more "Pythonic"

index = QgsSpatialIndex() 
map(index.insertFeature, grid.getFeatures())

2) Creation of a dictionary to collect the results (index of the grid squares where there are points inside and layers) and a function to complete it:

results = {}
# fonction to find the indexes of the grid squares inside which a point lies.
def mydict(layer, dictionary):
    for elem in layer.getFeatures():
        geom = elem.geometry().asPoint()
        nearestIds = index.nearestNeighbor(geom,1)
        results.setdefault(str(nearestIds[0]), set()).add(layer.name())

3) Iteration over point layers and application of the function:

canvas= qgis.utils.iface.mapCanvas()
for layer in canvas.layers():
    elem = layer.getFeatures().next()
    geom = elem.geometry()
    # only points layers
    if geom.wkbType() == QGis.WKBPoint:
         mydict(layer, results)

The result is a dictionary with the grid squares spatial index as key and the layers which have points in the square grid as values:

print results
{'11': set([u'point1', u'point3', u'point4']),
 '10': set([u'point5']),
 '13': set([u'point4']), 
 '12': set([u'point3', u'point5', u'point6']),
 '15': set([u'point4', u'point6']), 
 '14': set([u'point1', u'point2', u'point3', u'point4', u'point5', u'point6']),
  '1': set([u'point4', u'point5', u'point6']), 
  '0': set([u'point1', u'point2', u'point3', u'point4', u'point5', u'point6']), 
  '3': set([u'point4', u'point5']), 
  '2': set([u'point3', u'point6']), 
  '5': set([u'point4']),
  '4': set([u'point3', u'point5']),
  '7': set([u'point3', u'point6']), 
  '6': set([u'point1', u'point3']), 
  '9': set([u'point3']), 
  '8': set([u'point1', u'point2'])}

If you want to preserve the order of the indexes, you can use an OrderedDict

4) And if you only want the grid squares with six years of data:

for index_gridsquare, layers in results.iteritems():
    if len(values) == 6:
         print index_gridsquare,": ", values
14: set([u'point1', u'point2', u'point3', u'point4', u'point5', u'point6'])
0 : set([u'point1', u'point2', u'point3', u'point4', u'point5', u'point6'])

After that, I do not understand if you want to select:

• the square grids:

• the points of the six layers present in the square grids:

Answer 2

Here's an approach that doesn't require an programming.

I'm assuming that you have 6 point layers (year1 .. year6) and one polygon layer (grid).

Go to Vector > Research Tools > Select by location.

Select features in the grid that intersect (i.e. touch) features in the points layer.

Enable editing of the grid layer, open it's attribute table and open the field calculator. Create a new field named year1 and set the value to 0 for all records. Using the field calculator again, update the year1 field and set the value to 1 for selected features only.

An alternate approach would be to concatenate the years in a single field. e.g. in the field calculator: concat("years",'1,'), then concat("years",'2,'), where years started out as an empty field. You'll end up with field values like 1,3,4,6,, which should be easy to style/filter.