Working with Geometries

Geometries are a common cause of friction in geospatial data processing. This is especially true when working with geometries that cross the antimeridian or cover a pole. Tilebox is designed to take away most of the friction involved in this, but it’s still recommended to follow the best practices for handling geometries outlined below.

In doing so, you can ensure that no geometry related issues will arise even when interfacing with other libraries and tools that may not properly support non-linearities in geometries.

​

Best Practices

To ensure expected behavior when working with geometries, it’s recommended to follow these best practices:

1. Use language specific geometry libraries for creating, parsing and serializing geometries.
2. Polygon exteriors should always be defined in counter-clockwise winding order.
3. Polygon holes (interior rings) should always be defined in clockwise winding order.
4. Longitude values should always be within the [-180, 180] range, and latitude values should always be within the [-90, 90] range.
5. Geometries that cross the antimeridian should be cut along the antimeridian into two parts—one for the eastern hemisphere and one for the western hemisphere.
6. Geometries that cover a pole need to be defined in a specific way to ensure correct handling.
7. When downloading geometries from external sources, such as from STAC Catalogues always make sure to verify that those assumptions are met to avoid unexpected behavior down the road.

​

Geometry vs Geography

Some systems expose two similar, but different data types related to geometries: Geometry and Geography. The main difference lies in how edge interpolation along the antimeridian is handled. Geometry represents a geometry in an arbitrary 2D cartesian coordinate system, and does not perform any edge interpolation, while Geography typically does wrap around the antimeridian by performing edge interpolation, such as converting cartesian coordinates to a 3D spherical coordinate system. For Geography types, the x coordinate is typically limited to the [-180, 180] range, and the y coordinate is limited to the [-90, 90] range, while for Geometry types no such limitations are imposed.

Tilebox does not make such a distinction between Geometry and Geography types. Instead, it provides a query option to specify a coordinate reference system to control whether geometry intersections and containment checks are performed in a 3D spherical coordinate system (which correctly handles antimeridian crossings) or in a standard 2D cartesian lat/lon coordinate system.

​

Terminology

Get familiar with some key concepts related to geometries.

Point

A Point is a specific location on the Earth’s surface, represented by a longitude and latitude coordinate.

Ring

A Ring is a collection of points that form a closed loop. The order of the points within the ring matter, since that defines the winding order of the ring, which is either clockwise or counter-clockwise. Every ring should be explicitly closed, meaning that the first and last point should be the same.

Polygon

A Polygon is a collection of rings. The first ring is called the exterior ring, and represents the boundary of the polygon. Any other rings are called interior rings, and represent holes in the polygon.

MultiPolygon

A MultiPolygon is a collection of polygons.

​

Geometry libraries and common formats

Geometries can be expressed in many different ways, and many formats exist for representing geometries. To handle these, the Tilebox Client SDKs delegate geometry handling to external, well-known geometry libraries.

Here is an example of how to define a Polygon geometry, covering the area of the state of Colorado using Tilebox Client SDKs.

from shapely import Polygon

area = Polygon(
    ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)),
)

Through these libraries, Tilebox also supports some of the most common formats for representing geometries.

​

GeoJSON

GeoJSON is a geospatial data interchange format based on JavaScript Object Notation (JSON).

{
    "type": "Polygon",
    "coordinates": [
        [
            [-109.05, 41.0],
            [-109.045, 37.0],
            [-102.05, 37.0],
            [-102.05, 41.0],
            [-109.05, 41.0]
        ]
    ]
}

Reading such a GeoJSON geometry can be achieved as follows.

from shapely import from_geojson

with open("colorado.geojson", "r") as f:
    area = from_geojson(json.load(f))

​

Well-Known Text (WKT)

Well-Known Text (WKT) is a text markup language for representing vector geometry objects.

from shapely import from_wkt

wkt = "POLYGON ((-109.05 41, -109.045 37, -102.05 37, -102.05 41, -109.05 41))"
area = from_wkt(wkt)

​

Well-Known Binary (WKB)

Well-Known Binary (WKB) is a binary representation of vector geometry objects, a binary equivalent to WKT.

Given a colorado.wkb file containing a binary encoding of the Colorado polygon, it can be read as follows.

from shapely import from_wkb

with open("colorado.wkb", "rb") as f:
    area = from_wkb(f.read())

​

Winding Order

The winding order of a ring defines the orientation of the ring, and is either clockwise or counter-clockwise. The winding order of a ring is determined by the order of the points within the ring.

Geometries in Tilebox follow the right hand rule defined by the GeoJSON specification:

A linear ring MUST follow the right-hand rule with respect to the area it bounds, i.e., exterior rings are counterclockwise, and holes are clockwise.

This means that the exterior ring of a polygon has a counter-clockwise winding order, while interior rings have a clockwise winding order.

When this rule is not followed, querying may yield unexpected results, as can be seen in the example below.

Take the following Polygon consisting of the same exterior ring but in different winding orders.

POLYGON (
  (-5 56, -6 29, 14 28, 23 55, -5 56)
)

POLYGON (
  (-5 56, 23 55, 14 28, -6 29, -5 56)
)

This is how those two geometries would be interpreted on a sphere.

