diff --git a/Dockerfile b/Dockerfile index ade5f0f..3f78308 100644 --- a/Dockerfile +++ b/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 \ No newline at end of file diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..2d7b798 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1 @@ +include landsat_mosaic_tiler/data/wrs.db diff --git a/landsat_mosaic_tiler/cutline.py b/landsat_mosaic_tiler/cutline.py new file mode 100644 index 0000000..7e8233d --- /dev/null +++ b/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) diff --git a/landsat_mosaic_tiler/data/wrs.db b/landsat_mosaic_tiler/data/wrs.db new file mode 100644 index 0000000..d1d54a0 Binary files /dev/null and b/landsat_mosaic_tiler/data/wrs.db differ diff --git a/landsat_mosaic_tiler/handlers/tiles.py b/landsat_mosaic_tiler/handlers/tiles.py index f6c68a4..5ae073e 100644 --- a/landsat_mosaic_tiler/handlers/tiles.py +++ b/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" diff --git a/landsat_mosaic_tiler/landsat_tiler.py b/landsat_mosaic_tiler/landsat_tiler.py new file mode 100644 index 0000000..bdbbf0a --- /dev/null +++ b/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 diff --git a/landsat_mosaic_tiler/test.png b/landsat_mosaic_tiler/test.png new file mode 100644 index 0000000..914a7a6 Binary files /dev/null and b/landsat_mosaic_tiler/test.png differ diff --git a/landsat_mosaic_tiler/test2.png b/landsat_mosaic_tiler/test2.png new file mode 100644 index 0000000..2d6c94d Binary files /dev/null and b/landsat_mosaic_tiler/test2.png differ diff --git a/setup.py b/setup.py index 47ca44a..17ea086 100644 --- a/setup.py +++ b/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"], diff --git a/test.png b/test.png new file mode 100644 index 0000000..aa1a8f4 Binary files /dev/null and b/test.png differ diff --git a/test2.png b/test2.png new file mode 100644 index 0000000..348bac5 Binary files /dev/null and b/test2.png differ