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 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?