diff --git a/README.md b/README.md index 5f159aa..f46c9bf 100644 --- a/README.md +++ b/README.md @@ -2,14 +2,24 @@ The OversightML Imagery Toolkit is a Python package that contains image processing and photogrammetry routines commonly used during the analysis of imagery collected by satellites and unmanned aerial vehicles (UAVs). It builds upon GDAL -by providing additional support for images compliant with the National Imagery Transmission Format (NITF) and Sensor -Independent Complex Data (SICD) standards. +by providing additional support for images compliant with the National Imagery Transmission Format (NITF), Sensor +Independent Complex Data (SICD), and Sensor Independent Derived Data (SIDD) standards. +This library contains four core packages under the `aws.osml` namespace: +* **photogrammetry**: convert locations between the image (x, y) and geodetic (lon, lat, elev) coordinate systems +* **gdal**: utilities to work with datasets loaded by GDAL +* **image_processing**: common image manipulation routines +* **features**: common geospatial feature manipulation routines + +## Documentation + +* **APIs**: You can find API documentation for the OSML Imagery Toolkit hosted on our [GitHub project page](https://aws-solutions-library-samples.github.io/osml-imagery-toolkit/). +If you are working from the source code running `tox -e docs` will trigger the Sphinx documentation build. +* **Example Notebooks**: Example notebooks for some operations are in the `examples` directory ## Installation -The intent is to vend / distribute this software through a Python Package Index. -If your environment has a distribution, -you should be able to install it using pip: +This software is available through a Python Package Index. +If your environment has a distribution, you should be able to install it using pip: ```shell pip install osml-imagery-toolkit[gdal] ``` @@ -23,97 +33,6 @@ Note that GDAL is listed as an extra dependency for this package. This is done t don't want to use GDAL or those that have their own custom installation steps for that library. Future versions of this package will include image IO backbones that have fewer dependencies. - -## Documentation - -You can find documentation for this library in the `./doc` directory. Sphinx is used to construct a searchable HTML -version of the API documents. - -```shell -tox -e docs -``` - -## Example Usage - -This library contains four core packages under the `aws.osml` namespace. -* photogrammetry: convert locations between the image (x, y) and geodetic (lon, lat, elev) coordinate systems -* gdal: help load and manage datasets loaded by GDAL -* image_processing: common image manipulation routines -* formats: utilities for handling format specific information; normally not accessed directly - -```python -from aws.osml.gdal import GDALImageFormats, GDALCompressionOptions, load_gdal_dataset -from aws.osml.image_processing import GDALTileFactory -from aws.osml.photogrammetry import ImageCoordinate, GeodeticWorldCoordinate, SensorModel -``` - -### Tiling with Updated Image Metadata - -Many applications break large remote sensing images into smaller chips or tiles for distributed processing or -dissemination. GDAL's Translate function provides basic capabilities, but it does not correctly update geospatial -metadata to reflect the new image extent. These utilities provide those functions so tile consumers can correctly -interpret the pixel information they have been provided. For NITF imagery that includes the addition of a new ICHIPB -TRE. With SICD the XML ImageData elements are adjusted to identify the sub-image bounds. - -```python -# Load the image and create a sensor model -ds, sensor_model = load_gdal_dataset("./imagery/sample.nitf") -tile_factory = GDALTileFactory(ds, - sensor_model, - GDALImageFormats.NITF, - GDALCompressionOptions.NONE - ) - -# Bounds are [left_x, top_y, width, height] -nitf_encoded_tile_bytes = tile_factory.create_encoded_tile([0, 0, 1024, 1024]) -``` - -### Tiling for Display - -Some images, for example 11-bit panchromatic images or SAR imagery with floating point complex data, can not be -displayed directly without remapping the pixels into an 8-bit per pixel grayscale or RGB color model. The TileFactory -supports creation of tiles suitable for human review by setting both the output_type and range_adjustment options. - -```python -viz_tile_factory = GDALTileFactory(ds, - sensor_model, - GDALImageFormats.PNG, - GDALCompressionOptions.NONE, - output_type=gdalconst.GDT_Byte, - range_adjustment=RangeAdjustmentType.DRA) - -viz_tile = viz_tile_factory.create_encoded_tile([0, 0, 1024, 1024]) -``` - -### More Precise Sensor Models - -OversightML provides implementations of the Replacement Sensor Model (RSM), Rational Polynomial -Camera (RPC), and Sensor Independent Complex Data (SICD) sensor models to assist in geo positioning. -When loading a dataset, the toolkit will construct the most accurate sensor model -from the available image metadata. That sensor model can be used in conjunction with an optional -elevation model to convert between image and geodetic coordinates. - -```python -ds, sensor_model = load_gdal_dataset("./imagery/sample.nitf") -elevation_model = DigitalElevationModel( - SRTMTileSet(version="1arc_v3"), - GDALDigitalElevationModelTileFactory("./local-SRTM-tiles")) - -# Note the order of ImageCoordinate is (x, y) -geodetic_location_of_ul_corner = sensor_model.image_to_world( - ImageCoordinate([0, 0]), - elevation_model=elevation_model) - -lon_degrees = -77.404453 -lat_degrees = 38.954831 -meters_above_ellipsoid = 100.0 - -image_location = sensor_model.world_to_image( - GeodeticWorldCoordinate([radians(lon_degrees), - radians(lat_degrees), - meters_above_ellipsoid])) -``` - ## Contributing This project welcomes contributions and suggestions. If you would like to submit a pull request, see our diff --git a/doc/conf.py b/doc/conf.py index 2377b78..5544ef4 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -46,7 +46,7 @@ def setup(app): # -- Project information ----------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information -project = "OversightML Imagery Core" +project = "OversightML Imagery Toolkit" copyright = "{}, Amazon.com".format(datetime.datetime.now().year) author = "Amazon Web Services" diff --git a/doc/images/Photogrammetry-OODiagram.png b/doc/images/Photogrammetry-OODiagram.png new file mode 100644 index 0000000..27ed237 Binary files /dev/null and b/doc/images/Photogrammetry-OODiagram.png differ diff --git a/doc/images/SAR-HistogramStretchImage.png b/doc/images/SAR-HistogramStretchImage.png new file mode 100644 index 0000000..303f86c Binary files /dev/null and b/doc/images/SAR-HistogramStretchImage.png differ diff --git a/doc/images/SAR-QuarterPowerImage.png b/doc/images/SAR-QuarterPowerImage.png new file mode 100644 index 0000000..0874e8c Binary files /dev/null and b/doc/images/SAR-QuarterPowerImage.png differ diff --git a/doc/index.rst b/doc/index.rst index b6c5157..371ea9c 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -1,11 +1,11 @@ -aws-osml-imagery-core +osml-imagery-toolkit ===================== -The OversightML Imagery Core is a Python library that contains image processing and photogrammetry routines commonly +The OversightML Imagery Toolkit is a Python package that contains image processing and photogrammetry routines commonly used during the analysis of imagery collected by satellites and unmanned aerial vehicles (UAVs). It builds upon GDAL -by providing additional support for images compliant with the Sensor Independent Complex Data (SICD) and National -Imagery Transmission Format (NITF) standards. +by providing additional support for images compliant with the National Imagery Transmission Format (NITF), Sensor +Independent Complex Data (SICD), and Sensor Independent Derived Data (SIDD) standards. .. toctree:: :maxdepth: 1 diff --git a/src/aws/osml/features/__init__.py b/src/aws/osml/features/__init__.py index 5167245..7e04349 100644 --- a/src/aws/osml/features/__init__.py +++ b/src/aws/osml/features/__init__.py @@ -1,6 +1,10 @@ from .feature_index import Feature2DSpatialIndex, STRFeature2DSpatialIndex from .imaged_feature_property_accessor import ImagedFeaturePropertyAccessor +""" +The features package contains classes that assist with working with geospatial features derived from imagery. +""" + __all__ = [ "Feature2DSpatialIndex", "ImagedFeaturePropertyAccessor", diff --git a/src/aws/osml/gdal/__init__.py b/src/aws/osml/gdal/__init__.py index a5c7e08..19b5b64 100644 --- a/src/aws/osml/gdal/__init__.py +++ b/src/aws/osml/gdal/__init__.py @@ -3,6 +3,60 @@ # flake8: noqa """ The gdal package contains utilities that assist with loading imagery and metadata using the OSGeo GDAL library. + +Loading Imagery and Sensor Models with OSML +******************************************* + +OSML provides utilities to load a dataset and automatically construct an appropriate sensor model from metadata +available in the image. Metadata handled by GDAL (e.g. GeoTIFF tags or NITF segment metadata and TREs) is available +through the dataset accessors. + +.. code-block:: python + :caption: Example of loading a dataset and sensor model using OSML + + from aws.osml.gdal import load_gdal_dataset + + # Load the image and create a sensor model + dataset, sensor_model = load_gdal_dataset("./imagery/sample.nitf") + width = dataset.RasterXSize + height = dataset.RasterYSize + + print(f"Loaded image with dimensions: ({height}, {width}) (rows, cols)") + print(f"Using Sensor Model Implementation: {type(sensor_model).__name__}") + print(dataset.GetMetadata()) + + +Access to NITF Data Extension Segments +************************************** + +SICD and SIDD imagery contains additional metadata in a XML Data Extension Segment that is not currently processed +by GDAL. This information can be accessed with the help of the NITFDESAccessor. + +.. code-block:: python + :caption: Example of loading a dataset and sensor model using OSML + + import base64 + import xml.dom.minidom + from aws.osml.gdal import load_gdal_dataset, NITFDESAccessor + + dataset, sensor_model = load_gdal_dataset("./sample-sicd.nitf") + + des_accessor = NITFDESAccessor(dataset.GetMetadata("xml:DES")) + xml_data_content_segments = des_accessor.get_segments_by_name("XML_DATA_CONTENT") + if xml_data_content_segments is not None: + for xml_data_segment in xml_data_content_segments: + xml_bytes = des_accessor.parse_field_value(xml_data_segment, "DESDATA", base64.b64decode) + xml_str = xml_bytes.decode("utf-8") + if "SICD" in xml_str: + temp = xml.dom.minidom.parseString(xml_str) + new_xml = temp.toprettyxml() + print(new_xml) + break + +------------------------- + +APIs +**** """ from .gdal_config import GDALConfigEnv, set_gdal_default_configuration diff --git a/src/aws/osml/image_processing/__init__.py b/src/aws/osml/image_processing/__init__.py index fbb8d08..fc4be9d 100644 --- a/src/aws/osml/image_processing/__init__.py +++ b/src/aws/osml/image_processing/__init__.py @@ -3,9 +3,94 @@ # flake8: noqa """ The image_processing package contains various utilities for manipulating overhead imagery. + +Image Tiling: Tiling With Updated Image Metadata +************************************************ + +Many applications break large remote sensing images into smaller chips or tiles for distributed processing or +dissemination. GDAL's Translate function provides basic capabilities, but it does not correctly update geospatial +metadata to reflect the new image extent. These utilities provide those functions so tile consumers can correctly +interpret the pixel information they have been provided. + +.. code-block:: python + :caption: Example showing creation of a NITF tile from the upper left corner of an image + + # Load the image and create a sensor model + ds, sensor_model = load_gdal_dataset("./imagery/sample.nitf") + tile_factory = GDALTileFactory(ds, + sensor_model, + GDALImageFormats.NITF, + GDALCompressionOptions.NONE + ) + + # Bounds are [left_x, top_y, width, height] + nitf_encoded_tile_bytes = tile_factory.create_encoded_tile([0, 0, 1024, 1024]) + + +Image Tiling: Tiles for Display +******************************* + +Some images, for example 11-bit panchromatic images or SAR imagery with floating point complex data, can not be +displayed directly without remapping the pixels into an 8-bit per pixel grayscale or RGB color model. The TileFactory +supports creation of tiles suitable for human review by setting both the output_type and range_adjustment options. +Note that the output_size parameter can be used to generate lower resolution tiles. This operation will make use of +GDAL generated overviews if they are available to the dataset. + +.. code-block:: python + :caption: Example showing creation of a PNG tile scaled down from the full resolution image + + viz_tile_factory = GDALTileFactory(ds, + sensor_model, + GDALImageFormats.PNG, + GDALCompressionOptions.NONE, + output_type=gdalconst.GDT_Byte, + range_adjustment=RangeAdjustmentType.DRA) + + viz_tile = viz_tile_factory.create_encoded_tile([0, 0, 1024, 1024], output_size=(512, 512)) + + +Complex SAR Data Display +************************ + +There are a variety of different techniques to convert complex SAR data to a simple image suitable for human display. +The toolkit contains two helper functions that can convert complex image data into an 8-bit grayscle representation +The equations implemented are described in Sections 3.1 and 3.2 of SAR Image Scaling, Dynamic Range, Radiometric +Calibration, and Display (SAND2019-2371). + +.. code-block:: python + :caption: Example converting complex SAR data into a 8-bit per pixel image for display + + import numpy as np + from aws.osml.image_processing import histogram_stretch, quarter_power_image + + sicd_dataset, sensor_model = load_gdal_dataset("./sample-sicd.nitf") + complex_pixels = sicd_dataset.ReadAsArray() + + histo_stretch_pixels = histogram_stretch(complex_pixels) + quarter_power_pixels = quarter_power_image(complex_pixels) + + +.. figure:: ../images/SAR-HistogramStretchImage.png + :width: 400 + :alt: Histogram Stretch Applied to Sample SICD Image + + Example of applying histogram_stretch to a sample SICD image. + + +.. figure:: ../images/SAR-QuarterPowerImage.png + :width: 400 + :alt: Quarter Power Image Applied to Sample SICD Image + + Example of applying quarter_power_image to a sample SICD image. + + +------------------------- + +APIs +**** """ from .gdal_tile_factory import GDALTileFactory from .sar_complex_imageop import histogram_stretch, quarter_power_image -__all__ = ["GDALTileFactory"] +__all__ = ["GDALTileFactory", "histogram_stretch", "quarter_power_image"] diff --git a/src/aws/osml/photogrammetry/__init__.py b/src/aws/osml/photogrammetry/__init__.py index 247347b..b3191bd 100644 --- a/src/aws/osml/photogrammetry/__init__.py +++ b/src/aws/osml/photogrammetry/__init__.py @@ -2,8 +2,105 @@ # __init__.py file. # flake8: noqa """ -The photogrammetry package contains implementations of various sensor and elevation models used to convert between -the image (x, y) and geodetic (lon, lat, elev) coordinate systems. +Many users need to estimate the geographic position of an object found in a georeferenced image. The osml-imagery-toolkit +provides open source implementations of the image-to-world and world-to-image equations for some common replacement +sensor models. These sensor models work with many georeferenced imagery types and do not require orthorectification of +the image. In the current implementation support is provided for: + +* **Rational Polynomials**: Models based on rational polynomials specified using RSM and RPC metadata found in NITF TREs +* **SAR Sensor Independent Models**: Models as defined by the SICD and SIDD standards with metadata found in the NITF XML data segment. +* **Perspective and Affine Projections**: Simple matrix based projections that can be computed from geolocations of the 4 image corners or `tags found in GeoTIFF images `_. + +*Note that the current implementation does not support the RSM Grid based sensor models or the adjustable parameter +options. These features will be added in a future release.* + +.. figure:: ../images/Photogrammetry-OODiagram.png + :width: 400 + :alt: Photogrammetry Class Diagram + + Class diagram of the aws.osml.photogrammetry package showing public and private classes. + +Geolocating Image Pixels: Basic Examples +**************************************** + +Applications do not typically interact with a specific sensor model implementation directly. Instead, they let the +SensorModel abstraction encapsulate the details and rely on the image IO utilities to construct the appropriate +type given the available metadata. + +.. code-block:: python + :caption: Example showing calculation of an image location for a geodetic location + + dataset, sensor_model = load_gdal_dataset("./imagery/sample.nitf") + + lon_degrees = -77.404453 + lat_degrees = 38.954831 + meters_above_ellipsoid = 100.0 + + # Note the GeodeticWorldCoordinate is (longitude, latitude, elevation) with longitude and latitude in **radians** + # and elevation in meters above the WGS84 ellipsoid. The resulting ImageCoordinate is returned in (x, y) pixels. + image_location = sensor_model.world_to_image( + GeodeticWorldCoordinate([radians(lon_degrees), + radians(lat_degrees), + meters_above_ellipsoid])) + +.. code-block:: python + :caption: Example showing use of a SensorModel to geolocate 4 image corners + + dataset, sensor_model = load_gdal_dataset("./imagery/sample.nitf") + width = dataset.RasterXSize + height = dataset.RasterYSize + + image_corners = [[0, 0], [width, 0], [width, height], [0, height]] + geo_image_corners = [sensor_model.image_to_world(ImageCoordinate(corner)) + for corner in image_corners] + + # GeodeticWorldCoordinates have custom formatting defined that supports a variety of common output formats. + # The example shown below will produce a ddmmssXdddmmssY formatted coordinate (e.g. 295737N0314003E) + for geodetic_corner in geo_image_corners: + print(f"{geodetic_corner:%ld%lm%ls%lH%od%om%os%oH}") + +Geolocating Image Pixels: Addition of an External Elevation Model +***************************************************************** + +The image-to-world calculation can optionally use an external digital elevation model (DEM) when geolocating points +on an image. How the elevation model will be used varies by sensor model but examples include: + +* Using DEM elevations as a constraint during iterations of a rational polynomial camera's image-to-world calculation. +* Computing the intersection of a R/Rdot contour with a DEM for sensor independent SAR models. + +All of these approaches make the fundamental assumption that the pixel lies on the terrain surface. If a DEM is not +available we assume a constant surface with elevation provided in the image metadata. + +.. code-block:: python + :caption: Example showing use of an external SRTM DEM to provide elevation data for image center + + ds, sensor_model = load_gdal_dataset("./imagery/sample.nitf") + elevation_model = DigitalElevationModel( + SRTMTileSet(version="1arc_v3"), + GDALDigitalElevationModelTileFactory("./local-SRTM-tiles")) + + # Note the order of ImageCoordinate is (x, y) and the resulting geodetic coordinate is + # (longitude, latitude, elevation) with longitude and latitude in **radians** and elevation in meters. + geodetic_location_of_image_center = sensor_model.image_to_world( + ImageCoordinate([ds.RasterXSize/2, ds.RasterYSize/2]), + elevation_model=elevation_model) + + +External References +******************* + +* Manual of Photogrammetry: https://www.amazon.com/Manual-Photogrammetry-PhD-Chris-McGlone/dp/1570830991 +* NITF Compendium of Controlled Support Data Extensions: https://nsgreg.nga.mil/doc/view?i=5417 +* The Replacement Sensor Model (RSM): Overview, Status, and Performance Summary: https://citeseerx.ist.psu.edu/doc_view/pid/c25de8176fe790c28cf6e1aff98ecdea8c726c96 +* RPC Whitepaper: https://users.cecs.anu.edu.au/~hartley/Papers/cubic/cubic.pdf +* SICD Volume 3, Image Projections Description Document: https://nsgreg.nga.mil/doc/view?i=5383 +* WGS84 Standard: https://nsgreg.nga.mil/doc/view?i=4085 + +------------------------- + +APIs +**** + """ from .chipped_image_sensor_model import ChippedImageSensorModel @@ -48,27 +145,25 @@ __all__ = [ "ChippedImageSensorModel", "CompositeSensorModel", - "GeodeticWorldCoordinate", - "ImageCoordinate", - "WorldCoordinate", - "geocentric_to_geodetic", - "geodetic_to_geocentric", + "ConstantElevationModel", "DigitalElevationModel", "DigitalElevationModelTileFactory", "DigitalElevationModelTileSet", - "ConstantElevationModel", "ElevationModel", "ElevationRegionSummary", "GDALAffineSensorModel", "GenericDEMTileSet", - "SARImageCoordConverter", + "GeodeticWorldCoordinate", "INCAProjectionSet", - "PlaneProjectionSet", + "ImageCoordinate", "PFAProjectionSet", - "PolynomialXYZ", + "PlaneProjectionSet", "Polynomial2D", + "PolynomialXYZ", "ProjectiveSensorModel", "RGAZCOMPProjectionSet", + "RPCPolynomial", + "RPCSensorModel", "RSMContext", "RSMGroundDomain", "RSMGroundDomainForm", @@ -77,10 +172,12 @@ "RSMPolynomial", "RSMPolynomialSensorModel", "RSMSectionedPolynomialSensorModel", - "RPCPolynomial", - "RPCSensorModel", - "SensorModel", - "SensorModelOptions", + "SARImageCoordConverter", "SICDSensorModel", "SRTMTileSet", + "SensorModel", + "SensorModelOptions", + "WorldCoordinate", + "geocentric_to_geodetic", + "geodetic_to_geocentric", ] diff --git a/src/aws/osml/photogrammetry/coordinates.py b/src/aws/osml/photogrammetry/coordinates.py index 30fd187..d2718df 100644 --- a/src/aws/osml/photogrammetry/coordinates.py +++ b/src/aws/osml/photogrammetry/coordinates.py @@ -64,7 +64,7 @@ class GeodeticWorldCoordinate(WorldCoordinate): This class uses a custom format specification for a geodetic coordinate uses % directives similar to datetime. These custom directives can be combined as needed with literal values to produce a wide - range of output formats. For example, '%ld%lm%ls%lH%od%om%os%oh' will produce a ddmmssXdddmmssY formatted + range of output formats. For example, '%ld%lm%ls%lH%od%om%os%oH' will produce a ddmmssXdddmmssY formatted coordinate. The first half, ddmmssX, represents degrees, minutes, and seconds of latitude with X representing North or South (N for North, S for South). The second half, dddmmssY, represents degrees, minutes, and seconds of longitude with Y representing East or West (E for East, W for West), respectively.