find_landsat_and_sentinel.py

import sys
import json
import datetime

import click

import pandas as pd
import geopandas as gpd
from shapely.geometry import box, mapping

# LANDSAT_METADATA_URL = 'http://storage.googleapis.com/gcp-public-data-landsat/index.csv.gz'
# SENTINEL2_METADATA_URL = 'http://storage.googleapis.com/gcp-public-data-sentinel-2/index.csv.gz'

s2_path = "/Users/vincentsarago/Workspace/sentinel.csv"
l8_path = "/Users/vincentsarago/Workspace/landsat.csv"
crs = {'init': 'epsg:4326'}


def create_feature(scene, selection):
    return {
        'type': 'Feature',
        'properties': {
            'L8_SCENE_ID': scene.SCENE_ID,
            'L8_PRODUCT_ID': scene.PRODUCT_ID,
            'L8_SENSING_TIME': scene.SENSING_TIME.isoformat(),
            'S2_PRODUCT_ID': list(selection.PRODUCT_ID),
            'S2_SENSING_TIME': [x.isoformat() for x in selection.SENSING_TIME],
        },
        'geometry': mapping(scene.geometry)}


def main():
    l8_meta = pd.read_csv(l8_path)
    l8_meta.drop(columns=["SENSOR_ID", "COLLECTION_NUMBER", "COLLECTION_CATEGORY", "DATA_TYPE", "TOTAL_SIZE"], inplace=True)
    l8_meta = l8_meta.loc[l8_meta.SPACECRAFT_ID == "LANDSAT_8"]
    l8_meta['SENSING_TIME'] = pd.to_datetime(l8_meta['SENSING_TIME'])
    l8_meta = l8_meta.loc[
        l8_meta.SENSING_TIME > datetime.datetime(2015, 6, 23)
    ]
    l8_meta = l8_meta[(l8_meta.NORTH_LAT < 60) & (l8_meta.SOUTH_LAT > -60)]

    geometry = [
        box(*bbox) for bbox in zip(
            l8_meta.WEST_LON,
            l8_meta.SOUTH_LAT,
            l8_meta.EAST_LON,
            l8_meta.NORTH_LAT
        )
    ]
    l8_meta.drop(columns=["WEST_LON", "SOUTH_LAT", "EAST_LON", "NORTH_LAT"], inplace=True)
    l8_catalog = gpd.GeoDataFrame(l8_meta, crs=crs, geometry=geometry)
    del geometry, l8_meta

    s2_meta = pd.read_csv(s2_path)
    s2_meta.drop(columns=["DATATAKE_IDENTIFIER", "TOTAL_SIZE", "GEOMETRIC_QUALITY_FLAG", "GENERATION_TIME"], inplace=True)
    s2_meta['SENSING_TIME'] = pd.to_datetime(s2_meta['SENSING_TIME'])
    s2_meta = s2_meta.loc[(s2_meta.NORTH_LAT < 60) & (s2_meta.SOUTH_LAT > -60)]

    geometry = [
        box(*bbox) for bbox in zip(
            s2_meta.WEST_LON,
            s2_meta.SOUTH_LAT,
            s2_meta.EAST_LON,
            s2_meta.NORTH_LAT
        )
    ]
    s2_meta.drop(columns=["WEST_LON", "SOUTH_LAT", "EAST_LON", "NORTH_LAT"], inplace=True)
    s2_catalog = gpd.GeoDataFrame(s2_meta, crs=crs, geometry=geometry)
    del geometry, s2_meta

    geojson = {
        "type": "FeatureCollection",
        "features": []
    }
    with click.progressbar(l8_catalog.iterrows(), length=len(l8_catalog), file=sys.stderr, show_percent=True) as rows:
        for _, scene in rows:
            selection = s2_catalog[(scene['SENSING_TIME'] - s2_catalog['SENSING_TIME']).abs() < pd.Timedelta('10 seconds')]
            # intersects
            inter = selection[selection.intersects(scene.geometry)]
            if len(inter):
                geojson["features"].append(create_feature(scene, inter))

    with open('L8_S2_30seconds.geojson', 'w') as f:
        f.write(json.dumps(geojson))


if __name__ == '__main__':
    main()