Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question

1 Answer 1

up vote 5 down vote accepted

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) ).
share|improve this answer
    
Aaarghh!!! Absolutely right! In the particular case of a grid, it's straightforward. All I need from this is a mapping from the two integer grid location to the actual grid id –  Jose Jun 23 '11 at 16:35
1  
@Jose That's right. An easy way to make that mapping is to create a new id in the grid dataset that will readily match the ids you calculate in the program; e.g., take the x-index (between 0 and 179) and add it to 180 times the y-index. –  whuber Jun 23 '11 at 16:38
    
For reference, in this particular case it works well, but it wouldn't in another case where I no longer have a grid, but irregular polygons. In that second case, would there be any additional speed ups you can think of? –  Jose Jun 23 '11 at 16:41
1  
@Jose It depends on a lot of things. You will either be selecting "grid" features with point-in-polygon searches or you will be selecting shapefile features based on containing grid features. Sometimes one of those searches is substantially faster than the other. E.g., in case there are relatively few grid features and they use few vertices, you can build a good spatial index and cache it for the duration. If the shapefile features are numerous and complex, it will then be faster to search for grid features than to search for shapefile features. –  whuber Jun 23 '11 at 16:47
    
Choosing your answer as it is really exactly what I needed. Thanks! –  Jose Jun 23 '11 at 17:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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