rslearn.utils.geometry

geometry

Spatiotemporal geometry utilities.

Projection

A projection specifies a CRS, x resolution, and y resolution.

The coordinate reference system (CRS) defines the meaning of the coordinates. The resolutions specify the pixels per projection unit, and are used to map pixel coordinates to CRS coordinates.

Source code in rslearn/utils/geometry.py

class Projection:
    """A projection specifies a CRS, x resolution, and y resolution.

    The coordinate reference system (CRS) defines the meaning of the coordinates. The
    resolutions specify the pixels per projection unit, and are used to map pixel
    coordinates to CRS coordinates.
    """

    def __init__(self, crs: CRS, x_resolution: float, y_resolution: float) -> None:
        """Initialize a new Projection.

        Args:
            crs: the CRS
            x_resolution: the x resolution
            y_resolution: the y resolution
        """
        self.crs = crs
        self.x_resolution = x_resolution
        self.y_resolution = y_resolution

    def __eq__(self, other: Any) -> bool:
        """Returns whether this projection is the same as the other projection."""
        if not isinstance(other, Projection):
            return False
        if self.crs != other.crs:
            return False
        if not is_same_resolution(self.x_resolution, other.x_resolution):
            return False
        if not is_same_resolution(self.y_resolution, other.y_resolution):
            return False
        return True

    def __repr__(self) -> str:
        """Returns a string representation of this projection."""
        return (
            f"Projection(crs={self.crs}, "
            + f"x_resolution={self.x_resolution}, "
            + f"y_resolution={self.y_resolution})"
        )

    def __str__(self) -> str:
        """Returns a human-readable string summary of this projection."""
        return f"{self.crs}_{self.x_resolution}_{self.y_resolution}"

    def __hash__(self) -> int:
        """Returns a hash of this projection."""
        return hash((self.crs, self.x_resolution, self.y_resolution))

    def serialize(self) -> dict:
        """Serializes the projection to a JSON-encodable dictionary."""
        return {
            "crs": self.crs.to_string(),
            "x_resolution": self.x_resolution,
            "y_resolution": self.y_resolution,
        }

    @staticmethod
    def deserialize(d: dict) -> "Projection":
        """Deserializes a projection from a JSON-decoded dictionary."""
        return Projection(
            crs=CRS.from_string(d["crs"]),
            x_resolution=d["x_resolution"],
            y_resolution=d["y_resolution"],
        )

serialize

serialize() -> dict

Serializes the projection to a JSON-encodable dictionary.

Source code in rslearn/utils/geometry.py

def serialize(self) -> dict:
    """Serializes the projection to a JSON-encodable dictionary."""
    return {
        "crs": self.crs.to_string(),
        "x_resolution": self.x_resolution,
        "y_resolution": self.y_resolution,
    }

deserialize staticmethod

deserialize(d: dict) -> Projection

Deserializes a projection from a JSON-decoded dictionary.

Source code in rslearn/utils/geometry.py

@staticmethod
def deserialize(d: dict) -> "Projection":
    """Deserializes a projection from a JSON-decoded dictionary."""
    return Projection(
        crs=CRS.from_string(d["crs"]),
        x_resolution=d["x_resolution"],
        y_resolution=d["y_resolution"],
    )

ResolutionFactor

Multiplier for the resolution in a Projection.

The multiplier is either an integer x, or the inverse of an integer (1/x).

Factors greater than 1 increase the projection_units/pixel resolution, increasing the resolution (more pixels per projection unit). Factors less than 1 make it coarser (less pixels).

Source code in rslearn/utils/geometry.py

