1

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.

enter image description here

  • In your explanation of the problem, I'm afraid I lost you at the part "The result will be a series of layers...". Do you want to aggregate all 6 years of point data for each grid square? or do you want to filter out those points for which some years are missing? Should the resulting layers be point layers? – Micha Oct 31 '13 at 14:29
  • I see the confusion. I've added an example of what I mean. So in the example, years 1 and 2 are different point layers (each point is georeferenced to the grid). I'd like to filter year 1 and 2 for only points where they share the same grid. So the results will be 2 x point layers, where data in each grid square of each layer has the full years. For me, the problem is difficult because its 6 layers. hope that helps? – Paul Oct 31 '13 at 15:17
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:

enter image description here

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:

enter image description here

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

enter image description here

0

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.

Select by location

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.