Many datasets consist of granules that represent specific geographical areas on the Earth’s surface. Often, a polygon defining the outline of these areas—a footprint—accompanies other granule metadata in time series datasets. Tilebox provides built-in support for working with geometries.
Here’s an example that loads some granules from the ERS SAR Opendata dataset, which contains geometries.
from tilebox.datasets import Client
client = Client()
datasets = client.datasets()
ers_collection = datasets.open_data.asf.ers_sar.collection("ERS-2")
ers_data = ers_collection.load(("2008-02-10T21:00", "2008-02-10T22:00"))
Shapely
In the ers_data dataset, each granule includes a geometry field that represents the footprint of each granule as a polygon. Tilebox automatically converts geometry fields to Polygon or MultiPolygon objects from the Shapely library. By integrating with Shapely, you can use the rich set of libraries and tools it provides. That includes support for computing polygon characteristics such as total area, intersection checks, and conversion to other formats.
geometries = ers_data.geometry.values
print(geometries)
Each geometry is a shapely.Geometry.
[<POLYGON ((-150.753 74.25, -152.032 73.336, -149.184 73.002, -147.769 73.899...>
<POLYGON ((-151.29 73.878, -152.517 72.962, -149.722 72.634, -148.362 73.535...>
<POLYGON ((-152.281 73.146, -153.418 72.224, -150.721 71.908, -149.46 72.814...>
<POLYGON ((-153.218 72.393, -154.273 71.467, -151.672 71.161, -150.5 72.073,...>
<POLYGON ((-154.073 71.65, -155.056 70.721, -152.542 70.424, -151.449 71.341...>
<POLYGON ((-154.776 70.899, -155.699 69.967, -153.242 69.675, -152.214 70.59...>
<POLYGON ((-155.527 70.147, -156.393 69.212, -154.015 68.928, -153.05 69.852...>
<POLYGON ((-156.238 69.385, -157.052 68.447, -154.749 68.171, -153.842 69.09...>
<POLYGON ((-156.896 68.633, -157.664 67.693, -155.431 67.423, -154.575 68.35...>
<POLYGON ((-157.522 67.871, -158.248 66.929, -156.082 66.666, -155.273 67.6,...>
<POLYGON ((-158.113 67.109, -158.801 66.165, -156.697 65.908, -155.932 66.84...>
<POLYGON ((-158.628 66.348, -159.283 65.402, -157.226 65.149, -156.497 66.08...>
<POLYGON ((-159.669 64.814, -160.263 63.866, -158.317 63.624, -157.656 64.56...>]
Geometries are not always Polygon objects. More complex footprint geometries are represented as MultiPolygon objects.
Accessing Coordinates
You can select a polygon from the geometries and access the underlying coordinates and an automatically computed centroid point.
Accessing coordinates and computing a centroid point
polygon = geometries[0]
lon, lat = polygon.exterior.coords.xy
center, = list(polygon.centroid.coords)
print(lon)
print(lat)
print(center)
array('d', [-150.753244, -152.031574, -149.183655, -147.769339, -150.753244])
array('d', [74.250081, 73.336051, 73.001748, 73.899483, 74.250081])
(-149.92927907414239, 73.62538063474753)
Interactive environments such as Jupyter Notebooks can visualize Polygon shapes graphically. Just type polygon in an empty cell and execute it for a visual representation of the polygon shape.
Visualization on a Map
To visualize polygons on a map, you can use Folium. Below is a helper function that produces an OpenStreetMap with the Polygon overlaid.
visualize helper function
# pip install folium
from folium import Figure, Map, Polygon as FoliumPolygon, GeoJson, TileLayer
from folium.plugins import MiniMap
from shapely import Polygon, to_geojson
from collections.abc import Iterable
def visualize(poly: Polygon | Iterable[Polygon], zoom=4):
"""Visualize a polygon or a list of polygons on a map"""
if not isinstance(poly, Iterable):
poly = [poly]
fig = Figure(width=600, height=600)
map = Map(location=geometries[len(geometries)//2].centroid.coords[0][::-1], zoom_start=zoom, control_scale=True)
map.add_child(MiniMap())
fig.add_child(map)
for p in poly:
map.add_child(GeoJson(to_geojson(p)))
return fig
Here’s how to use it.
The visualize helper function supports a list of polygons, which can display the data layout of the ERS granules.
Visualizing multiple polygons
Shapely supports converting Polygons to some common formats, such as GeoJSON or Well-Known Text (WKT).
from shapely import to_geojson
print(to_geojson(polygon))
{"type":"Polygon","coordinates":[[[-150.753244,74.250081],[-152.031574,73.336051],[-149.183655,73.001748],[-147.769339,73.899483],[-150.753244,74.250081]]]}
from shapely import to_wkt
print(to_wkt(polygon))
POLYGON ((-150.753244 74.250081, -152.031574 73.336051, -149.183655 73.001748, -147.769339 73.899483, -150.753244 74.250081))
Checking intersections
One common task when working with geometries is checking if a given geometry falls into a specific area of interest. Shapely provides an intersects method for this purpose.
from shapely import box
# Box representing the rectangular area lon=(-160, -150) and lat=(69, 70)
area_of_interest = box(-160, 69, -150, 70)
for i, polygon in enumerate(geometries):
if area_of_interest.intersects(polygon):
print(f"{ers_data.granule_name[i].item()} intersects the area of interest!")
else:
print(f"{ers_data.granule_name[i].item()} doesn't intersect the area of interest!")
E2_66974_STD_F264 doesn't intersect the area of interest!
E2_66974_STD_F265 doesn't intersect the area of interest!
E2_66974_STD_F267 doesn't intersect the area of interest!
E2_66974_STD_F269 doesn't intersect the area of interest!
E2_66974_STD_F271 doesn't intersect the area of interest!
E2_66974_STD_F273 intersects the area of interest!
E2_66974_STD_F275 intersects the area of interest!
E2_66974_STD_F277 intersects the area of interest!
E2_66974_STD_F279 doesn't intersect the area of interest!
E2_66974_STD_F281 doesn't intersect the area of interest!
E2_66974_STD_F283 doesn't intersect the area of interest!
E2_66974_STD_F285 doesn't intersect the area of interest!
E2_66974_STD_F289 doesn't intersect the area of interest!
Combining polygons
As shown in the visualization of the granule footprints, the granules collectively form an orbit from pole to pole. Measurements are often combined during processing. You can do the same with geometries by combining them into a single polygon, which represents the hull around all individual footprints using shapely.unary_union.
Combining multiple polygons
from shapely.ops import unary_union
hull = unary_union(geometries)
visualize(hull)
The computed hull consists of two polygons due to a gap (probably a missing granule) in the geometries. Such geometries are represented as Multi Polygons.
Multi Polygons
A collection of one or more non-overlapping polygons combined into a single geometry is called a MultiPolygon. Footprint geometries can be of type MultiPolygon due to gaps or pole discontinuities. The computed hull in the previous example is a MultiPolygon.
Accessing individual polygons of a MultiPolygon
print(f"The computed hull of type {type(hull).__name__} consists of {len(hull.geoms)} sub polygons")
for i, poly in enumerate(hull.geoms):
print(f"Sub polygon {i} has an area of {poly.area}")
The computed hull of type MultiPolygon consists of 2 sub polygons
Sub polygon 0 has an area of 2.025230449898011
Sub polygon 1 has an area of 24.389998081651527
Antimeridian Crossings
A common issue with longitude / latitude geometries is crossings of the 180-degree meridian, or the antimeridian. For example, the coordinates of a LineString from Japan to the United States might look like this:
140, 141, 142, ..., 179, 180, -179, -178, ..., -125, -124
Libraries like Shapely are not designed to handle spherical coordinate systems, so caution is necessary with such geometries.
Here’s an ERS granule demonstrating this issue.
# A granule that crosses the antimeridian
granule = ers_collection.find("0119bb86-0260-5819-6aab-f99796417155")
polygon = granule.geometry.item()
print(polygon.exterior.coords.xy)
visualize(polygon)
array('d', [177.993407, 176.605009, 179.563047, -178.904076, 177.993407])
array('d', [74.983185, 74.074615, 73.727752, 74.61847, 74.983185])
This 2D visualization appears incorrect. Both the visualization and any calculations performed may yield inaccurate results. For instance, testing whether the granule intersects the 0-degree meridian provides a false positive.
Problems with calculating intersections
from shapely import LineString
null_meridian = LineString([(0, -90), (0, 90)])
print(polygon.intersects(null_meridian)) # True - but this is incorrect!
The GeoJSON specification offers a solution for this problem. In the section Antimeridian Cutting, it suggests always cutting lines and polygons into two parts—one for the eastern hemisphere and one for the western hemisphere.
In python, this can be achieved using the antimeridian package.
Cutting the polygon along the antimeridian
# pip install antimeridian
import antimeridian
fixed_polygon = antimeridian.fix_polygon(polygon)
visualize(fixed_polygon)
print(fixed_polygon.intersects(null_meridian)) # False - this is correct now
Since Shapely is unaware of the spherical nature of this data, the centroid of the fixed polygon is still incorrect. The antimeridian package also includes a function to correct this.
Calculating the centroid of a cut polygon crossing the antimeridian
print("Wrongly computed centroid coordinates (Shapely)")
print(list(fixed_polygon.centroid.coords))
print("Correct centroid coordinates (Antimeridian taken into account)")
print(list(antimeridian.centroid(fixed_polygon).coords))
Wrongly computed centroid coordinates (shapely)
[(139.8766350146937, 74.3747116658462)]
Correct centroid coordinates (antimeridian taken into account)
[(178.7782777050171, 74.3747116658462)]
Spherical Geometry
Another approach to handle the antimeridian issue is performing all coordinate-related calculations, such as polygon intersections, in a spherical coordinate system.
One useful library for this is spherical_geometry. Here’s an example.
# pip install spherical-geometry
from spherical_geometry.polygon import SphericalPolygon
from spherical_geometry.vector import lonlat_to_vector
lon, lat = polygon.exterior.coords.xy
spherical_poly = SphericalPolygon.from_lonlat(lon, lat)
# Let's check the x, y, z coordinates of the spherical polygon:
print(list(spherical_poly.points))
[array([[-0.25894363, 0.00907234, 0.96584983],
[-0.2651968 , -0.00507317, 0.96418096],
[-0.28019363, 0.00213687, 0.95994112],
[-0.27390375, 0.01624885, 0.96161984],
[-0.25894363, 0.00907234, 0.96584983]])]
Now, you can compute intersections or check if a particular point is within the polygon. You can compare the incorrect calculation using shapely with the correct version when using spherical_geometry.
Correct calculations using spherical geometry
# A point on the null-meridian, way off from our polygon
null_meridian_point = 0, 74.4
# A point actually inside our polygon
point_inside = 178.8, 74.4
print("Shapely results:")
print("- Null meridian point inside:", polygon.contains(Point(*null_meridian_point)))
print("- Actual inside point inside:", polygon.contains(Point(*point_inside)))
print("Spherical geometry results:")
print("- Null meridian point inside:", spherical_poly.contains_lonlat(*null_meridian_point))
print("- Actual inside point inside:", spherical_poly.contains_lonlat(*point_inside))
Shapely results:
- Null meridian point inside: True
- Actual inside point inside: False
Spherical geometry results:
- Null meridian point inside: False
- Actual inside point inside: True