class ResolutionFactor:
    """Multiplier for the resolution in a Projection.

    The multiplier is either an integer x, or the inverse of an integer (1/x).

    Factors greater than 1 increase the projection_units/pixel resolution, increasing
    the resolution (more pixels per projection unit). Factors less than 1 make it coarser
    (less pixels).
    """

    def __init__(self, numerator: int = 1, denominator: int = 1):
        """Create a new ResolutionFactor.

        Args:
            numerator: the numerator of the fraction.
            denominator: the denominator of the fraction. If set, numerator must be 1.
        """
        if numerator != 1 and denominator != 1:
            raise ValueError("one of numerator or denominator must be 1")
        if not isinstance(numerator, int) or not isinstance(denominator, int):
            raise ValueError("numerator and denominator must be integers")
        if numerator < 1 or denominator < 1:
            raise ValueError("numerator and denominator must be >= 1")
        self.numerator = numerator
        self.denominator = denominator

    def multiply_projection(self, projection: Projection) -> Projection:
        """Multiply the projection by this factor."""
        if self.denominator > 1:
            return Projection(
                projection.crs,
                projection.x_resolution * self.denominator,
                projection.y_resolution * self.denominator,
            )
        else:
            return Projection(
                projection.crs,
                projection.x_resolution / self.numerator,
                projection.y_resolution / self.numerator,
            )

    def multiply_bounds(self, bounds: PixelBounds) -> PixelBounds:
        """Multiply the bounds by this factor.

        When coarsening, the width and height of the given bounds must be a multiple of
        the denominator.
        """
        if self.denominator > 1:
            # Verify the width and height are multiples of the denominator.
            # Otherwise the new width and height is not an integer.
            width = bounds[2] - bounds[0]
            height = bounds[3] - bounds[1]
            if width % self.denominator != 0 or height % self.denominator != 0:
                raise ValueError(
                    f"width {width} or height {height} is not a multiple of the resolution factor {self.denominator}"
                )
            # TODO: an offset could be introduced by bounds not being a multiple
            # of the denominator -> will need to decide how to handle that.
            return (
                bounds[0] // self.denominator,
                bounds[1] // self.denominator,
                bounds[2] // self.denominator,
                bounds[3] // self.denominator,
            )
        else:
            return (
                bounds[0] * self.numerator,
                bounds[1] * self.numerator,
                bounds[2] * self.numerator,
                bounds[3] * self.numerator,
            )

multiply_projection

multiply_projection(projection: Projection) -> Projection

Multiply the projection by this factor.

Source code in rslearn/utils/geometry.py

def multiply_projection(self, projection: Projection) -> Projection:
    """Multiply the projection by this factor."""
    if self.denominator > 1:
        return Projection(
            projection.crs,
            projection.x_resolution * self.denominator,
            projection.y_resolution * self.denominator,
        )
    else:
        return Projection(
            projection.crs,
            projection.x_resolution / self.numerator,
            projection.y_resolution / self.numerator,
        )

multiply_bounds

multiply_bounds(bounds: PixelBounds) -> PixelBounds

Multiply the bounds by this factor.

When coarsening, the width and height of the given bounds must be a multiple of the denominator.

Source code in rslearn/utils/geometry.py

def multiply_bounds(self, bounds: PixelBounds) -> PixelBounds:
    """Multiply the bounds by this factor.

    When coarsening, the width and height of the given bounds must be a multiple of
    the denominator.
    """
    if self.denominator > 1:
        # Verify the width and height are multiples of the denominator.
        # Otherwise the new width and height is not an integer.
        width = bounds[2] - bounds[0]
        height = bounds[3] - bounds[1]
        if width % self.denominator != 0 or height % self.denominator != 0:
            raise ValueError(
                f"width {width} or height {height} is not a multiple of the resolution factor {self.denominator}"
            )
        # TODO: an offset could be introduced by bounds not being a multiple
        # of the denominator -> will need to decide how to handle that.
        return (
            bounds[0] // self.denominator,
            bounds[1] // self.denominator,
            bounds[2] // self.denominator,
            bounds[3] // self.denominator,
        )
    else:
        return (
            bounds[0] * self.numerator,
            bounds[1] * self.numerator,
            bounds[2] * self.numerator,
            bounds[3] * self.numerator,
        )

STGeometry

A spatiotemporal geometry.

Specifiec crs and resolution and corresponding shape in pixel coordinates. Also specifies an optional time range (time range is unlimited if unset).

Source code in rslearn/utils/geometry.py

