Add support for RPCs

[enomis-dev]

Currently Telluric doesn't support RPCs (https://en.wikipedia.org/wiki/Rational_polynomial_coefficient). This issue aims to centralize discussions related to this topic. I think that the following features should be added: 1) support for rpcs in constructors (in this context it should be decided which class should be used to represent RPCs) 2) property for accessing them 3) existing interfaces such as footprint(), corners() should continue to work for rpc-based rasters as well 4) reprojection should support rpc based reprojection 5) relevant methods (such as crop(), resize()) should update rpcs as it is done for affines

[enomis-dev]

My first pull request covers points 1,2 and 4. Points 3 and 5 still need to be addressed and require more discussion. For these last points I link a very interesting project that could be used in Telluric to support operations (such as footprint, corners(), crop() and resize()) with rasters that have only rpcs https://github.com/centreborelli/rpcm

[enomis-dev]

About my previous post I would like to add https://github.com/centreborelli/rpcm as dependency to Telluric to deal with rpcs. It allows to convert image coordinates plus altitude into geographic coordinates. This could be exploited in Telluric to retrieve raster footprint or corners without having to rewrite new code. What do you think? @AmitAronovitch @drnextgis

[drnextgis]

Interesting.

1. Could you please give an example of which sort of problems cannot be solved without that dependency?
2. Does it make sense to consider extracting needed piece of code from rpcm (if the license allows that) and enrich RPC class.

[enomis-dev]

I was thinking to use the rpcm module to implement footprint and corner methods for raster provided with rpc but without crs and with default affine. The project is under the license BSD 2-Clause License. Example:

from telluric import GeoRaster2
from telluric import GeoVector
from rpcm.rpc_model import rpc_from_geotiff
from rpcm import RPCModel
import numpy as np
from telluric.constants import WGS84_CRS
from telluric.util.raster_utils import calc_transform
from types import SimpleNamespace

# georaster without crs and with default affine
georaster = GeoRaster2.open("grayscale.tif")

# Both method will raise CRSError: CRS is invalid: None
# georaster.footprint()
# georaster.corners()

Possible implementation with RPCM modul

def convert_rasterio_rpcs_to_RPCM(rpcs):
    """
    Convert RPCs from rasterio.rpc.RPC to a RPCModel

    Args:
        rpcs (rasterio.rpc.RPC): raster rpcs

    Returns:
        rpcs (rpcm.rpc_model.RPCModel): raster rpcs
    """
    rpcs_dict = rpcs.to_dict()
    rpcs_dict = {k.upper(): v for k, v in rpcs_dict.items()}
    line_num_coeff = rpcs_dict['LINE_NUM_COEFF']
    rpcs_dict['LINE_NUM_COEFF'] = ' '.join(str(e) for e in line_num_coeff)
    line_den_coeff = rpcs_dict['LINE_DEN_COEFF']
    rpcs_dict['LINE_DEN_COEFF'] = ' '.join(str(e) for e in line_den_coeff)
    samp_num_coeff = rpcs_dict['SAMP_NUM_COEFF']
    rpcs_dict['SAMP_NUM_COEFF'] = ' '.join(str(e) for e in samp_num_coeff)
    samp_den_coeff = rpcs_dict['SAMP_DEN_COEFF']
    rpcs_dict['SAMP_DEN_COEFF'] = ' '.join(str(e) for e in samp_den_coeff)
    return RPCModel(rpcs_dict)

def rpc_raster_footprint(raster, rpc=None):
    """
    Retrieve raster footprint from rpcs

    Args:
        raster (GeoRaster2): raster
        rpcs (rasterio.rpc.RPC): raster rpcs

    Returns:
        GeoVector: raster footprint
    """
    if rpc is None:
        rpc = raster.rpcs
    rpc = convert_rasterio_rpcs_to_RPCM(rpc)
    corners = np.vstack([np.array(raster.image_corner(corner)) for corner in raster.corner_types()])
    # rpcm localization function
    coords = rpc.localization(corners[:, 0], corners[:, 1], rpc.alt_offset)
    return GeoVector.polygon(np.vstack(coords).T, crs=WGS84_CRS)

def rpc_raster_corner(raster, corner, rpc=None):
    """
    Retrieve raster corners from rpcs

    Args:
        raster (GeoRaster2): raster
        rpcs (rasterio.rpc.RPC): raster rpcs

    Returns:
        GeoVector: raster corners
    """
    if rpc is None:
        rpc = raster.rpcs
    rpc = convert_rasterio_rpcs_to_RPCM(rpc)
    img_corner = raster.image_corner(corner)
    # rpcm localization function
    lon, lat = rpc.localization(img_corner.x, img_corner.y, rpc.alt_offset)
    return GeoVector.point(lon, lat, crs=WGS84_CRS)

georaster_bounds = rpc_raster_footprint(georaster)
georaster_bounds

GeoVector(shape=POLYGON ((111.446359536235 34.90388330200413, 111.391204658326 34.91225086092173, 111.3934946814279 34.92250972381049, 111.4486565101207 34.91414140921886, 111.446359536235 34.90388330200413)), crs=EPSG:4326)

rpc_raster_corner(georaster, 'ul')

GeoVector(shape=POINT (111.446359536235 34.90388330200413), crs=EPSG:4326)

However actually today I thought to a possible implementation without RPCM model class that however has some different results:

def rpc_raster_footprint_telluric(raster, dst_crs):
    src = SimpleNamespace(width=raster.width, height=raster.height, transform=raster.transform, crs=raster.crs)
    dst_crs, dst_transform, _, _ = calc_transform(src, dst_crs=dst_crs, rpcs=raster.rpcs)
    new_raster = raster.copy_with(crs=dst_crs, affine=dst_transform)
    footprint = new_raster.footprint()
    return footprint

def rpc_raster_corner_telluric(raster, corner, dst_crs):
    src = SimpleNamespace(width=raster.width, height=raster.height, transform=raster.transform, crs=raster.crs)
    dst_crs, dst_transform, _, _ = calc_transform(src, dst_crs=dst_crs, rpcs=raster.rpcs)
    new_raster = raster.copy_with(crs=dst_crs, affine=dst_transform)
    return new_raster.corner(corner)

rpc_raster_footprint_telluric(georaster, WGS84_CRS)

GeoVector(shape=POLYGON ((111.3911674407107 34.92250008283582, 111.4459688804706 34.92250008283582, 111.4459688804706 34.91012694526505, 111.3911674407107 34.91012694526505, 111.3911674407107 34.92250008283582)), crs=EPSG:4326)

rpc_raster_corner_telluric(georaster, 'ul', WGS84_CRS)

GeoVector(shape=POINT (111.3911674407107 34.92250008283582), crs=EPSG:4326)

However, as you see we have different results. Wdyt?

[drnextgis]

Can you explain the difference in results? Which one is more accurate?

[enomis-dev]

The differences in results come from the fact that the second solution is based on Telluric and so then rasterio. The first solution is based on this rpcm model that uses an iterative algorithm. I compared the results with QGIS. This is the footprint of the same raster from QGIS:

[[[[111.39116755378235, 34.922500283177854],
   [111.44875281264665, 34.922501113211695],
   [111.44875239762972, 34.90383157669373],
   [111.39116755378238, 34.90383282174453],
   [111.39116755378235, 34.922500283177854]]]])])

As it can be easily seen the footprint is almost equal to the one calculated with "Telluric" solution. I propose to not add other dependencies to Telluric and implement the methods I suggested that rely only on Telluric functions.