diff --git a/rdwatch/core/tasks/__init__.py b/rdwatch/core/tasks/__init__.py index 12b07673..23f7c80d 100644 --- a/rdwatch/core/tasks/__init__.py +++ b/rdwatch/core/tasks/__init__.py @@ -331,7 +331,7 @@ def get_siteobservations_images( logger.info(f'Skipping Timestamp: {timestamp}') continue if found.exists() and not force: - found_timestamps[observation.timestamp] = True + found_timestamps[observation.timestamp.replace(microsecond=0)] = True continue results = fetch_boundbox_image( bbox, @@ -344,12 +344,12 @@ def get_siteobservations_images( logger.info(f'COULD NOT FIND ANY IMAGE FOR TIMESTAMP: {timestamp}') continue bytes = results['bytes'] - percent_black = get_percent_black_pixels(bytes) - cloudcover = results['cloudcover'] - found_timestamp = results['timestamp'] + found_timestamp = results['timestamp'].replace(microsecond=0) if bytes is None: logger.info(f'COULD NOT FIND ANY IMAGE FOR TIMESTAMP: {timestamp}') continue + percent_black = get_percent_black_pixels(bytes) + cloudcover = results['cloudcover'] if dayRange != -1 and percent_black < no_data_limit: found_timestamps[found_timestamp] = True elif dayRange == -1: @@ -523,6 +523,7 @@ def get_siteobservations_images( image_dimensions=[imageObj.width, imageObj.height], ) else: + logger.info('Skipping timestamp because image already found') count += 1 return downloaded_count diff --git a/rdwatch/core/utils/worldview_nitf/raster_tile.py b/rdwatch/core/utils/worldview_nitf/raster_tile.py index 21b756dd..638602f2 100644 --- a/rdwatch/core/utils/worldview_nitf/raster_tile.py +++ b/rdwatch/core/utils/worldview_nitf/raster_tile.py @@ -9,12 +9,39 @@ import boto3 import rasterio # type: ignore from rio_tiler.io.rasterio import Reader +from rio_tiler.models import ImageData +from rio_tiler.utils import pansharpening_brovey from rdwatch.core.utils.worldview_nitf.satellite_captures import WorldViewNITFCapture logger = logging.getLogger(__name__) +def find_rgb_channels(colorinterp): + rgb_indexes = {} + + # Define the channels to look for + channels = ['red', 'green', 'blue', 'gray'] + + for index, channel in enumerate(colorinterp): + if channel in channels: + rgb_indexes[channel] = index + + # Return indexes as a tuple, band indexes are 1 BASED + red = rgb_indexes.get('red', -1) + 1 + green = rgb_indexes.get('green', -1) + 1 + blue = rgb_indexes.get('blue', -1) + 1 + gray = rgb_indexes.get('gray', -1) + 1 + if red > 0 and green > 0 and blue > 0: + return (red, green, blue) + elif gray > 0: + return gray + else: + raise ValueError( + f'colorInterp: {colorinterp} for files does not have RGB or Gray channels' + ) + + @contextmanager def _download_nitf_image(uri: str) -> Generator[str, None, None]: s3 = boto3.client('s3') @@ -47,8 +74,10 @@ def get_worldview_nitf_tile( VSI_CACHE='TRUE', VSI_CACHE_SIZE=5000000, ): - with Reader(input=capture.uri) as img: - rgb = img.tile(x, y, z, tilesize=512) + with Reader(input=capture.uri) as vis_img: + information = vis_img.info() + rgb_channels = find_rgb_channels(information.colorinterp) + rgb = vis_img.tile(x, y, z, tilesize=512, indexes=rgb_channels) return rgb.render(img_format='WEBP') @@ -61,25 +90,87 @@ def get_worldview_nitf_bbox( with rasterio.Env( GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR', GDAL_HTTP_MERGE_CONSECUTIVE_RANGES='YES', - GDAL_CACHEMAX=200, - CPL_VSIL_CURL_CACHE_SIZE=20000000, - GDAL_BAND_BLOCK_CACHE='HASHSET', GDAL_HTTP_MULTIPLEX='YES', GDAL_HTTP_VERSION=2, - VSI_CACHE='TRUE', - VSI_CACHE_SIZE=5000000, + CPL_VSIL_CURL_CHUNK_SIZE=524288, ): - with _download_nitf_image(capture.uri) as nitf_filepath: - with Reader(input=nitf_filepath) as img: - rgb = img.part(bbox) - - if scale == 'default': - rgb.rescale(in_range=((0, 255),)) - elif scale == 'bits': - if capture.bits_per_pixel != 8: - max_bits = 2**capture.bits_per_pixel - 1 - rgb.rescale(in_range=((0, max_bits),)) - elif isinstance(scale, list) and len(scale) == 2: # scale is an integeter range - rgb.rescale(in_range=((scale[0], scale[1]),)) - - return rgb.render(img_format=format) + logger.info(f'Downloading WorldView NITF bbox: {bbox} scale: {scale}') + if capture.panuri is not None: + with Reader(input=capture.panuri) as pan_img: + logger.info(f'Downloading PanURI Chip: {capture.panuri}') + try: + pan_chip = pan_img.part(bbox=bbox, dst_crs='epsg:4326') + except ValueError as e: + logger.info(f'Value Error: {e} - {capture.panuri}') + return None + with Reader(input=capture.uri) as vis_img: + logger.info(f'Downloading URI Chip: {capture.uri}') + information = vis_img.info() + rgb_channels = find_rgb_channels(information.colorinterp) + if all( + channel is not None for channel in rgb_channels + ): # Ensure all RGB channels exist + try: + vis_chip = vis_img.part( + bbox=bbox, + indexes=rgb_channels, + dst_crs='epsg:4326', + width=pan_chip.width, + height=pan_chip.height, + ) + except ValueError as e: + logger.info(f'Value Error: {e} - {capture.uri}') + return None + final_chip = ImageData( + pansharpening_brovey( + vis_chip.data, pan_chip.data, 0.2, 'uint16' + ) + ) + statistics = list(final_chip.statistics().values()) + if scale == 'default': + in_range = tuple( + (item.min, item.max) for item in statistics + ) + elif scale == 'bits': + in_range = tuple( + (item['percentile_2'], item['percentile_98']) + for item in statistics + ) + elif ( + isinstance(scale, list) and len(scale) == 2 + ): # scale is an integeter range + in_range = tuple( + (scale[0], scale[1]) for item in statistics + ) + final_chip.rescale(in_range=in_range) + return final_chip.render(img_format=format) + else: + with Reader(input=capture.uri) as vis_img: + information = vis_img.info() + rgb_channels = find_rgb_channels(information.colorinterp) + if all( + channel is not None for channel in rgb_channels + ): # Ensure all RGB channels exist + try: + final_chip = vis_img.part( + bbox=bbox, + dst_crs='epsg:4326', + indexes=rgb_channels, + ) + except ValueError as e: + logger.info(f'Value Error: {e} - {capture.uri}') + return None + statistics = list(final_chip.statistics().values()) + if scale == 'default': + in_range = tuple((item.min, item.max) for item in statistics) + elif scale == 'bits': + in_range = tuple( + (item['percentile_2'], item['percentile_98']) + for item in statistics + ) + elif ( + isinstance(scale, list) and len(scale) == 2 + ): # scale is an integeter range + in_range = tuple((scale[0], scale[1]) for item in statistics) + final_chip.rescale(in_range=in_range) + return final_chip.render(img_format=format) diff --git a/rdwatch/core/utils/worldview_nitf/satellite_captures.py b/rdwatch/core/utils/worldview_nitf/satellite_captures.py index 88f22b6a..a5981f30 100644 --- a/rdwatch/core/utils/worldview_nitf/satellite_captures.py +++ b/rdwatch/core/utils/worldview_nitf/satellite_captures.py @@ -1,3 +1,4 @@ +import logging from dataclasses import dataclass from datetime import datetime, timedelta from typing import cast @@ -5,6 +6,8 @@ from rdwatch.core.utils.capture import URICapture from rdwatch.core.utils.worldview_nitf.stac_search import worldview_search +logger = logging.getLogger(__name__) + @dataclass() class WorldViewNITFCapture(URICapture): @@ -14,6 +17,8 @@ class WorldViewNITFCapture(URICapture): collection: str bits_per_pixel: int epsg: int + panuri: int | None + instruments: list[str] def get_features( @@ -34,13 +39,22 @@ def get_captures( timebuffer = timedelta(hours=1) features = [f for f in get_features(timestamp, bbox, timebuffer=timebuffer)] - captures = [] + vis_captures: list[WorldViewNITFCapture] = [] + pan_captures: list[WorldViewNITFCapture] = [] for feature in features: if 'data' in feature['assets']: cloudcover = 0 + instruments = [] if 'properties' in feature: if 'eo:cloud_cover' in feature['properties']: cloudcover = feature['properties']['eo:cloud_cover'] + if ( + 'nitf:compression' in feature['properties'] + and feature['properties']['nitf:compression'] != 'NC' + ): + continue + if 'instruments' in feature['properties']: + instruments = feature['properties']['instruments'] capture = WorldViewNITFCapture( timestamp=datetime.fromisoformat( feature['properties']['datetime'].rstrip('Z') @@ -50,8 +64,25 @@ def get_captures( bits_per_pixel=feature['properties']['nitf:bits_per_pixel'], cloudcover=cloudcover, collection=feature['collection'], + instruments=instruments, + panuri=None, epsg=feature['properties']['proj:epsg'], ) - captures.append(capture) + if 'panchromatic' in instruments: + pan_captures.append(capture) + elif 'vis-multi' in instruments: + vis_captures.append(capture) + else: + logger.info(f'Instrument type {instruments} is no supported') + # Attempt to add panuri elements to visual captures that exist + for pan_capture in pan_captures: + pan_timestamp = pan_capture.timestamp + for vis_capture in vis_captures: + vis_timestamp = vis_capture.timestamp + # Compare bbox and timestamp (within a tolerance for timestamp) + if abs((vis_timestamp - pan_timestamp).total_seconds()) < 60: + # If a match is found, update the vis-multi capture's panuri field + vis_capture.panuri = pan_capture.uri + break - return captures + return vis_captures diff --git a/rdwatch/scoring/tasks/__init__.py b/rdwatch/scoring/tasks/__init__.py index bc7c8acf..33008232 100644 --- a/rdwatch/scoring/tasks/__init__.py +++ b/rdwatch/scoring/tasks/__init__.py @@ -246,12 +246,12 @@ def get_siteobservations_images( logger.info(f'COULD NOT FIND ANY IMAGE FOR TIMESTAMP: {timestamp}') continue bytes = results['bytes'] - percent_black = get_percent_black_pixels(bytes) - cloudcover = results['cloudcover'] - found_timestamp = results['timestamp'] if bytes is None: logger.info(f'COULD NOT FIND ANY IMAGE FOR TIMESTAMP: {timestamp}') continue + percent_black = get_percent_black_pixels(bytes) + cloudcover = results['cloudcover'] + found_timestamp = results['timestamp'] if dayRange != -1 and percent_black < no_data_limit: found_timestamps[found_timestamp] = True elif dayRange == -1: diff --git a/rdwatch/scoring/tasks/animation_export.py b/rdwatch/scoring/tasks/animation_export.py index 253951e8..dfc0d7ba 100644 --- a/rdwatch/scoring/tasks/animation_export.py +++ b/rdwatch/scoring/tasks/animation_export.py @@ -601,24 +601,24 @@ def create_site_animation_export( celery_id=task_id, arguments=settings, ) - # try: - site_export.celery_id = task_id - site_export.save() - file_path, name = create_animation(self, site_evaluation_id, settings) - if file_path is False and name is False: - logger.warning('No Images were found') - site_export.delete() - return - site_export.name = name - with open(file_path, 'rb') as file: - site_export.export_file.save(name, File(file)) - site_export.save() - if os.path.exists(file_path): - os.remove(file_path) - # except Exception as e: - # logger.warning(f'Error when processing Animation: {e}') - # if site_export: - # site_export.delete() + try: + site_export.celery_id = task_id + site_export.save() + file_path, name = create_animation(self, site_evaluation_id, settings) + if file_path is False and name is False: + logger.warning('No Images were found') + site_export.delete() + return + site_export.name = name + with open(file_path, 'rb') as file: + site_export.export_file.save(name, File(file)) + site_export.save() + if os.path.exists(file_path): + os.remove(file_path) + except Exception as e: + logger.warning(f'Error when processing Animation: {e}') + if site_export: + site_export.delete() @app.task(bind=True)