class STGeometry:
    """A spatiotemporal geometry.

    Specifiec crs and resolution and corresponding shape in pixel coordinates. Also
    specifies an optional time range (time range is unlimited if unset).
    """

    def __init__(
        self,
        projection: Projection,
        shp: shapely.Geometry,
        time_range: tuple[datetime, datetime] | None,
    ):
        """Creates a new spatiotemporal geometry.

        Args:
            projection: the projection
            shp: the shape in pixel coordinates
            time_range: optional start and end time (default unlimited)
        """
        self.projection = projection
        self.shp = shp
        self.time_range = time_range

    def contains_time(self, time: datetime) -> bool:
        """Returns whether this box contains the time."""
        if self.time_range is None:
            return True
        return time >= self.time_range[0] and time < self.time_range[1]

    def distance_to_time(self, time: datetime) -> timedelta:
        """Returns the distance from this box to the specified time.

        Args:
            time: the time to compute distance from

        Returns:
            the distance, which is 0 if the box contains the time
        """
        if self.time_range is None:
            return timedelta()
        if time < self.time_range[0]:
            return self.time_range[0] - time
        if time > self.time_range[1]:
            return time - self.time_range[1]
        return timedelta()

    def distance_to_time_range(
        self, time_range: tuple[datetime, datetime] | None
    ) -> timedelta:
        """Returns the distance from this geometry to the specified time range.

        Args:
            time_range: the time range to compute distance from

        Returns:
            the distance, which is 0 if the time ranges intersect
        """
        if self.time_range is None or time_range is None:
            return timedelta()
        if time_range[1] < self.time_range[0]:
            return self.time_range[0] - time_range[1]
        if self.time_range[1] < time_range[0]:
            return time_range[0] - self.time_range[1]
        return timedelta()

    def intersects_time_range(
        self, time_range: tuple[datetime, datetime] | None
    ) -> bool:
        """Returns whether this geometry intersects the other time range."""
        if self.time_range is None or time_range is None:
            return True
        if self.time_range[1] <= time_range[0]:
            return False
        if time_range[1] <= self.time_range[0]:
            return False
        return True

    def is_global(self) -> bool:
        """Returns whether this geometry has global spatial coverage.

        Global coverage is indicated by a special geometry with WGS84 projection and
        corners at (-180, -90, 180, 90) (see get_global_geometry).
        """
        if self.projection != WGS84_PROJECTION:
            return False
        if self.shp != shapely.box(*WGS84_BOUNDS):
            return False
        return True

    def is_too_large(self) -> bool:
        """Returns whether this geometry's spatial coverage is too large.

        This means that it will likely have issues during re-projections and such.
        """
        wgs84_bounds = self.to_projection(WGS84_PROJECTION).shp.bounds
        if wgs84_bounds[2] - wgs84_bounds[0] > MAX_GEOMETRY_DEGREES:
            return True
        if wgs84_bounds[3] - wgs84_bounds[1] > MAX_GEOMETRY_DEGREES:
            return True
        return False

    def to_wgs84(self) -> "STGeometry":
        """Convert to WGS84 with antimeridian splitting handling.

        For geometries already in WGS84, this is a no-op.
        """
        if self.projection.crs == WGS84_PROJECTION.crs:
            if self.projection == WGS84_PROJECTION:
                return self
            return self.to_projection(WGS84_PROJECTION)
        wgs84 = self.to_projection(WGS84_PROJECTION)
        return split_at_antimeridian(wgs84)

    def intersects(self, other: "STGeometry") -> bool:
        """Returns whether this box intersects the other box."""
        # Check temporal.
        if not self.intersects_time_range(other.time_range):
            return False

        # Check spatial.
        if self.is_global() or other.is_global():
            # One of the geometries indicates global coverage.
            return True

        # If projection is the same, it is an easy check.
        if other.projection == self.projection:
            return shp_intersects(self.shp, other.shp)

        # OK the projections are different. If the CRS is the same, then rescaling is
        # pretty safe so we can just rescale and compare shapely geometries. Otherwise,
        # we check the intersection in WGS84 so we can be sure to handle the
        # anti-meridian properly. Re-projection is required anyway so it shouldn't
        # be any less reliable to re-project to WGS84 (versus re-projecting one
        # geometry to the other's projection), and this allows us to have uniform
        # anti-meridian handling.
        if other.projection.crs == self.projection.crs:
            return shp_intersects(self.shp, other.to_projection(self.projection).shp)
        else:
            return shp_intersects(self.to_wgs84().shp, other.to_wgs84().shp)

    def to_projection(self, projection: Projection) -> "STGeometry":
        """Transforms this geometry to the specified projection.

        Note that this does not handle antimeridian splitting. Use to_wgs84 when
        re-projecting to WGS84 to get antimeridian splitting handling.
        """

        def apply_resolution(
            array: np.ndarray,
            x_resolution: float,
            y_resolution: float,
            forward: bool = True,
        ) -> np.ndarray:
            if forward:
                return np.stack(
                    [array[:, 0] / x_resolution, array[:, 1] / y_resolution], axis=1
                )
            else:
                return np.stack(
                    [array[:, 0] * x_resolution, array[:, 1] * y_resolution], axis=1
                )

        # Undo resolution.
        shp = shapely.transform(
            self.shp,
            lambda array: apply_resolution(
                array,
                self.projection.x_resolution,
                self.projection.y_resolution,
                forward=False,
            ),
        )
        # Change crs.
        # We only apply transform_geom if the CRS doesn't match, because even if we
        # call transform_geom with the same source and destination CRS, it takes
        # several milliseconds.
        if self.projection.crs != projection.crs:
            shp = rasterio.warp.transform_geom(self.projection.crs, projection.crs, shp)
            shp = shapely.geometry.shape(shp)

            # Sometimes the geometry becomes invalid after re-projection, e.g. if there
            # were near-duplicate vertices that crossed over or became the same after
            # transform_geom. So we apply make_valid to fix simple problems.
            if not shp.is_valid:
                shp = shapely.make_valid(shp)

        # Apply new resolution.
        shp = shapely.transform(
            shp,
            lambda array: apply_resolution(
                array, projection.x_resolution, projection.y_resolution, forward=True
            ),
        )

        return STGeometry(projection, shp, self.time_range)

    def __repr__(self) -> str:
        """Returns a string representation of this STGeometry."""
        return (
            f"STGeometry(projection={self.projection}, shp={self.shp}, "
            + f"time_range={self.time_range})"
        )

    def serialize(self) -> dict:
        """Serializes the geometry to a JSON-encodable dictionary."""
        return {
            "projection": self.projection.serialize(),
            "shp": self.shp.wkt,
            "time_range": (
                [self.time_range[0].isoformat(), self.time_range[1].isoformat()]
                if self.time_range
                else None
            ),
        }

    @staticmethod
    def deserialize(d: dict) -> "STGeometry":
        """Deserializes a geometry from a JSON-decoded dictionary."""
        return STGeometry(
            projection=Projection.deserialize(d["projection"]),
            shp=shapely.wkt.loads(d["shp"]),
            time_range=(
                (
                    datetime.fromisoformat(d["time_range"][0]),
                    datetime.fromisoformat(d["time_range"][1]),
                )
                if d["time_range"]
                else None
            ),
        )

