Search SWOT
This notebook demonstrates how to find half orbits intersecting a given geographical area.
Tutorial Objectives
Retrieve the half-orbit numbers intersecting a given geographical area.
Note
Required environment to run this notebook:
geopandas+shapelymatplotlib+cartopyfcollections: available here.
Import + code
import warnings
warnings.filterwarnings('ignore')
import geopandas as gpd
from shapely import geometry
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.patches as patches
def get_half_orbits_intersect(bbox):
"""Get half orbits that intersect a bounding box.
Parameters
----------
bbox:
the bounding box
Returns
-------
gpd.GeoDataFrame:
A Geopandas dataframe containing intersecting half orbits numbers and geometries
"""
swath_geometries = gpd.read_file(GEOMETRIES_FILE)
bbox_polygon = geometry.box(*bbox)
def _filter_intersect(row, polygon):
half_orbit_polygon = row.geometry
return polygon.intersects(half_orbit_polygon)
select = swath_geometries.apply(_filter_intersect, polygon=bbox_polygon, axis=1)
return swath_geometries[select]
def plot_geometries(geometries, lon_range, lat_range, title, plot_extent=None):
fig, ax = plt.subplots(ncols=1, figsize=(12, 8), subplot_kw={'projection': ccrs.PlateCarree()})
gpd.GeoSeries(geometries.geometry).plot(ax=ax,transform=ccrs.PlateCarree(),alpha=1)
square = patches.Rectangle((lon_range[0], lat_range[0]), lon_range[1]-lon_range[0], lat_range[1]-lat_range[0], edgecolor='orange', facecolor='none', transform=ccrs.PlateCarree())
ax.add_patch(square)
ax.set_title(title)
ax.coastlines()
if plot_extent:
ax.set_extent(plot_extent, crs=ccrs.PlateCarree())
Parameters
Define a geographical area
# Mediterranean sea
lon_range = 6, 11
lat_range = 38, 43
bbox = [lon_range[0], lat_range[0], lon_range[1], lat_range[1]]
plot_extent = [lon_range[0], lon_range[1], lat_range[0], lat_range[1]]
# science, calval
phase = "science"
Getting Swot geometries file using fcollections.sad
import fcollections.sad
fcollections.sad.summary()
┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ Type ┃ Available ┃ Keys ┃ Lookup Folders ┃ ┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ gshhg │ 0/15 │ GSHHS_c,GSHHS_f,GSHHS_h,GSHHS_i,GSHHS_l,border_c,bo │ /home/atonneau/.config/sad │ │ │ │ rder_f,border_h,border_i,border_l,river_c,river_f,r │ │ │ │ │ iver_h,river_i,river_l │ │ │ karin_footprints │ 2/2 │ calval,science │ /home/atonneau/.config/sad │ └──────────────────┴───────────┴─────────────────────────────────────────────────────┴────────────────────────────┘
aux = fcollections.sad.KarinFootprints()
aux.keys
{'calval', 'science'}
GEOMETRIES_FILE = aux[phase]
Search for matching half orbits
swath_geoms = get_half_orbits_intersect(bbox)
swath_geoms
| pass_number | geometry | |
|---|---|---|
| 0 | 1 | POLYGON ((-83.3713 -77.05377, -82.41365 -77.05... |
| 28 | 29 | POLYGON ((-85.83706 -77.05377, -84.8794 -77.05... |
| 235 | 236 | POLYGON ((-65.48095 78.27207, -64.06949 78.268... |
| 250 | 251 | POLYGON ((-79.67267 -77.05377, -78.71502 -77.0... |
| 263 | 264 | POLYGON ((-67.94671 78.27207, -67.15265 78.270... |
| 278 | 279 | POLYGON ((-82.13843 -77.05377, -81.18077 -77.0... |
| 291 | 292 | POLYGON ((-70.41246 78.27207, -69.001 78.26856... |
| 306 | 307 | POLYGON ((-84.60418 -77.05377, -83.64652 -77.0... |
| 513 | 514 | POLYGON ((-64.24808 78.27207, -63.45402 78.270... |
| 541 | 542 | POLYGON ((-66.71383 78.27207, -65.30237 78.268... |
| 556 | 557 | POLYGON ((-80.90555 -77.05377, -79.94789 -77.0... |
| 569 | 570 | POLYGON ((-69.17958 78.27207, -67.76812 78.268... |
half_orbits = list(swath_geoms['pass_number'])
half_orbits
[1, 29, 236, 251, 264, 279, 292, 307, 514, 542, 557, 570]
Plot geometries
plot_geometries(swath_geoms, lon_range, lat_range, f"Swot {phase} Geometries")
half_orbit_num = 264
plot_geometries(swath_geoms[swath_geoms['pass_number']==half_orbit_num], lon_range, lat_range, f"Swot {phase} Geometry, half orbit={half_orbit_num}", plot_extent)
