#!/usr/bin/env python
# coding: utf8
#
# Copyright (c) 2023 Centre National d'Etudes Spatiales (CNES).
#
# This file is part of shareloc
# (see https://github.com/CNES/shareloc).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
"""
This module contains the GeoModel abstract class
"""
# Standard imports
from abc import abstractmethod
from typing import Union
import numpy as np
from affine import Affine
import bindings_cpp
from shareloc.geofunctions.dtm_intersection import DTMIntersection
from shareloc.proj_utils import coordinates_conversion, transform_index_to_physical_point
[docs]
class GeoModelTemplate:
"""
Class for general specification of a geometric model
declined in rpc.py and grid.py and rpc_optim.py
"""
@abstractmethod
def __init__(self):
"""
Return the geomodel object associated with the geomodel_type
given in the configuration
"""
# geomodel type. Set by the subclass
# geomodel epsg projection code
# Define GeoModelTemplate functions interface
@abstractmethod
[docs]
def direct_loc_h(self, row, col, alt, fill_nan=False):
"""
direct localization at constant altitude
:param row: line sensor position
:type row: float or 1D numpy.ndarray dtype=float64
:param col: column sensor position
:type col: float or 1D numpy.ndarray dtype=float64
:param alt: altitude
:param fill_nan: fill numpy.nan values with lon and lat offset if true (same as OTB/OSSIM), nan is returned
otherwise
:type fill_nan: boolean
:return: ground position (lon,lat,h)
:rtype: numpy.ndarray 2D dimension with (N,3) shape, where N is number of input coordinates
"""
@abstractmethod
[docs]
def direct_loc_dtm(self, row, col, dtm):
"""
direct localization on dtm
:param row: line sensor position
:type row: float
:param col: column sensor position
:type col: float
:param dtm: dtm intersection model
:type dtm: shareloc.geofunctions.dtm_intersection
:return: ground position (lon,lat,h) in dtm coordinates system
:rtype: numpy.ndarray 2D dimension with (N,3) shape, where N is number of input coordinates
"""
@abstractmethod
[docs]
def inverse_loc(self, lon, lat, alt):
"""
Inverse localization
:param lon: longitude position
:type lon: float or 1D numpy.ndarray dtype=float64
:param lat: latitude position
:type lat: float or 1D numpy.ndarray dtype=float64
:param alt: altitude
:type alt: float
:return: sensor position (row, col, alt)
:rtype: tuple(1D np.array row position, 1D np.array col position, 1D np.array alt)
"""
@classmethod
@abstractmethod
[docs]
def load(cls, geomodel_path: str):
"""
load function with class specific args
:param geomodel_path: filename of geomodel
"""
[docs]
def check_lonlat(self, lonlath: np.ndarray) -> np.ndarray:
"""
returns filtered latlonh. abs(longitude) should be <= 180, abs(latitude) should be <=90,
and no NaN in coordinates.
:param latlonh: np.array of latlonh
:return: filtered array
"""
filtering = np.logical_or(
np.logical_or(np.abs(lonlath[:, 0]) > 180.0, np.abs(lonlath[:, 1]) > 90.0),
np.sum(np.isnan(lonlath), axis=1) > 0,
)
lonlath[filtering, :] = np.nan
return lonlath
[docs]
def get_dtm_alt_offset(self, corners: np.ndarray, dtm: Union[DTMIntersection, bindings_cpp.DTMIntersection]):
"""
returns min/max altitude offset between dtm coordinates system and Geomodel one
:param corners: corners of the DTM's footprint
:type corners: np.ndarray (4x2)
:param dtm: DTM to get alt offset from
:type dtm: DTMIntersection or bindings_cpp.DTMIntersection
:return: min/max altimetric difference between RPC's epsg minus dtm alti expressed in dtm epsg
:rtype: list of float (1x2)
"""
alti_moy = (dtm.get_alt_min() + dtm.get_alt_max()) / 2.0
ground_corners = np.zeros([3, 4])
ground_corners[2, :] = alti_moy
transform = dtm.get_transform()
if not isinstance(transform, Affine):
transform = Affine.from_gdal(*transform)
ground_corners[1::-1, :] = transform_index_to_physical_point(transform, corners[:, 0], corners[:, 1])
converted_corners = coordinates_conversion(ground_corners.transpose(), dtm.get_epsg(), self.epsg)
return [np.min(converted_corners[:, 2]) - alti_moy, np.max(converted_corners[:, 2]) - alti_moy]