Combined area of all features within features of another layer

Question

I have a grid shapefile (basically, a half-degree grid of the world, where each cell has a unique id number), and a set of files with features. I want to record for each grid cell, all the areas and/or values of one field of the second set of files that fall within the gridcell. The "fall within" is whether the centroid is within the cell. My code looks like this (Python+OGR, I can't get shapely installed on the work system)

    # grid_layer is my grid layer, it has been opened elsewhere
    # outut_result is a dict indexed by grid ID#
    feat_shape = ogr.Open ( shapefile )
    feat_l = feat_shape.GetLayer ( 0 )
    extent = feat_l.GetExtent()
    # Apply the extent filter to the grid layer
    grid_layer.SetSpatialFilterRect ( *extent )
    # Now loop over each grid cell within the feat region
    for grid in grid_layer:
        grid_geom = grid.GetGeometryRef()
        grid_no = grid.GetFieldAsInteger ('ID')
        # Loop over features within this grid id
        # Start by applying a spatial filter
        feat_l.SetSpatialFilterRect ( *grid_geom.GetEnvelope())
        for feat in feat_l:
            # Get the feat's geometry and calculate its centroid
            feat_geom = feat.GetGeometryRef()
            centroid = feat_geom.Centroid()
                # Is the centroid within the grid_geom?
                if centroid.Within ( grid_geom ):
                    value = feat.GetFieldAsInteger ( 'field_of_interest' )
                    # Or feature's area, or whatever
                    # Store it in our output dict
                    output_result.setdefault(grid_no, []).append( value )


            del centroid
            del feat_geom

        del grid_geom

At the end of this process, I am left with output_result and my information: output_result[29345] results in a list with all the values of field_of_interest, which I then can process further with e.g. numpy or scipy.

My question is that this processing takes forever, as there are a lot of shapefiles. Then writing it out is rather messy (doing it with a pickle at the moment, but that needs even more code afterwards). Am I missing some really obvious "point-and-clicky" way of solving this problem?

Answer 1

It is inefficient to be looping over the entire grid for each shapefile. You can take advantage of the regular structure of your grid to speed up the processing tremendously: from the coordinates of a feature's centroid you can mathematically compute the identifier of the grid cell in which it falls. Output a tuple of this identifier and the feature's attribute, then proceed to the next feature. This is a fast, constant-time operation for each feature. Later on you can sort by grid cell identifier or perform summaries keyed on the cell identifier: that's quick and straightforward.

The grid cell containing a point (x,y) (with -180 <= y < 180 and -90 < x < 90) is identified with the ordered pair of integers

( Floor((x+90)/0.5), Floor((y+180)/0.5) )

In general, cells in a grid with mesh (h,k) and origin (x0,y0) are identified by the tuples

( Floor((x-x0)/h), Floor((y-y0)/k) ).