WIP: Try to use cutline to remove broken edges at low zooms

Dockerfile

@@ -6,6 +6,8 @@ ENV PYTHONUSERBASE=/var/task
 
 COPY setup.py setup.py
 COPY landsat_mosaic_tiler/ landsat_mosaic_tiler/
+# Need to copy MANIFEST.in so that it correctly includes package data
+COPY MANIFEST.in MANIFEST.in
 
 RUN pip install . --user
 RUN rm -rf landsat_mosaic_tiler setup.py

MANIFEST.in

@@ -0,0 +1 @@
+include landsat_mosaic_tiler/data/wrs.db

landsat_mosaic_tiler/cutline.py

@@ -0,0 +1,49 @@
+"""landsat_mosaic_tiler.cutline: Generate cutline for image."""
+
+from rasterio.crs import CRS
+from rasterio.warp import transform_geom
+from shapely import wkt
+from shapely.geometry import Polygon, box, mapping, shape
+from shapely.ops import transform
+
+# import rasterio
+# src = 's3://landsat-pds/c1/L8/037/034/LC08_L1TP_037034_20200707_20200707_01_RT/LC08_L1TP_037034_20200707_20200707_01_RT_B1.TIF'
+# r = rasterio.open(src)
+# pathrow = '0'
+# sceneid = 'LC08_L1TP_037034_20200707_20200707_01_RT'
+# r.bounds
+# transform_bounds
+# Visualize(box(*transform_bounds(r.crs, CRS.from_epsg(4326), *r.bounds)))
+
+
+def reproject_geom(geom, from_crs, to_crs):
+ return shape(transform_geom(from_crs, to_crs, mapping(geom)))
+
+
+def to_local_coords(geom, r):
+ return reproject_geom(geom, from_crs=CRS.from_epsg(4326), to_crs=r.crs)
+
+
+def to_image_coords(geom, r):
+ return transform(r.index, geom)
+
+
+def get_cutline(r, geom: Polygon):
+ """Get cutline to invalid edge pixels from image
+
+ Cutline is in **image** coordinates.
+
+ Args:
+ - r: opened rasterio dataset
+ - geom: Polygon in WGS84 of valid part of the image
+ """
+ # Convert geom to source's local CRS
+ geom_local = to_local_coords(geom, r)
+
+ # Convert geom to image coordinates
+ geom_image_coords = to_image_coords(geom_local, r)
+
+ # Intersect with image
+ image_bounds = box(0, 0, r.width, r.height)
+ cutline = geom_image_coords.intersection(image_bounds)
+ return wkt.dumps(cutline)

landsat_mosaic_tiler/handlers/tiles.py

@@ -6,18 +6,20 @@
 
 import mercantile
 import numpy
+from lambda_proxy.proxy import API
+from landsat_mosaic_tiler.landsat_tiler import tile as landsatTiler
 from landsat_mosaic_tiler.pixel_methods import pixSel
 from landsat_mosaic_tiler.utils import get_tilejson, post_process_tile
-from cogeo_mosaic.backends import MosaicBackend
-from lambda_proxy.proxy import API
 from PIL import Image
 from rasterio.transform import from_bounds
 from rio_tiler.colormap import get_colormap
-from rio_tiler.io.landsat8 import tile as landsatTiler
 from rio_tiler.profiles import img_profiles
 from rio_tiler.utils import expression as expressionTiler
 from rio_tiler.utils import render
 from rio_tiler_mosaic.mosaic import mosaic_tiler
+# from rio_tiler.io.landsat8 import tile as origLandsatTiler
+
+from cogeo_mosaic.backends import MosaicBackend
 
 app = API(name="landsat-mosaic-tiler-tiles", debug=False)
 
