Handling Sentinel-5P Footprints

Leider ist der Eintrag nur auf Amerikanisches Englisch verfügbar. Der Inhalt wird unten in einer verfügbaren Sprache angezeigt. Klicken Sie auf den Link, um die aktuelle Sprache zu ändern.

In the frame of our Earth Observation projects at Catalysts, we frequently stumble over small things, like efficiently handling satellite footprints. In this blog-post we want to share some of our experience when it comes to handling large area footprints that stretch over the dateline or the poles. It’s not rocket science per se, but maybe not that obvious either.

The Sentinel-5 Precursor mission flying the TROPOMI instrument was launched last October. It is observing our earth with a high spatial and spectral resolution to monitor various atmospheric aspects. You can read more about it here.

Since a few weeks ago, the data is being published on the Copernicus Open Access Hub for everyone to access. My colleague Alejandro and I had a quick look to see how well it can be used for our GRASP projects.

Total Column Carbon Monoxide Concentration,
S5P/TROPOMI Orbit 3954

There is not much going on in that orbit, there is a bit of CO emission near Jakutsk in Sibiria, Russia, likely due to vast wildfires in that region.

Output of S5P/TROPOMI Orbit 3954

Looking at the Footprint Boundaries

Just like most other earth observation datasets, Sentinel-5P offline products are disseminated in terms of orbits, typically 14 per day. When developing algorithms, and testing on a specific region, it is essential to locate the relevant data quickly. Going through all 14 orbits, loading all data and discarding what is not needed is tremendously slow.

A much quicker way is to only test the footprint boundaries for a match. This can typically be done with a handful of tests.

While the S5P products are pretty big in size they luckily have the footprint boundaries stored in an attribute:

footprint = [float(f) for f in getattr(ds['METADATA/EOP_METADATA/om:featureOfInterest/eop:multiExtentOf/gml:surfaceMembers/gml:exterior'], 'gml:posList').split(' ')]
flon, flat = footprint[1::2], footprint[0::2]

m = Basemap()
m.drawcoastlines()
m.scatter(flon, flat, marker=".", color="green", label="Outline")
S5P/TROPOMI Orbit 3954 Footprint

Given the wide swath of TROPOMI, 2600km, and the fact that the products are distributed as full (half) orbit files, the footprint boundaries will almost always cross the dateline. On the Earth’s surface this is fine, but on a 2D map the wrap-around is troublesome: … 175, 178, -179, -176, …

If libraries are not built to handle these cases, they’ll likely interpret the shape wrongly and mess up the covered area.

m = Basemap()
m.drawcoastlines()
fill(flon, flat, label="Fill")
S5P/TROPOMI Orbit 3954 Footprint

Testing Footprint Intersections

Shapely is usually a good tool to handle geometric operations, e.g. to test for polygon-polygon intersections. But when using it to handle the S5P orbits, it seems to have issues like the other 2D libraries:

from shapely import geometry
fp_poly = geometry.Polygon(list(zip(flon, flat)))
fp_poly
Output of Shapely

If we run our suite of test locations, we get this:

tests = [
    ("New Zealand", 170,-45, True),
    ("North Pole",    0, 90, True),
    ("Argentina",   -70,-45, False),
    ("Iceland",     -20, 66, True),
]

for (l,x,y,t) in tests:
    a = fp_poly.contains(geometry.Point(x,y))
    print (l, "is", "inside" if a else "outside", "(correct)" if a == t else "(*wrong*)")

New Zealand is outside (*wrong*)
North Pole is outside (*wrong*)
Argentina is inside (*wrong*)
Iceland is inside (correct)

All tests but one fail.

So, how do we fix this?

Well, I can think of three ways to handle it:

  1. Cut the one big footprint polygon in several smaller, so that each individual one has no dateline crossing anymore
  2. Define the footprint by a point and offsets from it
  3. Tinker with the coordinates by shifting and wrapping and using some fancy checks to test the intersection
  4. Reproject the coordinates into another space, where there is no discontinuity

For sake of simplicity, it would be best if ESA could change the products to include (1). But for the time being (4) seems like the best option.

Working with Spherical Geometries

One space in which no such discontinuity exists is the spherical domain. Luckily there is a ready-to-use Spherical Geometry package around, which can be installed via pip.

A SphericalPolygon can be constructed either via transforming lat/lon coordinates into cartesian, or via the from_lonlat constructor functions.

from spherical_geometry.polygon import SphericalPolygon

sp = SphericalPolygon.from_lonlat(lon=flon, lat=flat)
for (site,x,y,expect) in tests:
    actual = sp.contains_lonlat(x,y)
    print (site, "is", "inside" if a else "outside", "(correct)" if actual == expect else "(*wrong*)")

New Zealand is inside (*wrong*)
North Pole is inside (*wrong*)
Argentina is inside (*wrong*)
Iceland is inside (*wrong*)

Hmmm, this is not working at all. Did we get any of the parameters wrong?

A simple test

Let’s have a look at a simple, non-problematic test.

sp = SphericalPolygon.from_lonlat([-10, -10, 10, 10, -10], [-10, 10, 10, -10, -10])