contains_time

contains_time(time: datetime) -> bool

Returns whether this box contains the time.

Source code in rslearn/utils/geometry.py

def contains_time(self, time: datetime) -> bool:
    """Returns whether this box contains the time."""
    if self.time_range is None:
        return True
    return time >= self.time_range[0] and time < self.time_range[1]

distance_to_time

distance_to_time(time: datetime) -> timedelta

Returns the distance from this box to the specified time.

Parameters:

Name	Type	Description	Default
time	datetime	the time to compute distance from	required

Returns:

Type	Description
timedelta	the distance, which is 0 if the box contains the time

Source code in rslearn/utils/geometry.py

def distance_to_time(self, time: datetime) -> timedelta:
    """Returns the distance from this box to the specified time.

    Args:
        time: the time to compute distance from

    Returns:
        the distance, which is 0 if the box contains the time
    """
    if self.time_range is None:
        return timedelta()
    if time < self.time_range[0]:
        return self.time_range[0] - time
    if time > self.time_range[1]:
        return time - self.time_range[1]
    return timedelta()

distance_to_time_range

distance_to_time_range(time_range: tuple[datetime, datetime] | None) -> timedelta

Returns the distance from this geometry to the specified time range.

Parameters:

Name	Type	Description	Default
time_range	tuple[datetime, datetime] | None	the time range to compute distance from	required

Returns:

Type	Description
timedelta	the distance, which is 0 if the time ranges intersect

Source code in rslearn/utils/geometry.py

def distance_to_time_range(
    self, time_range: tuple[datetime, datetime] | None
) -> timedelta:
    """Returns the distance from this geometry to the specified time range.

    Args:
        time_range: the time range to compute distance from

    Returns:
        the distance, which is 0 if the time ranges intersect
    """
    if self.time_range is None or time_range is None:
        return timedelta()
    if time_range[1] < self.time_range[0]:
        return self.time_range[0] - time_range[1]
    if self.time_range[1] < time_range[0]:
        return time_range[0] - self.time_range[1]
    return timedelta()

intersects_time_range

intersects_time_range(time_range: tuple[datetime, datetime] | None) -> bool

Returns whether this geometry intersects the other time range.

Source code in rslearn/utils/geometry.py

def intersects_time_range(
    self, time_range: tuple[datetime, datetime] | None
) -> bool:
    """Returns whether this geometry intersects the other time range."""
    if self.time_range is None or time_range is None:
        return True
    if self.time_range[1] <= time_range[0]:
        return False
    if time_range[1] <= self.time_range[0]:
        return False
    return True

is_global

is_global() -> bool

Returns whether this geometry has global spatial coverage.

Global coverage is indicated by a special geometry with WGS84 projection and corners at (-180, -90, 180, 90) (see get_global_geometry).

Source code in rslearn/utils/geometry.py

def is_global(self) -> bool:
    """Returns whether this geometry has global spatial coverage.

    Global coverage is indicated by a special geometry with WGS84 projection and
    corners at (-180, -90, 180, 90) (see get_global_geometry).
    """
    if self.projection != WGS84_PROJECTION:
        return False
    if self.shp != shapely.box(*WGS84_BOUNDS):
        return False
    return True

is_too_large

is_too_large() -> bool

Returns whether this geometry's spatial coverage is too large.

This means that it will likely have issues during re-projections and such.

Source code in rslearn/utils/geometry.py