@@ -237,6 +239,8 @@ def tiles(
  sio.seek(0)
  return ("OK", f"image/{ext}", sio.getvalue(), return_kwargs)
 
+ # tile, mask = origLandsatTiler(sceneid, x, y, z, bands)
+ # tile, mask = landsatTiler(sceneid, x, y, z, bands)
  rtile = post_process_tile(tile, mask, rescale=rescale, color_formula=color_ops)
 
  if ext == "bin":
@@ -246,7 +250,8 @@ def tiles(
 
  driver = "jpeg" if ext == "jpg" else ext
  options = img_profiles.get(driver, {})
-
+ # with open('test2.png', 'wb') as f:
+ # f.write(render(rtile, mask, img_format=driver, colormap=color_map, **options))
  if ext == "tif":
  ext = "tiff"
  driver = "GTiff"

landsat_mosaic_tiler/landsat_tiler.py

@@ -0,0 +1,175 @@
+"""
+landsat_mosaic_tiler:landsat_tiler.py
+
+Custom landsat tiler to use cutline to remove broken overview edges
+"""
+
+import os
+import sqlite3
+from concurrent import futures
+from pathlib import Path
+from typing import Any, Dict, Sequence, Tuple, Union
+
+import numpy
+import rasterio
+from landsat_mosaic_tiler.cutline import get_cutline
+from rio_tiler import constants, reader
+from rio_tiler.errors import InvalidBandName, TileOutsideBounds
+from rio_tiler.io.landsat8 import (LANDSAT_BANDS, _convert, _landsat_get_mtl,
+ landsat_parser)
+from rio_tiler.utils import pansharpening_brovey, tile_exists
+from rio_toa import toa_utils
+from shapely import wkb
+
+
+def find_wrs_db_path():
+ lambda_root = os.getenv('LAMBDA_TASK_ROOT')
+ if not lambda_root:
+ return str((Path(__file__) / 'data' / 'wrs.db').resolve())
+
+ return f'{lambda_root}/landsat_mosaic_tiler/data/wrs.db'
+
+
+def get_wrs_geometry(db_path, pathrow):
+ """Get scene geometry from SQLite DB of WRS geometries
+ """
+ with sqlite3.connect(db_path) as conn:
+ c = conn.cursor()
+ t = (pathrow, )
+ c.execute('SELECT geometry FROM wrs2 WHERE pathrow=?', t)
+ wkb_geom = c.fetchone()[0]
+ return wkb.loads(wkb_geom)
+
+
+# ['LC08_L1TP_036035_20190628_20190706_01_T1',
+# 'LC08_L1TP_036034_20190815_20190821_01_T1',
+# 'LC08_L1TP_037034_20190721_20190801_01_T1',
+# 'LC08_L1TP_037035_20190705_20190719_01_T1']
+#
+# sceneid = 'LC08_L1TP_036035_20190628_20190706_01_T1'
+# tile_x = x
+# tile_y = y
+# tile_z = z
+# bands
+# tilesize
+# pan
+# kwargs= {}
+# wrs_db_path = 'data/wrs.db'
+
+def tile(
+ sceneid: str,
+ tile_x: int,
+ tile_y: int,
+ tile_z: int,
+ bands: Union[Sequence[str], str] = ["4", "3", "2"],
+ tilesize: int = 256,
+ pan: bool = False,
+ **kwargs: Any,
+) -> Tuple[numpy.ndarray, numpy.ndarray]:
+ """
+ Create mercator tile from Landsat-8 data.
+
+ Attributes
+ ----------
+ sceneid : str
+ Landsat sceneid. For scenes after May 2017,
+ sceneid have to be LANDSAT_PRODUCT_ID.
+ tile_x : int
+ Mercator tile X index.
+ tile_y : int
+ Mercator tile Y index.
+ tile_z : int
+ Mercator tile ZOOM level.
+ bands : tuple, str, optional (default: ("4", "3", "2"))
+ Bands index for the RGB combination.
+ tilesize : int, optional (default: 256)
+ Output image size.
+ pan : boolean, optional (default: False)
+ If True, apply pan-sharpening.
+ kwargs: dict, optional
+ These will be passed to the 'rio_tiler.utils._tile_read' function.
+
+ Returns
+ -------
+ data : numpy ndarray
+ mask: numpy array
+
+ """
+ if isinstance(bands, str):
+ bands = (bands, )
+
+ for band in bands:
+ if band not in LANDSAT_BANDS:
+ raise InvalidBandName(
+ "{} is not a valid Landsat band name".format(band))
+
+ scene_params = landsat_parser(sceneid)
+
+ meta: Dict = _landsat_get_mtl(sceneid)["L1_METADATA_FILE"]
+
+ landsat_prefix = "{scheme}://{bucket}/{prefix}/{scene}".format(
+ **scene_params)
+
+ bounds = toa_utils._get_bounds_from_metadata(meta["PRODUCT_METADATA"])
+ if not tile_exists(bounds, tile_z, tile_x, tile_y):
+ raise TileOutsideBounds(
+ "Tile {}/{}/{} is outside image bounds".format(
+ tile_z, tile_x, tile_y))
+
+ # Find geometry from wrs2 grid, to later create cutline
+ wrs_db_path = find_wrs_db_path()
+ pathrow = scene_params['path'] + scene_params['row']
+ geom = get_wrs_geometry(wrs_db_path, pathrow)
+
+ def worker(band: str):
+ asset = f"{landsat_prefix}_B{band}.TIF"
+
+ if band == "QA":
+ nodata = 1
+ resamp = "nearest"
+ else:
+ nodata = 0
+ resamp = "bilinear"
+
+ with rasterio.open(asset) as src_dst:
+ # Create cutline
+ # The cutline removes broken pixels along the edges of the overviews
+ cutline = get_cutline(src_dst, geom)
+
+ tile, mask = reader.tile(
+ src_dst,
+ tile_x,
+ tile_y,
+ tile_z,
+ tilesize=tilesize,
+ nodata=nodata,
+ resampling_method=resamp,
+ warp_vrt_option={'cutline': cutline},
+ **kwargs)
+
+ return tile, mask
+
+ results = [worker(band) for band in bands]
+ data = [x[0] for x in results]
+ masks = [x[1] for x in results]
+ data = numpy.concatenate(data)
+ mask = numpy.all(masks, axis=0).astype(numpy.uint8) * 255
+ return data, mask
+
+
+
+ with futures.ThreadPoolExecutor(
+ max_workers=constants.MAX_THREADS) as executor:
+ data, masks = zip(*list(executor.map(worker, bands)))
+ data = numpy.concatenate(data)
+ mask = numpy.all(masks, axis=0).astype(numpy.uint8) * 255
+
+ if pan:
+ pan_data, mask = worker("8")
+ data = pansharpening_brovey(data, pan_data, 0.2, pan_data.dtype)
+
+ if bands[0] != "QA" or len(bands) != 1:
+ for bdx, band in enumerate(bands):
+ data[bdx] = _convert(data[bdx], band, meta)
+
+ return data, mask

setup.py

@@ -13,8 +13,11 @@
  "loguru",
  "mercantile",
  "rio-color",
- "rio-tiler>=2.0a6",
  "rio_tiler_mosaic>=0.0.1dev3",
+ "rio-tiler==2.0a11",
+ # To turn off threading
+ "loguru",
+ "shapely"
 ]
 extra_reqs = {
  "test": ["pytest", "pytest-cov", "mock"],