Proof of concept for My AOI analysis

aoi.geojson

data.vrt

<OGRVRTDataSource>
    <OGRVRTUnionLayer name="data">
        <OGRVRTLayer name="tile1">
            <SrcDataSource>/vsis3/gfw2-data/alerts-tsv/country-pages/ten-by-ten-gadm36/wdpa__10S_060W.tsv</SrcDataSource>
            <SrcLayer>wdpa__10S_060W</SrcLayer>
            <GeometryType>wkbPolygon</GeometryType>
            <GeometryField encoding="WKT" field='field_1'/>
        </OGRVRTLayer>
        <OGRVRTLayer name="tile2">
            <SrcDataSource>/vsis3/gfw2-data/alerts-tsv/country-pages/ten-by-ten-gadm36/wdpa__00N_060W.tsv</SrcDataSource>
            <SrcLayer>wdpa__00N_060W</SrcLayer>
            <GeometryType>wkbPolygon</GeometryType>
            <GeometryField encoding="WKT" field='field_1'/>
        </OGRVRTLayer>
        <OGRVRTLayer name="tile3">
            <SrcDataSource>/vsis3/gfw2-data/alerts-tsv/country-pages/ten-by-ten-gadm36/wdpa__00N_070W.tsv</SrcDataSource>
            <SrcLayer>wdpa__00N_070W</SrcLayer>
            <GeometryType>wkbPolygon</GeometryType>
            <GeometryField encoding="WKT" field='field_1'/>
        </OGRVRTLayer>
        <OGRVRTLayer name="tile4">
            <SrcDataSource>/vsis3/gfw2-data/alerts-tsv/country-pages/ten-by-ten-gadm36/wdpa__10S_070W.tsv</SrcDataSource>
            <SrcLayer>wdpa__10S_070W</SrcLayer>
            <GeometryType>wkbPolygon</GeometryType>
            <GeometryField encoding="WKT" field='field_1'/>
        </OGRVRTLayer>
    </OGRVRTUnionLayer>
</OGRVRTDataSource>

intersect.py

import argparse
import subprocess
import json

import fiona
from shapely.geometry import shape


def main():

    parser = argparse.ArgumentParser(description='Intersect an input AOI with a polyname of interest')
    parser.add_argument('--aoi', '-a', help='The input AOI as a geojson file', required=True)
    parser.add_argument('--polyname', '-p', help='The polyname of interest', required=True)
    args = parser.parse_args()

    with open(args.aoi) as src:
        geom = shape(json.load(src)['features'][0]['geometry'])

    # will ultimately want to build a dynamic Union VRT from these tiles
    # for now i'll just hardcode this
    # VRT info: https://www.gdal.org/drv_vrt.html
    print find_tiles(geom)

    subprocess.check_call(['ogr2ogr', '-f', 'GeoJSON', 'output.geojson',
                           'data.vrt', 'data', '-clipsrc', args.aoi])


def find_tiles(geom):

    tiles = '/vsis3/gfw2-data/alerts-tsv/country-pages/loss_data_footprint.geojson'
    int_tiles = []

    with fiona.open(tiles, 'r', 'GeoJSON') as tiles:
        for tile in tiles:
            if shape(tile['geometry']).intersects(geom):
                tile_dict = tile['properties']
                tile_name = tile_dict['id']

                int_tiles.append(tile_name)

    return int_tiles


if __name__ == '__main__':
    main()

output.geojson