def is_too_large(self) -> bool:
    """Returns whether this geometry's spatial coverage is too large.

    This means that it will likely have issues during re-projections and such.
    """
    wgs84_bounds = self.to_projection(WGS84_PROJECTION).shp.bounds
    if wgs84_bounds[2] - wgs84_bounds[0] > MAX_GEOMETRY_DEGREES:
        return True
    if wgs84_bounds[3] - wgs84_bounds[1] > MAX_GEOMETRY_DEGREES:
        return True
    return False

to_wgs84

to_wgs84() -> STGeometry

Convert to WGS84 with antimeridian splitting handling.

For geometries already in WGS84, this is a no-op.

Source code in rslearn/utils/geometry.py

def to_wgs84(self) -> "STGeometry":
    """Convert to WGS84 with antimeridian splitting handling.

    For geometries already in WGS84, this is a no-op.
    """
    if self.projection.crs == WGS84_PROJECTION.crs:
        if self.projection == WGS84_PROJECTION:
            return self
        return self.to_projection(WGS84_PROJECTION)
    wgs84 = self.to_projection(WGS84_PROJECTION)
    return split_at_antimeridian(wgs84)

intersects

intersects(other: STGeometry) -> bool

Returns whether this box intersects the other box.

Source code in rslearn/utils/geometry.py

def intersects(self, other: "STGeometry") -> bool:
    """Returns whether this box intersects the other box."""
    # Check temporal.
    if not self.intersects_time_range(other.time_range):
        return False

    # Check spatial.
    if self.is_global() or other.is_global():
        # One of the geometries indicates global coverage.
        return True

    # If projection is the same, it is an easy check.
    if other.projection == self.projection:
        return shp_intersects(self.shp, other.shp)

    # OK the projections are different. If the CRS is the same, then rescaling is
    # pretty safe so we can just rescale and compare shapely geometries. Otherwise,
    # we check the intersection in WGS84 so we can be sure to handle the
    # anti-meridian properly. Re-projection is required anyway so it shouldn't
    # be any less reliable to re-project to WGS84 (versus re-projecting one
    # geometry to the other's projection), and this allows us to have uniform
    # anti-meridian handling.
    if other.projection.crs == self.projection.crs:
        return shp_intersects(self.shp, other.to_projection(self.projection).shp)
    else:
        return shp_intersects(self.to_wgs84().shp, other.to_wgs84().shp)

to_projection

to_projection(projection: Projection) -> STGeometry

Transforms this geometry to the specified projection.

Note that this does not handle antimeridian splitting. Use to_wgs84 when re-projecting to WGS84 to get antimeridian splitting handling.

Source code in rslearn/utils/geometry.py

def to_projection(self, projection: Projection) -> "STGeometry":
    """Transforms this geometry to the specified projection.

    Note that this does not handle antimeridian splitting. Use to_wgs84 when
    re-projecting to WGS84 to get antimeridian splitting handling.
    """

    def apply_resolution(
        array: np.ndarray,
        x_resolution: float,
        y_resolution: float,
        forward: bool = True,
    ) -> np.ndarray:
        if forward:
            return np.stack(
                [array[:, 0] / x_resolution, array[:, 1] / y_resolution], axis=1
            )
        else:
            return np.stack(
                [array[:, 0] * x_resolution, array[:, 1] * y_resolution], axis=1
            )

    # Undo resolution.
    shp = shapely.transform(
        self.shp,
        lambda array: apply_resolution(
            array,
            self.projection.x_resolution,
            self.projection.y_resolution,
            forward=False,
        ),
    )
    # Change crs.
    # We only apply transform_geom if the CRS doesn't match, because even if we
    # call transform_geom with the same source and destination CRS, it takes
    # several milliseconds.
    if self.projection.crs != projection.crs:
        shp = rasterio.warp.transform_geom(self.projection.crs, projection.crs, shp)
        shp = shapely.geometry.shape(shp)

        # Sometimes the geometry becomes invalid after re-projection, e.g. if there
        # were near-duplicate vertices that crossed over or became the same after
        # transform_geom. So we apply make_valid to fix simple problems.
        if not shp.is_valid:
            shp = shapely.make_valid(shp)

    # Apply new resolution.
    shp = shapely.transform(
        shp,
        lambda array: apply_resolution(
            array, projection.x_resolution, projection.y_resolution, forward=True
        ),
    )

    return STGeometry(projection, shp, self.time_range)

serialize

serialize() -> dict

Serializes the geometry to a JSON-encodable dictionary.

Source code in rslearn/utils/geometry.py

