bounds-dense.py

#!/usr/bin/env python
import json
import os

import click
import numpy as np
import rasterio
from rasterio.warp import transform


def transform_dense_bounds(src_crs, dst_crs, left, bottom, right, top,
                           densify_pts=21):
    """
    like warp.transform_bounds but returns the densified coordinate array
    """

    if densify_pts < 0:
        raise ValueError('densify parameter must be >= 0')

    in_xs = []
    in_ys = []

    if densify_pts > 0:
        densify_factor = 1.0 / float(densify_pts + 1)

        # Add points along outer edges, starting at the lower left

        # left edge
        in_xs.extend([left] * (densify_pts + 2))
        in_ys.extend(
            bottom + np.arange(0, densify_pts + 2, dtype=np.float32) *
            ((top - bottom) * densify_factor)
        )

        # top
        in_xs.extend(
            left + np.arange(0, densify_pts + 2, dtype=np.float32) *
            ((right - left) * densify_factor)
        )
        in_ys.extend([top] * (densify_pts + 2))

        # right
        in_xs.extend([right] * (densify_pts + 2))
        in_ys.extend(
            top - np.arange(0, densify_pts + 2, dtype=np.float32) *
            ((top - bottom) * densify_factor)
        )

        # bottom
        in_xs.extend(
            right - np.arange(0, densify_pts + 2, dtype=np.float32) *
            ((right - left) * densify_factor)
        )
        in_ys.extend([bottom] * (densify_pts + 2))

    else:
        in_xs = [left, left, right, right, left]
        in_ys = [bottom, top, top, bottom, bottom]

    xs, ys = transform(src_crs, dst_crs, in_xs, in_ys)

    # close the polygon
    xs.append(xs[0])
    ys.append(ys[0])

    return list(zip(xs, ys))


# Bounds command.
@click.command(name="bounds-dense",
               short_help="Write densified bounding boxes to stdout as GeoJSON.")
@click.argument('INFILE', nargs=1, type=click.Path())
def bounds_dense(infile):
    with rasterio.drivers(CPL_DEBUG=True):
        with rasterio.open(infile) as src:
            coords = transform_dense_bounds(
                src.crs, {'init': 'epsg:4326'}, *src.bounds, densify_pts=21)

        click.echo(json.dumps({
            'type': 'Feature',
            'geometry': {
                'type': 'Polygon',
                'coordinates': [coords]},
            'properties': {
                'title': infile,
                'filename': os.path.basename(infile)}}))


if __name__ == "__main__":
    bounds_dense()

query_footprints.py

#!/usr/bin/env python
import json

import click
import cligj
from shapely.geometry import asShape

@click.command()
@click.argument("footprints", type=click.File('r'))
@cligj.features_in_arg
def query(footprints, features):
    fc = json.loads(footprints.read())
    footprints = [(f, asShape(f['geometry'])) for f in fc['features']]

    matches = []
    for qfeat in features:
        qshape = asShape(qfeat['geometry'])
        # TODO need spatial index
        for feat, shape in footprints:
            if qshape.intersects(shape):
                matches.append(feat)

    for match in matches:
        click.echo(json.dumps(match))


if __name__ == "__main__":
    query()