Since such unexpected behavior can cause issues and be hard to debug, Tilebox does detect incorrect winding orders and automatically fixes them when ingesting data. This means that Geometries that really cover almost the whole globe except a small area cannot be ingested into Tilebox by simply specifying them as a Polygon with an exterior ring in clockwise winding order.

Instead, such geometries should be defined as a Polygon consisting of two rings:

• an exterior ring covering the whole globe in counter-clockwise winding order
• an interior ring specifying the hole in clockwise winding order

Such a definition, which Tilebox will interpret as intended, and also is valid according to the GeoJSON specification, is shown below.

POLYGON (
  (-180 90, -180 -90, 180 -90, 180 90, -180 90),
  (-5 56, 23 55, 14 28, -6 29, -5 56)
)

To verify the winding order of a geometry, you can use the following code snippets.

from shapely import Polygon

polygon = Polygon(
    ((-109.05, 41.00), (-109.045, 37.0), (-102.05, 37.0), (-102.05, 41.00), (-109.05, 41.00)),
)
print(polygon.exterior.is_ccw)  # True

​

Antimeridian Crossings

Geometries that cross the antimeridian are a common occurrence in satellite data, but can cause issues when not handled correctly.

Take the following Polygon that crosses the antimeridian.

Below are a couple of different ways to express such a geometry.

POLYGON ((173 5, -175 3, -172 20, 177 23, 173 5))

While this is a valid representation, it causes issues for a lot of libraries that do calculations or visualization in a cartesian coordinate system, since the line from lon=173 to lon=-175 will be interpreted as a line that is 348 degrees long, crossing the null meridian and almost all the way around the globe. Often times, visualizations of such geometries will look like in the image below.

To avoid such issues, some tools working in cartesian coordinate space sometimes extend the possible longitude range beyond to usual bounds of [-180, 180]. Expressing the geometry in such a way, may look like this.

POLYGON ((173 5, 185 3, 188 20, 177 23, 173 5))

Or, the exact same area on the globe can also be expressed as a geometry by extending the longitude range below -180.

POLYGON ((-187 5, -175 3, -172 20, -183 23, -187 5))

While most visualization tools will probably handle such geometries correctly, special care needs to be taken when working with such geometries when it comes to intersection and containment checks.

The below code snippet illustrates this, by constructing a Polygon that covers the exact same area in the two different methods shown, and checking whether they intersect each other (which they should, since they represent the same area on the globe).

from shapely import from_wkt

a = from_wkt("POLYGON ((173 5, 185 3, 188 20, 177 23, 173 5))")
b = from_wkt("POLYGON ((-187 5, -175 3, -172 20, -183 23, -187 5))")

print(a.intersects(b))  # False

​

Antimeridian Cutting

Additionally, none of the representations shown are valid according to the GeoJSON specification. The GeoJSON specification does offers a solution for this problem, though, which is Antimeridian Cutting. It suggests always cutting lines and polygons that cross the antimeridian into two parts—one for the eastern hemisphere and one for the western hemisphere.

To achieve this, check out the antimeridian python package, and an implementation of it in Go.

MULTIPOLYGON (
    ((180 22.265994, 177 23, 173 5, 180 3.855573, 180 22.265994)),
    ((-180 3.855573, -175 3, -172 20, -180 22.265994, -180 3.855573))
)

By cutting geometries along the antimeridian:

• they conform to the GeoJSON specification
• there is no ambiguity about how to correctly interpret the geometry
• they are correctly handled by most visualization tools
• intersection and containment checks yield correct results no matter the coordinate reference system used
• all longitude values are within the valid [-180, 180] range, making it easier to work with them

​

Pole Coverings

Geometries that cover a pole can be especially tricky to handle correctly. No one algorithm that can handle all possible cases correctly, and different libraries and tools may handle them differently. A lot of times, injecting prior knowledge, such as the fact that a geometry is supposed to cover a pole, can help to resolve such issues. For an example of this, check out the relevant section in the antimeridian documentation.

Generally speaking though, there are two approaches that work well:

​

Approach 1: Cutting out a hole

Define a geometry with an exterior ring that covers the whole globe. Then cut out a hole by defining an interior ring of the area that is not covered by the geometry.

POLYGON (
    (-180 90, -180 -90, 180 -90, 180 90, -180 90), 
    (-133 72, 153 77, 130 -64, -160 -66, -133 72)
)

​

Approach 2: Splitting into multiple parts

Another approach, that works well for circular caps covering a pole, involves cutting the geometry into multiple parts, but instead of splitting the geometry only at the antimeridian, also splitting it at the null meridian, the 90 and -90 meridians. For a circular cap covering the north pole, this results in four triangular parts, which can be combined into a MultiPolygon. Visualization libraries as well as intersection and containment checks done in cartesian coordinate space will then typically handle such a geometry correctly.

MULTIPOLYGON (
    ((-90 75, 0 75, 0 90, -90 90, -90 75)),
    ((0 75, 90 75, 90 90, 0 90, 0 75)),
    ((90 75, 180 75, 180 90, 90 90, 90 75)),
    ((-180 75, -90 75, -90 90, -180 90, -180 75))
)