def serialize(self) -> dict:
    """Serializes the geometry to a JSON-encodable dictionary."""
    return {
        "projection": self.projection.serialize(),
        "shp": self.shp.wkt,
        "time_range": (
            [self.time_range[0].isoformat(), self.time_range[1].isoformat()]
            if self.time_range
            else None
        ),
    }

deserialize staticmethod

deserialize(d: dict) -> STGeometry

Deserializes a geometry from a JSON-decoded dictionary.

Source code in rslearn/utils/geometry.py

@staticmethod
def deserialize(d: dict) -> "STGeometry":
    """Deserializes a geometry from a JSON-decoded dictionary."""
    return STGeometry(
        projection=Projection.deserialize(d["projection"]),
        shp=shapely.wkt.loads(d["shp"]),
        time_range=(
            (
                datetime.fromisoformat(d["time_range"][0]),
                datetime.fromisoformat(d["time_range"][1]),
            )
            if d["time_range"]
            else None
        ),
    )

is_same_resolution

is_same_resolution(res1: float, res2: float) -> bool

Returns whether the two resolutions are the same.

Source code in rslearn/utils/geometry.py

def is_same_resolution(res1: float, res2: float) -> bool:
    """Returns whether the two resolutions are the same."""
    return (max(res1, res2) / min(res1, res2) - 1) < RESOLUTION_EPSILON

shp_intersects

shp_intersects(shp1: Geometry, shp2: Geometry) -> bool

Returns whether the two shapes intersect.

Tries shp.intersects but falls back to shp.intersection which can be more reliable.

Source code in rslearn/utils/geometry.py

def shp_intersects(shp1: shapely.Geometry, shp2: shapely.Geometry) -> bool:
    """Returns whether the two shapes intersect.

    Tries shp.intersects but falls back to shp.intersection which can be more
    reliable.
    """
    try:
        return shp1.intersects(shp2)
    except shapely.GEOSException:
        return shp1.intersection(shp2).area > 0

get_global_raster_bounds

get_global_raster_bounds(projection: Projection) -> PixelBounds

Get very large pixel bounds for a global raster in the given projection.

This is useful for data sources that cover the entire world and don't want to compute exact bounds in arbitrary projections (which can fail for projections like UTM that only cover part of the world).

Parameters:

Name	Type	Description	Default
projection	Projection	the projection to get bounds in.	required

Returns:

Type	Description
PixelBounds	Pixel bounds that will intersect with any reasonable window. We assume that the
PixelBounds	absolute value of CRS coordinates is at most 2^32, and adjust it based on the
PixelBounds	resolution in the Projection in case very fine-grained resolutions are used.

Source code in rslearn/utils/geometry.py

def get_global_raster_bounds(projection: Projection) -> PixelBounds:
    """Get very large pixel bounds for a global raster in the given projection.

    This is useful for data sources that cover the entire world and don't want to
    compute exact bounds in arbitrary projections (which can fail for projections
    like UTM that only cover part of the world).

    Args:
        projection: the projection to get bounds in.

    Returns:
        Pixel bounds that will intersect with any reasonable window. We assume that the
        absolute value of CRS coordinates is at most 2^32, and adjust it based on the
        resolution in the Projection in case very fine-grained resolutions are used.
    """
    crs_bound = 2**32
    pixel_bound_x = int(crs_bound / abs(projection.x_resolution))
    pixel_bound_y = int(crs_bound / abs(projection.y_resolution))
    return (-pixel_bound_x, -pixel_bound_y, pixel_bound_x, pixel_bound_y)

get_global_geometry

get_global_geometry(time_range: tuple[datetime, datetime] | None) -> STGeometry

Gets a geometry that indicates global spatial coverage for the given time range.

Parameters:

Name	Type	Description	Default
time_range	tuple[datetime, datetime] | None	the time range for the STGeometry.	required

Returns:

Type	Description
STGeometry	STGeometry with global spatial coverage and specified time range.

Source code in rslearn/utils/geometry.py

def get_global_geometry(time_range: tuple[datetime, datetime] | None) -> STGeometry:
    """Gets a geometry that indicates global spatial coverage for the given time range.

    Args:
        time_range: the time range for the STGeometry.

    Returns:
        STGeometry with global spatial coverage and specified time range.
    """
    return STGeometry(WGS84_PROJECTION, shapely.box(*WGS84_BOUNDS), time_range)

flatten_shape

flatten_shape(shp: Geometry) -> list[Geometry]

Flatten the shape into a list of primitive shapes (Point, LineString, and Polygon).

Parameters:

Name	Type	Description	Default
shp	Geometry	the shape, which could be a primitive shape like polygon or a collection.	required

Returns:

Type	Description
list[Geometry]	list of primitive shapes.

Source code in rslearn/utils/geometry.py

