feat: add preliminary map tile creation

setup.cfg

@@ -41,6 +41,7 @@ include_package_data = True
 install_requires =
     numpy>=1.23.0
     scikit-optimize>=0.9.0
+    opencv-python-headless>=4.9.0
     cachetools>=5.3.2
     geojson>=3.1.0
     pyproj>=3.6.1

src/aws/osml/gdal/dynamic_range_adjustment.py

@@ -72,12 +72,12 @@ def from_counts(
         max_counts = cumulative_counts[-1]
         low_threshold = min_percentage * max_counts
         e_min = 0
-        while cumulative_counts[e_min] < low_threshold:
+        while cumulative_counts[e_min] < low_threshold and e_min < len(cumulative_counts) - 1:
             e_min += 1
 
         high_threshold = max_percentage * max_counts
         e_max = num_histogram_bins - 1
-        while cumulative_counts[e_max] > high_threshold:
+        while cumulative_counts[e_max] > high_threshold and e_max > 0:
             e_max -= 1
 
         min_value = max([actual_min_value, e_min - a * (e_max - e_min)])

src/aws/osml/image_processing/__init__.py

@@ -50,6 +50,35 @@
 
     viz_tile = viz_tile_factory.create_encoded_tile([0, 0, 1024, 1024], output_size=(512, 512))
 
+Image Tiling: Map Tiles / Orthophotos
+*************************************
+
+The TileFactory supports creation of tiles suitable for use by geographic information systems (GIS) or map-based
+visualization tools. Given a north-east aligned bounding box in geographic coordinates the tile factory can use
+the sensor models to orthorectify imagery to remove the perspective and terrain effects.
+
+.. code-block:: python
+    :caption: Example showing creation of a map tile from the WebMercatorQuad tile set
+
+    # Look up the tile boundary for a tile in a well known tile set
+    tile_set = MapTileSetFactory.get_for_id("WebMercatorQuad")
+    tile_id = MapTileId(tile_matrix=16, tile_row=37025, tile_col=54816)
+    tile = tile_set.get_tile(tile_id)
+
+    # Create an orthophoto for this tile
+    image_bytes = viz_tile_factory.create_orthophoto_tile(geo_bbox=tile.bounds, tile_size=tile.size)
+
+.. figure:: ../images/MapTileExample-BeforeAfter.png
+    :width: 600
+    :alt: Original image with perspective effects and same area after orthorectification
+
+    Example showing original image with perspective effects and same area after orthorectification.
+
+.. figure:: ../images/MapTileExample-MapOverlay.png
+    :width: 400
+    :alt: Orthophoto tile overlaid on Google Maps
+
+    Example showing map tile overlaid on Google Maps
 
 Complex SAR Data Display
 ************************
@@ -93,6 +122,17 @@
 """
 
 from .gdal_tile_factory import GDALTileFactory
+from .map_tileset import MapTile, MapTileId, MapTileSet
+from .map_tileset_factory import MapTileSetFactory, WellKnownMapTileSet
 from .sar_complex_imageop import histogram_stretch, quarter_power_image
 
-__all__ = ["GDALTileFactory", "histogram_stretch", "quarter_power_image"]
+__all__ = [
+    "GDALTileFactory",
+    "MapTile",
+    "MapTileId",
+    "MapTileSet",
+    "MapTileSetFactory",
+    "WellKnownMapTileSet",
+    "histogram_stretch",
+    "quarter_power_image",
+]

src/aws/osml/image_processing/gdal_tile_factory.py

@@ -5,10 +5,14 @@
 from secrets import token_hex
 from typing import Any, Dict, List, Optional, Tuple
 
+import cv2
+import numpy as np
 from osgeo import gdal, gdalconst
+from scipy.interpolate import RectBivariateSpline
 
 from aws.osml.gdal import GDALCompressionOptions, GDALImageFormats, NITFDESAccessor, RangeAdjustmentType, get_type_and_scales
-from aws.osml.photogrammetry import ImageCoordinate, SensorModel
+from aws.osml.gdal.dynamic_range_adjustment import DRAParameters
+from aws.osml.photogrammetry import GeodeticWorldCoordinate, ImageCoordinate, SensorModel
 
 from .sicd_updater import SICDUpdater
 from .sidd_updater import SIDDUpdater
@@ -95,7 +99,7 @@ def create_encoded_tile(
         # Create a new IGEOLO value based on the corner points of this tile
         if self.sensor_model is not None and self.tile_format == GDALImageFormats.NITF:
             gdal_translate_kwargs["creationOptions"].append("ICORDS=G")
-            gdal_translate_kwargs["creationOptions"].append("IGEOLO=" + self.create_new_igeolo(src_window))
+            gdal_translate_kwargs["creationOptions"].append("IGEOLO=" + self._create_new_igeolo(src_window))
 
         if self.sar_updater is not None and self.tile_format == GDALImageFormats.NITF:
             # If we're outputting a SICD or SIDD tile we need to update the XML metadata to include the new chip
@@ -135,7 +139,199 @@ def create_encoded_tile(
                 gdal.VSIFCloseL(vsifile_handle)
             gdal.GetDriverByName(self.tile_format).Delete(temp_ds_name)
 
-    def create_new_igeolo(self, src_window: List[int]) -> str:
+    def create_orthophoto_tile(
+        self, geo_bbox: Tuple[float, float, float, float], tile_size: Tuple[int, int]
+    ) -> Optional[bytearray]:
+        """
+        This method creates an orthorectified tile from an image assuming there is overlap in the coverage.
+
+        IMPORTANT: This is an experimental API that may change in future minor releases of the toolkit. This
+        early release is subject to the following limitations:
+        - All tiles are returned in PNG format
+        - An 8-bit conversion and dynamic range mapping is automatically applied using the statistics of the tile
+
+        :param geo_bbox: the geographic bounding box of the tile in the form (min_lon, min_lat, max_lon, max_lat)
+        :param tile_size: the shape of the output tile (width, height)
+        :return: the encoded image tile or None if one could not be produced
+        """
+        min_lon, min_lat, max_lon, max_lat = geo_bbox
+
+        # Setup 2 grids of the same size, the first is for longitude/latitude coordinates and the second is pixel
+        # coordinates. These grids are evenly spaced across the map tile. Note that the latitude and pixel row
+        # grids have been adjusted because the 0, 0 pixel is in the upper left corner of the map tile and as the
+        # image row increases the latitude should decrease. That is why the world_y grid ranges from max to min
+        # while all other grids range from min to max.
+        nx, ny = (3, 3)
+        world_x = np.linspace(min_lon, max_lon, nx)
+        world_y = np.linspace(max_lat, min_lat, ny)
+        world_xv, world_yv = np.meshgrid(world_x, world_y)
+
+        pixel_x = np.linspace(0, tile_size[0] - 1, nx)
+        pixel_y = np.linspace(0, tile_size[1] - 1, ny)
+
+        # Use the sensor model to compute the image pixel location that corresponds to each
+        # world coordinate in the map tile grid. Separate those results into arrays of the x, y
+        # component.
+        def world_to_image_func(lon, lat):
+            # TODO: Assign the elevation from a DEM
+            return self.sensor_model.world_to_image(GeodeticWorldCoordinate([lon, lat, 0.0]))
+
+        world_to_image_func_vectorized = np.vectorize(world_to_image_func)
+        src_coords = world_to_image_func_vectorized(world_xv, world_yv)
+        src_x = np.vectorize(lambda image_coord: image_coord.x)(src_coords)
+        src_y = np.vectorize(lambda image_coord: image_coord.y)(src_coords)
+
+        # Find min/max x and y for this grid and check to make sure it actually overlaps the image.
+        src_bbox = (
+            int(np.floor(np.min(src_x))),
+            int(np.floor(np.min(src_y))),
+            int(np.ceil(np.max(src_x))),
+            int(np.ceil(np.max(src_y))),
+        )
+        if (
+            src_bbox[0] > self.raster_dataset.RasterXSize
+            or src_bbox[1] > self.raster_dataset.RasterYSize
+            or src_bbox[2] < 0
+            or src_bbox[3] < 0
+        ):
+            # Source bbox does not intersect the image, no tile to create
+            return None
+
+        # Clip the source to the extent of the image and then select an overview level of similar resolution to the
+        # desired map tile. This will ensure we only read the minimum number of pixels necessary and warp them as
+        # little as possible.
+        src_bbox = (
+            max(src_bbox[0], 0),
+            max(src_bbox[1], 0),
+            min(src_bbox[2], self.raster_dataset.RasterXSize),
+            min(src_bbox[3], self.raster_dataset.RasterYSize),
+        )
+        logger.debug(f"After Clip to Image Bounds src_bbox = {src_bbox}")
+
+        def find_appropriate_r_level(src_bbox, tile_width) -> int:
+            src_dim = np.min([src_bbox[2] - src_bbox[0], src_bbox[3] - src_bbox[1]])
+            return int(np.max([0, int(np.floor(np.log2(src_dim / tile_width)))]))
+
+        num_overviews = self.raster_dataset.GetRasterBand(1).GetOverviewCount()
+        r_level = min(find_appropriate_r_level(src_bbox, tile_size[0]), num_overviews)
+        overview_bbox = tuple([int(x / 2**r_level) for x in src_bbox])
+        logger.debug(f"Using r-level: {r_level}")
+        logger.debug(f"overview_bbox = {overview_bbox}")
+        logger.debug(f"Dataset size = {self.raster_dataset.RasterXSize},{self.raster_dataset.RasterYSize}")
+
+        # Read pixels from the selected resolution level that match the region of the image needed to create the
+        # map tile. This data becomes the "src" in the cv2.remap transformation.
+        src = self._read_from_rlevel_as_array(overview_bbox, r_level)
+        logger.debug(f"src.shape = {src.shape}")
+
+        # Update the src_x and src_y coordinates because we cropped the image and pulled it from a different
+        # resolution level. The original coordinates assumed the image origin at 0,0 in a full resolution
+        # image
+        src_x = (src_x - src_bbox[0]) / 2**r_level
+        src_y = (src_y - src_bbox[1]) / 2**r_level
+
+        # Create 2D linear interpolators that map the pixels in the map tile to x and y values in the source image.
+        # This will allow us to efficiently generate the maps needed by the opencv::remap function for every pixel
+        # in the destination image.
+        src_x_interpolator = RectBivariateSpline(pixel_x, pixel_y, src_x, kx=1, ky=1)
+        src_y_interpolator = RectBivariateSpline(pixel_x, pixel_y, src_y, kx=1, ky=1)
+
+        # Create the map1 and map2 arrays that capture the non-linear relationship between each pixel in the map tile
+        # (dst) to pixels in the original image (src). See opencv::remap documentation for definitions of these
+        # parameters.
+        dst_x = np.linspace(0, tile_size[0] - 1, tile_size[0])
+        dst_y = np.linspace(0, tile_size[1] - 1, tile_size[1])
+        map1 = src_x_interpolator(dst_x, dst_y).astype(np.float32)
+        map2 = src_y_interpolator(dst_x, dst_y).astype(np.float32)
+
+        logger.debug(
+            f"Sanity check remap array sizes. " f"They should match the desired map tile size {tile_size[0]}x{tile_size[1]}"
+        )
+        logger.debug(f"map1.shape = {map1.shape}")
+        logger.debug(f"map2.shape = {map2.shape}")
+
+        dst = cv2.remap(src, map1, map2, cv2.INTER_LINEAR)
+        output_tile_pixels = self._create_display_image(dst)
+
+        # TODO: Formats other than PNG?
+        is_success, image_bytes = cv2.imencode(".png", output_tile_pixels)
+        if is_success:
+            return image_bytes
+        else:
+            return None
+
+    def _read_from_rlevel_as_array(
+        self, scaled_bbox: Tuple[int, int, int, int], r_level: int, band_numbers: Optional[List[int]] = None
+    ) -> np.array:
+        """
+        This method reads a region of the image from a specific image resolution level (r-level). The bounding box
+        must be scaled to match the resolution level.
+
+        :param scaled_bbox: a [minx, miny, maxx, maxy] bbox in the pixel coordinate system of the r_level
+        :param r_level: the selected resolution level, r0 = full resolution image, r1 = first overview, ...
+        :param band_numbers: the bands to select, if None all bands will be read. Note band numbers start at 1
+        :return: a NumPy array of shape [r, c, b] for images with multiple bands or [r, c] for images with just 1 band
+        """
+
+        # If no bands are specified we will read them all.
+        if not band_numbers:
+            band_numbers = [n + 1 for n in range(0, self.raster_dataset.RasterCount)]
+
+        # Loop through the bands of this image and retrieve the relevant pixels
+        band_pixels = []
+        for band_num in band_numbers:
+            ds_band = self.raster_dataset.GetRasterBand(band_num)
+            if r_level > 0:
+                overview = ds_band.GetOverview(r_level - 1)
+            else:
+                overview = ds_band
+
+            band_pixels.append(
+                overview.ReadAsArray(
+                    scaled_bbox[0],
+                    scaled_bbox[1],
+                    scaled_bbox[2] - scaled_bbox[0],
+                    scaled_bbox[3] - scaled_bbox[1],
+                )
+            )
+
+        # If the image has multiple bands then we can stack the results with the band being the 3rd dimension
+        # in the array. This aligns to how OpenCV wants to work with imagery. If the image doesn't have multiple
+        # bands then return a 2-dimensional grayscale image.
+        if self.raster_dataset.RasterCount > 1:
+            result = np.stack(band_pixels, axis=2)
+        else:
+            result = band_pixels[0]
+
+        return result
+
+    def _create_display_image(self, pixel_array: np.array) -> np.array:
+        """
+        This method selects the first 3 bands of a multi-band image (or a single band of grayscale) and performs a
+        simple dynamic range adjustment to map those values to the 0-255 range of an 8-bit per pixel image for
+        visualization.
+
+        :param pixel_array: the input image pixels
+        :return: a range adjusted 8-bit per pixel image of up to 3 bands
+        """
+        if pixel_array.ndim == 3 and pixel_array.shape[2] > 3:
+            pixel_array = pixel_array[:, :, 0:3]
+
+        max_channel_value = np.max(pixel_array)
+        hist, bin_edges = np.histogram(pixel_array.ravel(), bins=max_channel_value + 1, range=(0, max_channel_value))
+        dra_parameters = DRAParameters.from_counts(hist.tolist(), max_percentage=0.97)
+        for_display = (
+            255
+            * (pixel_array - dra_parameters.suggested_min_value)
+            / (dra_parameters.suggested_max_value - dra_parameters.suggested_min_value)
+        )
+
+        for_display = np.clip(for_display, 0.0, 255.0)
+        for_display = for_display.astype(np.uint8)
+
+        return for_display
+
+    def _create_new_igeolo(self, src_window: List[int]) -> str:
         """
         Create a new 60 character string representing the corner coordinates of this tile. The string conforms to
         the geographic "G" choice for the ICORDS field as specified in MIL-STD-2500C Table A-3 NITF Image Subheader.

src/aws/osml/image_processing/map_tileset.py

@@ -0,0 +1,87 @@
+#  Copyright 2024 Amazon.com, Inc. or its affiliates.
+
+from abc import ABC, abstractmethod
+from collections import namedtuple
+from dataclasses import dataclass
+
+from aws.osml.photogrammetry import GeodeticWorldCoordinate
+
+MapTileId = namedtuple("TildId", "tile_matrix tile_row tile_col")
+MapTileId.__doc__ = """
+This type represents the unique id of a map tile within a map tile set. It is implemented as a named tuple to make
+use of that constructs immutability and hashing features.
+"""
+
+MapTileSize = namedtuple("MapTileSize", "width height")
+MapTileSize.__doc__ = """
+This type represents the size of a map tile (width, height). These dimensions are typically the same (i.e. square
+tiles) but this is not required.
+"""
+
+MapTileBounds = namedtuple("MapTileBounds", "min_lon min_lat max_lon max_lat")
+MapTileBounds.__doc__ = """
+This type represents the geodetic bounds of a map tile (min_lon, min_lat, max_lon, max_lat).
+"""
+
+
+@dataclass
+class MapTile:
+    """
+    This dataclass provides a description of a map tile that is part of a well known tile set.
+    """
+
+    id: MapTileId
+    size: MapTileSize
+    bounds: MapTileBounds
+
+
+class MapTileSet(ABC):
+    """
+    This class provides an abstract interface to a well known set of map tiles.
+    """
+
+    @property
+    @abstractmethod
+    def tile_matrix_set_id(self) -> str:
+        """
+        Get the identifier for this map tile set. This is the tile matrix set ID in the OGC definitions.
+
+        :return: the tile matrix set ID
+        """
+
+    @abstractmethod
+    def get_tile(self, tile_id: MapTileId) -> MapTile:
+        """
+        Get a description of the tile identified by a specific map tile ID.
+
+        :param tile_id: the tile ID
+        :return: the tile description
+        """
+
+    @abstractmethod
+    def get_tile_for_location(self, world_coordinate: GeodeticWorldCoordinate, tile_matrix: int) -> MapTile:
+        """
+        Get a description of the tile containing a specific world location.
+
+        :param world_coordinate: the location in the world
+        :param tile_matrix: the tile_matrix or zoom level of interest
+        :return: the tile description
+        """
+
+    def get_tile_matrix_limits_for_area(
+        self, boundary_coordinates: list[GeodeticWorldCoordinate], tile_matrix: int
+    ) -> tuple[int, int, int, int]:
+        """
+        Get a list of all tiles that intersect a specific area.
+
+        :param boundary_coordinates: the boundary of the area
+        :param tile_matrix: the tile_matrix or zoom level of interest
+        :return: the (min_col, min_row, max_col, max_row) limits of tiles containing all points
+        """
+        map_tiles_for_corners = [
+            self.get_tile_for_location(world_corner, tile_matrix=tile_matrix) for world_corner in boundary_coordinates
+        ]
+        tile_rows = [tile_id.id.tile_row for tile_id in map_tiles_for_corners]
+        tile_cols = [tile_id.id.tile_col for tile_id in map_tiles_for_corners]
+
+        return min(tile_cols), min(tile_rows), max(tile_cols), max(tile_rows)