sp.draw(m, latlon=True, color="red")
m.drawcoastlines()
m.drawmeridians(range(-180,180,45), labels=[0,0,0,1]); m.drawparallels(range(-90,90,45), labels=[1,0,0,0]);
Output
sp.contains_lonlat(0,0), sp.contains_lonlat(-45,-45)

(True, False)

All good. We are at least using it the right way and it is working in general.

Let’s look at an example near the dateline now.

sp = SphericalPolygon.from_lonlat([170, 170, -170, -170, 170], [-10, 10, 10, -10, -10])

sp.draw(m, latlon=True, color="red")
m.drawcoastlines()
m.drawmeridians(range(-180,180,45), labels=[0,0,0,1]); m.drawparallels(range(-90,90,45), labels=[1,0,0,0]);
Output
sp.contains_lonlat(175,0), sp.contains_lonlat(0,0)

(True, True)

The plot does not work – which we kind of expected – but the test doesn’t work either. That’s a bummer. After all we want to use it for exactly that purpose.

Now, if one looks closer at the documentation, there is actually an inside or center parameter to specify a point within the polygon:

sp = SphericalPolygon.from_lonlat([170, 170, -170, -170, 170], [-10, 10, 10, -10, -10], center=(175,0))
sp.contains_lonlat(175,0), sp.contains_lonlat(0,0)

(True, False)

Excellent, that seems to do the trick.

Back to the full footprint

Let’s see how things look with that little hint:

sp = SphericalPolygon.from_lonlat(lon=flon, lat=flat, center=(170,-45))
for (site,x,y,expect) in tests:
    actual = sp.contains_lonlat(x,y)
    print (site, "is", "inside" if a else "outside", "(correct)" if actual == expect else "(*wrong*)")

New Zealand is inside (correct)
North Pole is inside (correct)
Argentina is inside (correct)
Iceland is inside (correct)

Awesome, 4 out of 4!

Now, it is a bit odd to specify the center manually. For our processing we chose an arbitrary point of the (lon,lat) values from the band’s georeferencing data.

sp = SphericalPolygon.from_lonlat(lon=flon, lat=flat, center=(data_lon[10,10], data_lat[10,10]))
for (site,x,y,expect) in tests:
    actual = sp.contains_lonlat(x,y)
    print (site, "is", "inside" if a else "outside", "(correct)" if actual == expect else "(*wrong*)")

New Zealand is inside (correct)
North Pole is inside (correct)
Argentina is inside (correct)
Iceland is inside (correct)

You should, however, not choose any of the footprint boundaries itself:

sp = SphericalPolygon.from_lonlat(lon=flon, lat=flat, center=(flon[0], flat[0]))
for (site,x,y,expect) in tests:
    actual = sp.contains_lonlat(x,y)
    print (site, "is", "inside" if a else "outside", "(correct)" if actual == expect else "(*wrong*)")

New Zealand is inside (correct)
North Pole is inside (*wrong*)
Argentina is inside (correct)
Iceland is inside (*wrong*)

If you don’t want to touch the geocoding variables for the test, you could try to infer the point e.g. by finding the equatorial crossing. A product’s swath will always be smaller than the non-covered part of the equator.

Speed Comparison
%timeit [ sp.contains_lonlat(x,y) for x in range(-50,50) for y in range(-50,50) ]

4.16 s ± 156 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit [ fp_poly.contains(geometry.Point(x,y)) for x in range(-50,50) for y in range(-50,50) ]

2.15 s ± 44 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In terms of processing speed, it seems to be around half as fast as Shapely, with some 41ms per point-test.

A Proper Plot with Split Polygons

At the end of the day, it would be nice to see how the polygon looks filled on the globe.

To do this, we first setup two split polygons, one for the Western Hemisphere and one for the Eastern, and intersect each with the footprint Polygon:

e = 1e-6 # the spherical geometry library runs into assertions if using 90 or 180 exactly
east = SphericalPolygon.from_lonlat(lon=[0,0,180-e,180-e], lat=[90-e,-90+e,-90+e,90-e], center=(15,0))
west = SphericalPolygon.from_lonlat(lon=[-180+e,-180+e,0,0], lat=[90-e,-90+e,-90+e,90-e], center=(-15,0))
west_int, east_int = sp.intersection(west), sp.intersection(east)
split = list(west_int.to_lonlat()) + list(east_int.to_lonlat())

Note, I had to add/remove a tiny epsilon near the poles and the dateline to not run into some assertions.

Output

The orbit is now split in three polygons, as the intersection for the Western Hemisphere is not continuous. The SphericalPolygon.to_lonlat export yields longitude values in the [0..360] domain, and hence need to be converted.

Notes

There is a frequently a warning popping up about great_circle_arc.py:259: RuntimeWarning: invalid value encountered in intersects. With some footprints, this lead to a failure of the intersection testing. This seems to origin in the optimized C implementation, but is not around in the numpy version (which is actually the fall-back). I have not yet had the time to investigate this.

The actual notebook used to create all the plots can be found here.

Related Posts

No results found

1 Kommentar. Hinterlasse eine Antwort

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.

Bitte füllen Sie dieses Feld aus
Bitte füllen Sie dieses Feld aus
Bitte gib eine gültige E-Mail-Adresse ein.
Sie müssen den Bedingungen zustimmen, um fortzufahren

Menü