def flatten_shape(shp: shapely.Geometry) -> list[shapely.Geometry]:
    """Flatten the shape into a list of primitive shapes (Point, LineString, and Polygon).

    Args:
        shp: the shape, which could be a primitive shape like polygon or a collection.

    Returns:
        list of primitive shapes.
    """
    if isinstance(
        shp,
        shapely.MultiPoint
        | shapely.MultiLineString
        | shapely.MultiPolygon
        | shapely.GeometryCollection,
    ):
        flat_list: list[shapely.Geometry] = []
        for component in shp.geoms:
            flat_list.extend(flatten_shape(component))
        return flat_list

    else:
        return [shp]

split_shape_at_antimeridian

split_shape_at_antimeridian(shp: Geometry, epsilon: float = 1e-06) -> Geometry

Split the given shape at the antimeridian.

The shape must be in WGS84 coordinates.

See split_at_antimeridian for details.

Parameters:

Name	Type	Description	Default
shp	Geometry	the shape to split.	required
epsilon	float	the padding in degrees.	1e-06

Returns:

Type	Description
Geometry	the split shape, in WGS84 projection.

Source code in rslearn/utils/geometry.py

def split_shape_at_antimeridian(
    shp: shapely.Geometry, epsilon: float = 1e-6
) -> shapely.Geometry:
    """Split the given shape at the antimeridian.

    The shape must be in WGS84 coordinates.

    See split_at_antimeridian for details.

    Args:
        shp: the shape to split.
        epsilon: the padding in degrees.

    Returns:
        the split shape, in WGS84 projection.
    """
    # We assume the shape is fine if:
    # 1. It doesn't need padding (no coordinates close to +/- 180).
    # 2. And all coordinates are either less than 90 or more than -90 (meaning the
    #    shape approaches the antimeridian on at most one side).
    bounds = shp.bounds
    if bounds[0] > -180 + epsilon and bounds[2] < 90:
        return shp
    if bounds[0] > -90 and bounds[2] < 180 - epsilon:
        return shp

    if isinstance(
        shp,
        shapely.MultiPoint
        | shapely.MultiLineString
        | shapely.MultiPolygon
        | shapely.GeometryCollection,
    ):
        return _collect_shapes(
            [split_shape_at_antimeridian(component) for component in shp.geoms]
        )

    if isinstance(shp, shapely.Point):
        # Points only need padding.
        lon = shp.x
        if lon < -180 + epsilon:
            lon = -180 + epsilon
        if lon > 180 - epsilon:
            lon = 180 - epsilon
        return shapely.Point(lon, shp.y)

    if isinstance(shp, shapely.LineString | shapely.Polygon):
        # We add 360 to the negative coordinates and then separate the parts above and
        # below 180.
        def add360(array: npt.NDArray[np.float32]) -> npt.NDArray[np.float32]:
            new_array = array.copy()
            new_array[new_array[:, 0] < 0, 0] += 360
            return new_array

        shp = shapely.transform(shp, add360)

        positive_part = shapely.box(0, -90, 180 - epsilon, 90)
        negative_part = shapely.box(180 + epsilon, -90, 360, 90)
        positive_shp = shp.intersection(positive_part)
        negative_shp = shp.intersection(negative_part)
        negative_shp = shapely.transform(negative_shp, lambda coords: coords - [360, 0])
        return _collect_shapes([positive_shp, negative_shp])

    raise TypeError("Unsupported shape type")

split_at_antimeridian

split_at_antimeridian(geometry: STGeometry, epsilon: float = 1e-06) -> STGeometry

Split lines and polygons in the given geometry at the antimeridian.

The returned geometry will always be in WGS84 projection.

Small padding is also introduced to ensure coordinates are a bit more than -180 or a bit less than 180.

For example, if the input is a polygon:

Polygon([[-180, 10], [180, 11], [-179, 11], [-179, 10]])

Then it would be converted to:

Polygon([[-179.999999, 10], [-179,999999, 11], [-179, 11], [-179, 10]])

This function may produce unexpected results if the geometries span more than 90 degrees on either dimension.

Parameters:

Name	Type	Description	Default
geometry	STGeometry	the geometry to split.	required
epsilon	float	the padding in degrees. It is equivalent to about 1 m at the equator. We ensure no longitude coordinates are within this padding of +/- 180.	1e-06

Returns:

Type	Description
STGeometry	the padded geometry, in WGS84 projection.

Source code in rslearn/utils/geometry.py

