Split Polygon along IDL to go into CMR

shape_split.py

import logging
import statistics
from math import sqrt

from shapely import wkt
from shapely.ops import linemerge, unary_union, polygonize
from shapely.geometry import Polygon

# Move the box out of earth coordinate system
def shift_polygon(bbox):
    shifted = []
    for p in bbox:
        shifted.append([(360.0 + float(p[1])) % 360, float(p[0])])
    return Polygon(shifted)

# Split a shape along IDL
def split_up_polygon(bbox):
    polygon = Polygon(bbox)
    idl_line = wkt.loads("LINESTRING(180 90, 180 -90)")
    merged = linemerge([polygon.boundary, idl_line])
    borders = unary_union(merged)
    return polygonize(borders)

# Check for IDL Crossing polygons
def does_polygon_cross_idl(bbox):
    lon = [float(p[1]) for p in bbox]
    max_lon, min_lon = max(lon), min(lon)
    if min_lon >= -160 or max_lon <= 160:
        logging.info("Does not touch IDL")
        return False
    logging.info("Granule intersects IDL")
    return True


# Turn a polgyon into 2 polygons
# INPUT: is a list of polygon points, in [LAT, LON] format WITH closure point!
#    bbox = [ [ LAT1, LON1 ], [ LAT2, LON2 ], [ LAT3, LON3 ], [ LAT4, LON5 ], [ LAT1, LON1 ] ]
# OUTPUT:  A list of polyons in [LAT, LON] WITHOUT closer points!
#  return ( [ [ [ LAT1, LON1 ], [ LAT2, LON2 ], ... ], [ LAT1, LON1 ], [ LAT2, LON2 ] ... ] ] )

def split_polygon_on_idl(bbox):
    # Shift away from IDL
    shifted_polygon = shift_polygon(bbox)
    polygons = split_up_polygon(shifted_polygon)

    # Shift back
    new_bboxes = []
    for p in polygons:
        new_bbox = []
        last_point = None

        # Calculate the center of the shape, this is used for fixing points ON The IDL.
        shape_avg = statistics.mean([p[0] for p in list(p.boundary.coords)]) - 180
        
        # Move each point back to proper earth coordinates
        for pt in list(p.boundary.coords):
            if last_point and (sqrt((last_point[0] - pt[0]) ** 2 + (last_point[1] - pt[1]) ** 2) < 0.01):
                # CMR will reject points the are too close. If splitting result in VERY short line segments, drop one
                logging.warning(f"points {last_point} and {pt} are too close. {pt} will be dropped")
                continue

            if pt[0] > 180.0:
                # Points > 180º are actaully on the other side of the globe
                new_bbox.append([pt[1], pt[0] - 360])
                
            elif pt[0] == 180.0 and shape_avg > 0:
                # if the point it = 180.0, and we're mostly in positive, use -180
                # This seems counter-intuitive, but it works
                new_bbox.append([pt[1], -179.999])
            elif pt[0] == 180.0 and shape_avg < 0:
                new_bbox.append([pt[1], 179.999])
          
            else:
                # this point is right where it should be
                new_bbox.append([pt[1], pt[0]])
                
            # Watch for points that are too close
            last_point = pt

        # CMR does not line small pieces. 
        if len(new_bbox) >= 4:
            del new_bbox[-1]  # CMR does not want closure point
            new_bboxes.append(new_bbox)
        else:
            # Shapes with less than 4 points are straight, arealess lines.
            # This happens when shapes only BARELY cross the IDL in one corner.
            # In this case, we just throw away the splinter
            logging.debug("IDL crossing splinter is being ignored")

    return new_bboxes