rs_server_common/footprint_facility.md

<< Back to index

Copyright 2024 - GAEL Systems

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.

check_cross_antimeridian(geometry)

Checks if the geometry pass over ±180 longitude position. The detection of antimeridian is performed according to the distance between longitudes positions between consecutive points of the geometry points. The distance shall be greater than 180 to avoid revert longitude signs around greenwich meridian 0. It is also considered crossing antimeridian when one of the polygon longitude is exactly to +-180 degrees. :parameter geometry: The geometry to be controlled. :return: True if the geometry pass over antimeridian, False otherwise.

Source code in docs/rs-server/services/common/rs_server_common/footprint_facility.py

def check_cross_antimeridian(geometry: Geometry) -> bool:
    """
    Checks if the geometry pass over ±180 longitude position.
    The detection of antimeridian is performed according to the distance
    between longitudes positions between consecutive points of the geometry
    points. The distance shall be greater than 180 to avoid revert longitude
    signs around greenwich meridian 0.
    It is also considered crossing antimeridian when one of the polygon
    longitude is exactly to +-180 degrees.
    :parameter geometry: The geometry to be controlled.
    :return: True if the geometry pass over antimeridian, False otherwise.
    """
    # Case of Collection of geometries (i.e. Multipolygons)
    if hasattr(geometry, "geoms"):
        for geom in geometry.geoms:
            if check_cross_antimeridian(geom):
                return True

    # Path of points shall exist (Polygon or Linestring)
    boundary = np.array(shapely.get_coordinates(geometry))
    i = 0
    while i < boundary.shape[0] - 1:
        if boundary[i, 0] == 180 or boundary[i, 0] == -180 or abs(boundary[i + 1, 0] - boundary[i, 0]) > 180:
            return True
        i += 1
    return False

rework_to_polygon_geometry(geometry)

Rework the geometry to manage polar and antimeridian singularity. This process implements the Polar inclusive algorithm. The objective of this algorithm is to add the North/South Pole into the list of coordinates of geometry polygon at the antimeridian cross.

When the geometry contains the pole the single polygon geometry including the pole in its border point list is properly interpreted by displays systems. When the geometry does not contain the pole, the geometry shall be split among the antimeridian line.

:param geometry: the geometry crossing the antimeridian. :return: the modified geometry with the closest pole included at antimeridian crossing.

Source code in docs/rs-server/services/common/rs_server_common/footprint_facility.py

def rework_to_polygon_geometry(geometry: Geometry):
    """Rework the geometry to manage polar and antimeridian singularity.
    This process implements the **Polar inclusive algorithm**.
    The objective of this algorithm is to add the North/South Pole into
    the list of coordinates of geometry polygon at the antimeridian cross.

    When the geometry contains the pole the single polygon geometry including
    the pole in its border point list is properly interpreted by displays
    systems. When the geometry does not contain the pole, the geometry shall be
    split among the antimeridian line.

    :param geometry: the geometry crossing the antimeridian.
    :return: the modified geometry with the closest pole included at
     antimeridian crossing.
    """
    if not check_cross_antimeridian(geometry):
        return geometry

    count, mixed, latitudes = _num_cross_antimeridian(geometry)
    if mixed and count > 2:
        logger.debug(f"WARN: Crossing antimeridian {count} times ...")
        logger.debug(f" And crossing equator {_num_cross_equator(geometry)} times ...")
        logger.debug("Split footprint at equator to support polar inclusion")
        geometry = _split_crude_polygon_to_equator(geometry)
        rwrk = []
        for geom in geometry.geoms:
            reworked = rework_to_polygon_geometry(geom)
            if isinstance(reworked, BaseMultipartGeometry):
                for geo in reworked.geoms:
                    rwrk.append(geo)
            else:
                rwrk.append(reworked)
        rwrk = _merge_polygon_to_equator(rwrk)
        return rwrk

    # Manage case of multipolygon input
    if isinstance(geometry, shapely.geometry.MultiPolygon):
        rwrk = []
        for geom in geometry.geoms:
            reworked = rework_to_polygon_geometry(geom)
            if isinstance(reworked, BaseMultipartGeometry):
                raise AlreadyReworkedPolygonError("Algorithm not supported already reworked inputs.")
            rwrk.append(reworked)
        return shapely.geometry.MultiPolygon(rwrk)

    if not isinstance(geometry, shapely.geometry.Polygon):
        raise ValueError("Polygon and MultiPolygon features only are " f"supported ({type(geometry).__name__})")

    boundaries = np.array(shapely.get_coordinates(geometry))
    i = 0
    vertical_set = False
    vsign = 1
    while i < boundaries.shape[0] - 1:
        if abs(boundaries[i + 1, 0] - boundaries[i, 0]) > 180:
            if not vertical_set:
                vsign = -1 if boundaries[i, 1] < 0 else 1
                vertical_set = True
            hsign = -1 if boundaries[i, 0] < 0 else 1
            lat = _lat_cross_antimeridian(boundaries[i], boundaries[i + 1])
            boundaries = np.insert(
                boundaries,
                i + 1,
                [
                    [hsign * 180, lat],
                    [hsign * 180, vsign * 90],
                    [-hsign * 180, vsign * 90],
                    [-hsign * 180, lat],
                ],
                axis=0,
            )
            i += 5
        else:
            i += 1
    geometry_type = type(geometry)
    reworked = geometry_type(boundaries)

    # When the geometry does not contain pole: Cuts the geometry among the
    # antimeridian line.
    if not _check_contains_pole(reworked):
        reworked = _split_polygon_to_antimeridian(reworked)
    else:
        # Case of footprint crossing equator, antimeridian and polar zone
        # Split at the equator
        # Warn:
        if _num_cross_equator(reworked) > 1:
            reworked = _split_polygon_to_equator(reworked)
        # When footprint contains overlapping, it happens at polar location.
        # Polygon containing overlapping are considered invalid in shapely
        # library. It includes validity check and correction methods.
        # The shapely correction method extrude the overlap areas and fails
        # to generate patchwork of polygons at polar area. This is probably
        # due to the antimeridian crossing.
        # Shapely "buffer" method fixe"s the geometry merging overlapping
        # regions.
        if not shapely.is_valid(reworked):
            reworked = (
                shapely.make_valid(reworked, method="structure", keep_collapsed=False)
                if _SHAPELY_SUPPORTS_NEW_MAKE_VALID
                else shapely.make_valid(reworked)
            )

    return shapely.buffer(reworked, 0)