def split_at_antimeridian(geometry: STGeometry, epsilon: float = 1e-6) -> STGeometry:
    """Split lines and polygons in the given geometry at the antimeridian.

    The returned geometry will always be in WGS84 projection.

    Small padding is also introduced to ensure coordinates are a bit more than -180 or
    a bit less than 180.

    For example, if the input is a polygon:

        Polygon([[-180, 10], [180, 11], [-179, 11], [-179, 10]])

    Then it would be converted to:

        Polygon([[-179.999999, 10], [-179,999999, 11], [-179, 11], [-179, 10]])

    This function may produce unexpected results if the geometries span more than 90
    degrees on either dimension.

    Args:
        geometry: the geometry to split.
        epsilon: the padding in degrees. It is equivalent to about 1 m at the equator.
            We ensure no longitude coordinates are within this padding of +/- 180.

    Returns:
        the padded geometry, in WGS84 projection.
    """
    # Convert to WGS84.
    geometry = geometry.to_projection(WGS84_PROJECTION)
    new_shp = split_shape_at_antimeridian(geometry.shp, epsilon=epsilon)
    return STGeometry(geometry.projection, new_shp, geometry.time_range)

safely_reproject_within_valid_area

safely_reproject_within_valid_area(src_geoms: Sequence[STGeometry], valid_geom: STGeometry) -> Sequence[STGeometry | None]

Re-project src_geoms into the projection of valid_geom.

Unlike direct to_projection(), this clips each source geometry in WGS84 to a buffered area around valid_geom before reprojecting. This minimizes distortions in case valid_geom is small but src_geoms may be large. It works best if src_geoms are also natively in WGS84; otherwise, there could be distortion issues re-projecting them to WGS84.

Returns a list with one entry per source geometry. The entry is either the re-projected geometry, or it may be (but is not guaranteed to be) None if the source geometry doesn't intersect valid_geom.

Source code in rslearn/utils/geometry.py

def safely_reproject_within_valid_area(
    src_geoms: Sequence[STGeometry], valid_geom: STGeometry
) -> Sequence[STGeometry | None]:
    """Re-project src_geoms into the projection of valid_geom.

    Unlike direct to_projection(), this clips each source geometry in WGS84 to a
    buffered area around valid_geom before reprojecting. This minimizes distortions in
    case valid_geom is small but src_geoms may be large. It works best if src_geoms are
    also natively in WGS84; otherwise, there could be distortion issues re-projecting
    them to WGS84.

    Returns a list with one entry per source geometry. The entry is either the
    re-projected geometry, or it may be (but is not guaranteed to be) None if the
    source geometry doesn't intersect valid_geom.
    """

    @functools.cache
    def get_valid_geom_wgs84() -> STGeometry:
        return valid_geom.to_wgs84()

    @functools.cache
    def get_clip_region() -> shapely.Geometry:
        """Buffered WGS84 region around valid_geom for clipping large sources.

        Buffers each component of the geometry individually and takes the union (This
        avoids issues with MultiPolygons coming from antimeridian splitting, where
        using shp.bounds directly would result in a big box).
        """
        shp = get_valid_geom_wgs84().shp
        parts = shp.geoms if hasattr(shp, "geoms") else [shp]
        boxes = []
        for part in parts:
            b = part.bounds
            # We buffer it by WGS84_CLIP_BUFFER_DEGREES, plus half its size (so we
            # minimize issues for larger windows).
            half_size = max(b[2] - b[0], b[3] - b[1]) / 2
            boxes.append(
                shapely.box(
                    b[0] - WGS84_CLIP_BUFFER_DEGREES - half_size,
                    b[1] - WGS84_CLIP_BUFFER_DEGREES - half_size,
                    b[2] + WGS84_CLIP_BUFFER_DEGREES + half_size,
                    b[3] + WGS84_CLIP_BUFFER_DEGREES + half_size,
                )
            )
        return shapely.unary_union(boxes)

    results: list[STGeometry | None] = []
    for src_geom in src_geoms:
        if src_geom.projection == valid_geom.projection:
            results.append(src_geom)
            continue
        # If the CRS is the same, to_projection is just rescaling which is safe.
        if src_geom.projection.crs == valid_geom.projection.crs:
            results.append(src_geom.to_projection(valid_geom.projection))
            continue

        src_wgs84 = src_geom.to_wgs84()
        valid_wgs84 = get_valid_geom_wgs84()

        if not shp_intersects(src_wgs84.shp, valid_wgs84.shp):
            results.append(None)
            continue

        # Clip in WGS84 to the buffered valid area so we only reproject a small piece.
        clipped_shp = src_wgs84.shp.intersection(get_clip_region())
        if clipped_shp.is_empty:
            results.append(None)
            continue

        clipped_geom = STGeometry(WGS84_PROJECTION, clipped_shp, src_geom.time_range)
        reprojected = clipped_geom.to_projection(valid_geom.projection)
        if reprojected.shp.is_empty:
            results.append(None)
            continue

        results.append(reprojected)

    return results