automating forest cover calculations with Python

Question

I have a layer of forest cover comprising ~ 30,000 individual forest patches, and I need to calculate the amount of forest cover within a 5 km buffer of each patch. The buffers all overlap each other, which stopped me from calculating forest cover with a simple intersect of the buffer layer with the forest layer. In a small test area I have been able to calculate amount of forest by doing a Union on the forest cover and buffer layers, using a definition query on the resulting layer to only show features where buffers and patches intersect, then Spatial Joining this layer to the original buffer layer with a 'sum' merge rule to give the total amount of forest per buffer. However, the largest area I have been able to apply this process to successfully is about 1/1000th the total area that I need to do it to, and even this takes about 1hr of processing to complete. Beyond this the initial Union fails - I think that this is because all the overlapping buffers create an enormous number of features when Union runs.

So, I’m wondering if it will be possible to break this process down into smaller chunks with a Python script? I’m a raw beginner with Python (aside from a few online tutorials I’ve just done) but guessing it might be possible to individually intersect each of the 30,000 buffers with the native forest layer, sequentially write the output of each of these intersections into e.g. a text file, then joining this to the buffer layer’s attribute table. Before I invest too much effort into learning scripting, can anyone tell me if this approach is likely to be doable for a Python newbie or am I barking up the wrong tree?

Thanks in advance,

Jay

Answer 1

EDIT: Just realized you probably wanted an ArcPy script, or at least something Arc related. I haven't used ArcPy in a long time, but it could be possible following the same logic as what I've said below to do this. Have a look at some of the ArcPy documentation on ESRIs website, it's very informative. Not sure how computationally efficient it will be in Arc however.

I can't offer a complete script but hopefully I can point you in the right direction.

Firstly, I think this could be possible with GDAL python bindings. Assuming you're using windows you can get a precompiled windows and python bindings in an easy installer at http://vbkto.dyndns.org/sdk/ or alternatively you can just download OSGeo4W. If you're not familiar with GDAL python bindings, http://www.gis.usu.edu/~chrisg/python/2009/ is a great intro.

So having done this i'd assume you'd want to do something like this:

Firstly

import ogr from osgeo

1) Create a driver and load the forest layer in

driver = ogr.GetDriverByName('ESRI Shapefile')
forestcover = driver.Open('forest.shp', 1)

2) Create a new layer which is a union of all 30'000 forest features. Load this layer in also.

driver = ogr.GetDriverByName('ESRI Shapefile')
forestunioned = driver.Open('forestunioned.shp', 1)

3) Loop through each feature in the original file with a for loop. Create a buffer 5m for each feature, whilst in the loop and intersect this buffer with the unioned forest features.

4) Bind the resultant intersection geometry to a variable

5) find the area of the resultant geometry

6) write the resultant area into a field for the feature

7) Loop through all the features using this method

Here is a sample of what the code might look like:

unionedlayer = forestunioned.GetLayer()
unionedgeom = unionedlayer.GetGeometryRef()
layer = forestcover.GetLayer()
feature = layer.GetNextFeature()
while feature:
    featgeom = feature.GetGeometryRef()
    bufferedfeat = featgeom.Buffer(5000) #not sure what the unit is for distance in OGR
    intersection = bufferedfeat.Intersection(unionedgeom)
    #Here you need to get the area somehow not sure if intersection.GetArea() works
    #Get the value as a integer
    #Write it the features attribute table 
    feature = layer.GetNextFeature()
layer.ResetReading()