Finding Greenspaces

OS Open Greenspace is an open dataset of greenspaces in Britain. You can use it to understand the location of public parks, playing fields, sports facilities and allotments along with access points for entering and exiting urban and rural greenspaces.

The data is available as 100km2 tiles. For LOCAL, I need to find greenspace cover at ward level within local authorities in the UK.

Method

• Create polygons from OS Open Greenspace shapefile
• Create polygons representing ward areas using ONS ward boundary data
• For each ward, find the intersection between ward boundary and greenspace polygons.
• The intersection area represents greenspace in the ward.

Here the intersection between entire greenspace data and 9k+ ward boundaries will take days or weeks to complete(per my test run on an AMD Ryzen 5 quad-core server). We can speed up the process by replacing the slow sequential scan with a spatial index lookup.

• Build an STRtree index of greenspaces.
• Find the candidate greenspaces envelopes containing ward boundary.
• Perform intersection between candidate greenspaces and ward boundaries. 

With the introduction of STRtree index, the process runs within minutes!

Python code is below. Uses Fiona library for reading shapefiles and shapely for geometric operations.

def greens_index(shapefile):
    """
    Build an index of green areas for looking up possible intersection
    candiates with ward.
    """
    c = fiona.open(shapefile, 'r')
    polygons = [shape(x['geometry']) for x in c]
    return STRtree(polygons), c.crs

def wards_iter(shapefile):
    """
    Yield ward records.

    Parameters
    ----------
    shapefile: Input ward boundary file

    Returns
    -------
    Yields tuples in the form (osward, osward_desc, bounday_polygon) where
    osward is the ons code for ward. 
    """
    wards = fiona.open(shapefile, 'r')

    for ward in wards:
        osward = ward['properties']['wd19cd']
        osward_desc = ward['properties']['wd19nm']

        polygons = shape(ward['geometry'])
        boundary = unary_union(polygons).buffer(0)
        yield osward, osward_desc, boundary

def greenspace_intersect(tree, boundary):
    """
    Find greenspace coverage in the given area - ward, geohash5 bbox etc.

    Parameters
    ----------
    tree: STRtree index of green area polygons
    boundary: boundary shape of a ward.

    Returns
    -------
    A shape representing green space in given boundary.

    Below we explicitly check for intersects because tree.query gives possible
    envelopes that may overlap. To check if the shape instead of envelope
    intersects, you need to explicitly check for intersection.

    """
    greens = [g for g in tree.query(boundary) if g.intersects(boundary)]
    # without unary union to merge, we will count overlapping areas multiple
    # times.
    # buffer(0) is the standard approach to draw an envelope around all shapes
    # to correct invalid geometries
    greens_merged = unary_union(greens).buffer(0)
    greens_ward = greens_merged.intersection(boundary)
    return greens_ward

# main 

tree = greens_index(shp_greenspaces)

for osward, osward_desc, boundary_ward in wards_iter(shp_ward):
    greens_in_ward = greenspace_intersect(tree, boundary_ward)
    
    ...... other ops, truncated