kylebarron / landsat-cogeo-mosaic

Create mosaicJSON for Landsat imagery

landsat-cogeo-mosaic

Tools to create MosaicJSONs for Landsat imagery.

Overview

A MosaicJSON is a file that defines how to combine multiple (satellite) imagery assets across time and space into a web mercator tile. This repository is designed to be used to create such files, so that they can be used for on-the-fly satellite tile generation, such as with landsat-mosaic-tiler.

Documentation Website

Install

pip install landsat-cogeo-mosaic

For CLI commands and some modules, you need additional dependencies:

pip install 'landsat-cogeo-mosaic[cli]'

Create Mosaics

In order to create the mosaic, you need the metadata on where and when Landsat took images, and also the cloud cover of each image. There are a couple ways to get this data. One is to ping the Sat API hosted by Development Seed. This allows you to find the metadata from specific regions and times, and is easier when you're creating a mosaic for a specific region.

Alternatively, for large mosaics it's easier to do bulk processing. In order to not overload the Sat API, I download a file from AWS with metadata from all images.

landsat-cogeo-mosaic's Issues

Avoid hitting recursion limit

Set recursion limit higher, just for calling sat api? Avoid with while?

Requirements.txt is missing some requirements

After running pip3 install landsat-cogeo-mosaic i expect to be able to run landsat-cogeo-mosaic, but i also had to install geopandas and keplergl_cli before it would work.

Use same-year imagery when possible

When you optimize for fewest assets and most overlap, sometimes the imagery can look like this:

Mosaic validation

• Find quadkey tiles over land
• Find geometries of mosaic by finding the pathrow of every asset and looking up from the pathrow dict

Plot

Work with any STAC, not just landsat

Would just require a change here:

landsat-cogeo-mosaic/landsat_cogeo_mosaic/mosaic.py

Lines 128 to 131 in 795459d

	# Add to mosaic definition
	mosaic_definition["tiles"][quadkey] = [
	scene["properties"]["landsat:product_id"] for scene in assets
	]

So that the correct product ID is extracted for whatever the STAC API is

Optimize better than unique pathrow

The blue grid is zoom 8 mercator tiles, the brown is landsat pathrows.

For example with the middle tile, you only need 2 pathrows to cover the tile. but with only optimizing by pathrow, you get 6. At moderate latitudes, including most of the US, you can optimize way better (at zoom 8 quadkeys) than just taking the unique pathrow.

You should make this optimization for all prebuilt mosaics.

Repro code:

import json
import mercantile
from keplergl_cli import Visualize
from geojson import GeometryCollection, FeatureCollection, Feature
from shapely.geometry import Polygon, asShape, mapping

with open('data/pr2coords.json') as f:
    pr_xw = json.load(f)

list_depth = lambda L: isinstance(L, list) and max(map(list_depth, L)) + 1

# Make geometries of landsat pathrows
pr_features = []
for pathrow, coords in pr_xw.items():
    coord_depth = list_depth(coords)
    if coord_depth == 3:
        geometry = {'type': 'Polygon', 'coordinates': coords}
    elif coord_depth == 4:
        geometry = {'type': 'MultiPolygon', 'coordinates': coords}

    pr_features.append(Feature(geometry=geometry))

pr_fc = FeatureCollection(pr_features)

# Make geometry of grid 8 tiles
# To go faster reduce bounds
tiles = list(mercantile.tiles(-180, -90, 180, 90, 8))
geoms = [mercantile.feature(tile)['geometry'] for tile in tiles]
gc = GeometryCollection(geoms)
feature = Feature(geometry=gc)
fc = FeatureCollection(features=[feature])

Visualize([fc, pr_fc])

sort by percent overlap

landsat-cogeo-mosaic/README.md

Lines 4 to 9 in 00386a7

	TODO: sort by percent overlap
	<!--
	landsat-cogeo-mosaic create \
	--bounds '-127.64,23.92,-64.82,52.72' \
	--max-cloud 05 \

---

This issue was generated by todo based on a TODO comment in 00386a7. It's been assigned to @kylebarron because they committed the code.

"Last-ditch" effort to find pathrow

Apparently, in some seasons there wasn't a single scene collected in some areas. For example, finding a mosaic of Winter 2015:

landsat-cogeo-mosaic create-from-db \
    --sqlite-path data/scene_list.db \
    --pathrow-index data/pr_index.json.gz \
    --min-zoom 7 \
    --max-zoom 12 \
    --min-date "2014-12-21" \
    --max-date "2015-03-21" \
    --max-cloud 100 \
    --sort-preference min-cloud \
    > test_winter.json

and visualizing it:

landsat-cogeo-mosaic visualize -p data/WRS2_descending_0/WRS2_descending.shp test_winter.json

you can see that some areas were never photographed

You should support a "last-ditch" retrieval for a pathrow. So for example, if there are no items found between min-date and max-date, then remove both and switch to nearest date of the average of max-date and min-date. That way you still have something for every mosaic for every pathrow in the index.

Compute quadkeys over land

Could be useful as a check to make sure that all quadkeys over land have assets

Add load_index_data function

Don't just have a index_data_path function, have a load_index_data function that optionally takes a path

landsat-cogeo-mosaic index fails with error

I'm trying to run

landsat-cogeo-mosaic index \
	    --wrs-path data/WRS2_descending.shp \
	    --scene-path data/scene_list.gz \
	    --bounds '-180,-85,180,85' \
	    --quadkey-zoom 8 \
	    > $@

Using the latest geopands, 0.8.1 it fails with the following error

Traceback (most recent call last):
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/bin/landsat-cogeo-mosaic", line 8, in <module>
    sys.exit(main())
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 829, in __call__
    return self.main(*args, **kwargs)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 782, in main
    rv = self.invoke(ctx)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 1259, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 1066, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 610, in invoke
    return callback(*args, **kwargs)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/landsat_cogeo_mosaic/cli.py", line 296, in index
    _index = create_index(
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/landsat_cogeo_mosaic/index.py", line 38, in create_index
    joined = gpd.sjoin(pathrows, tiles, op='intersects')
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/geopandas/tools/sjoin.py", line 56, in sjoin
    _crs_mismatch_warn(left_df, right_df, stacklevel=3)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/geopandas/array.py", line 93, in _crs_mismatch_warn
    left_srs = left.crs.to_string()
AttributeError: 'str' object has no attribute 'to_string'

Using geopandas 0.7.0 it fails with a different error

Traceback (most recent call last):
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/bin/landsat-cogeo-mosaic", line 8, in <module>
    sys.exit(main())
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 829, in __call__
    return self.main(*args, **kwargs)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 782, in main
    rv = self.invoke(ctx)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 1259, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 1066, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/click/core.py", line 610, in invoke
    return callback(*args, **kwargs)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/landsat_cogeo_mosaic/cli.py", line 296, in index
    _index = create_index(
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/landsat_cogeo_mosaic/index.py", line 41, in create_index
    gdf = optimize_index(joined)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/landsat_cogeo_mosaic/index.py", line 73, in optimize_index
    return pd.concat(res_df)
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/pandas/core/reshape/concat.py", line 274, in concat
    op = _Concatenator(
  File "/Users/jesse/venvs/landsat-cogeo-mosaic/lib/python3.8/site-packages/pandas/core/reshape/concat.py", line 331, in __init__
    raise ValueError("No objects to concatenate")
ValueError: No objects to concatenate

Both produces thousands of lines of

Shell is not a LinearRing
Shell is not a LinearRing
Shell is not a LinearRing
Shell is not a LinearRing

Script to upload/update dynamodb

try to use the new cogeo mosaic classes

Incorrect path-row positions

path-row positions are sometimes incorrect historically.

These are all defined with a row of 33 and a column of 35, but one of them is far off from the others:

So you should double check geometries of the assets when inserting. The mostly-overlapping set of assets have a maximum translation of 5km.

Missing tiles on east coast

Provide alternative sorts

For now, optimizing assets takes the one with most overlap first, but you could easily extend the sorting to use any landsat property, e.g. sun_azimuth, sun_elevation, etc

Missing areas

Missing data in corners?

Why do you sometimes get missing data in the corners? Here in montana...

https://kylebarron.dev/serverless-aerial-imagery/#9.67/47.9518/-113.851

Step through the data with features.geojson to try to reproduce.