Crossover analysis - xopr Documentation

Radar sounder crossover analysis refers to finding crossing points between radar flight lines and analyzing differences between the radar data at the same point. Differences can arise from many factors. Often, crossing flight paths are at (nearly) perpendicular angles. Because the along-track and across-track beamwidth of most radar systems is very different, this different imaging geometry can lead to differences in the radar data. These differences are most pronouned over rough terrain, where off-nadir clutter may show up differently in coicident data collected along different angles. Temporal changes, changes in radar systems, and errors in picking the location of the surface or bed are other sources of differences in crossover analysis.

In this notebook, we demonstratate how to automatically find and analyze radar crossovers. To do this, we will

Find STAC items representing radar data within a geographic region.
Use the STAC item geometry to identify points where the radar flight paths cross (crossover points).
Selectively load layer information to make a map of the differences in WGS84 elevation between the bed picks at each crossover point. (This step is parallelized using Dask.)
Interactively load and plot radar data round selected crossover points to see what’s happening.

%load_ext autoreload
%autoreload 2

import numpy as np
import xarray as xr
import geoviews as gv
import geoviews.feature as gf
import cartopy.crs as ccrs
import cartopy.io.shapereader as shapereader
import matplotlib.pyplot as plt
import shapely.geometry
import scipy.constants
import geopandas as gpd
from tqdm import tqdm

import requests
import time

import xopr

import holoviews as hv
import hvplot.xarray
import hvplot.pandas
hvplot.extension('bokeh')

# Useful background features for maps
background_map = gf.ocean.opts(projection=ccrs.SouthPolarStereo(), scale='50m') * gf.coastline.opts(projection=ccrs.SouthPolarStereo(), scale='50m')

# Establish an OPR session
# You'll probably want to set a cache directory if you're running this locally to speed
# up subsequent requests. You can do other things like customize the STAC API endpoint,
# but you shouldn't need to do that for most use cases.
opr = xopr.OPRConnection(cache_dir="/tmp")

# Or you can open a connection without a cache directory (for example, if you're parallelizing
# this on a cloud cluster without persistent storage).
#opr = xopr.OPRConnection()

Step 1: Finding radar lines¶

For this notebook, we’ll experiment with finding data by geographic region. xOPR includes a helper module xopr.geometry with some useful utilities. You can call xopr.geometry.get_antarctic_regions to select one or more regions from the MEaSUREs Antarctic Boundaries dataset. Below, we plot all of the regions to give you some options for what to try. Mouse over each region to see its name.

regions_df = xopr.geometry.get_antarctic_regions(merge_regions=False).to_crs('EPSG:3031')
regions_df.hvplot(frame_width=600, aspect='equal', hover_cols=['NAME'], c='TYPE')

For our demonstration, we’ll use David Glacier, but feel free to trade this out for other locations and experiment.

region = xopr.geometry.get_antarctic_regions(name="David", merge_regions=True, simplify_tolerance=100)
region_projected = xopr.geometry.project_geojson(region, source_crs='EPSG:4326', target_crs="EPSG:3031")

region_hv = hv.Polygons([region_projected]).opts(
    color='green',
    line_color='black',
    fill_alpha=0.5)

(background_map * region_hv).opts(aspect='equal')

stac_items_df = opr.query_frames(geometry=region, date_range="2021-01-01T00:00:00Z/2025-06-01T00:00:00Z") # Can also add a date range: date_range="2020-01-01T00:00:00Z/2024-01-01T00:00:00Z"
stac_items_df = stac_items_df.to_crs('EPSG:3031')

print(f"Found {len(stac_items_df)} frames across {stac_items_df['collection'].nunique()} collections:")
stac_items_df.groupby('collection').size()

Found 34 frames across 2 collections:

collection
2022_Antarctica_BaslerMKB    17
2023_Antarctica_BaslerMKB    17
dtype: int64

flight_lines = stac_items_df.hvplot(by='collection')
(background_map * region_hv * flight_lines).opts(frame_width=500, aspect='equal', active_tools=['pan', 'wheel_zoom'])

Step 2: Identify crossover points¶

xOPR includes a helper function to identify crossover points. Feel free to poke into the code and take a look. It takes advantage of GeoPandas’s spatial join (sjoin) function, so it’s actually quite concise.

intersections = xopr.find_intersections(stac_items_df, calculate_crossing_angles=True)

intersections = intersections[intersections['crossing_angle'] > 2] # Filter out nearly coincident crossings

print(f"Found {len(intersections)} crossover points between flight lines.")
(background_map * region_hv * flight_lines * intersections.hvplot(label='Intersection Points', c='crossing_angle')).opts(frame_width=500, aspect='equal', active_tools=['pan', 'wheel_zoom'])

Zoom in and check out the plot above. Every identified crossover point has a blue dot over it. The blue dots are shaded by the crossing angle. You’ll see that nearly coincident lines will have a large number of low crossing angles wheras most other crossovers are at angles closer to 90 degrees.

The cell below shows you what the result of xopr.find_intersections() looks like. It returns a GeoDataFrame where every column from the intersecting frames is preserved, with suffixes _1 and _2 to distinguish them.

intersections.head()

Step 3: Load layer data and find the different in bed elevation¶

Depending on how many intersection points you have, loading all of the layer data associated with each crossover may take quite some time. Luckily, this process can be parallelized fairly efficiently, so we’ll use Dask to make this run a bit faster. Take a look at the “Search and Scaling” notebook for some tips on working with Dask.

This is a good time to check that you have a reasonable number of intersection points if you’ve modified anything. We recommend starting with about 100 or fewer intersection points so that this won’t take more than about 5 minutes to run.

import dask
import dask.delayed as delayed
from dask import compute
from dask.distributed import LocalCluster

client = LocalCluster().get_client()
client

By default, Dask will start a server at http://127.0.0.1:8787/status. If it can’t use that port for some reason, the cell above should print out an alternative port. Visiting this status page is optional, but it will tell you about the progress of your job.

@delayed
def safe_get_layers_db(stac_item, opr=xopr.opr_access.OPRConnection()):
    try:
        retries = 1
        backoff_time = 5
        backoff_jitter = 30
        while retries > 0:
            try:
                return opr.get_layers_db(stac_item)
            except requests.exceptions.RequestException as e:
                sleep_time = backoff_time + np.random.uniform(0, backoff_jitter)
                print(f"Request error fetching layers for {stac_item['id']}: {e}. Retrying in {sleep_time:.1f} seconds...")
                time.sleep(sleep_time)
                retries -= 1
                backoff_time *= 2  # Exponential backoff
    except Exception as e:
        print(f"Error fetching layers for {stac_item['id']}: {e}")
        return None

def get_basal_layer_wgs84(stac_item, preloaded_layer=None, opr=xopr.opr_access.OPRConnection()):
    if (preloaded_layer is None) or len(preloaded_layer) < 2:
        layers = opr.get_layers_files(stac_item)
    else:
        layers = preloaded_layer
    
    basal_layer = layers[2]
    surface_layer = layers[1]

    surface_wgs84 = layers[1]['elev'] - (layers[1]['twtt'] * (scipy.constants.c / 2))
    delta_twtt = basal_layer['twtt'] - surface_layer['twtt']
    basal_wgs84 = surface_wgs84 - (delta_twtt * ((scipy.constants.c / np.sqrt(3.15)) / 2))

    basal_layer['wgs84'] = basal_wgs84
    return basal_layer

@delayed
def compute_crossover_error_impl(stac_item_1, stac_item_2, intersection_geometry, layer_1, layer_2):
    """Implementation that receives actual layer values"""
    try:
        bed_1 = get_basal_layer_wgs84(stac_item_1, preloaded_layer=layer_1).rename({'lat': 'Latitude', 'lon': 'Longitude'})
        bed_2 = get_basal_layer_wgs84(stac_item_2, preloaded_layer=layer_2).rename({'lat': 'Latitude', 'lon': 'Longitude'})

        bed_1 = xopr.geometry.project_dataset(bed_1, "EPSG:3031")
        bed_2 = xopr.geometry.project_dataset(bed_2, "EPSG:3031")

        x, y = intersection_geometry.coords[0]

        dist_1 = np.sqrt((bed_1['x'] - x)**2 + (bed_1['y'] - y)**2)
        dist_2 = np.sqrt((bed_2['x'] - x)**2 + (bed_2['y'] - y)**2)

        min_idx_1 = dist_1.argmin().item()
        min_idx_2 = dist_2.argmin().item()

        dist_between_pts = np.sqrt((bed_1['x'][min_idx_1] - bed_2['x'][min_idx_2])**2 + (bed_1['y'][min_idx_1] - bed_2['y'][min_idx_2])**2)

        elev_1 = bed_1['wgs84'][min_idx_1].item()
        elev_2 = bed_2['wgs84'][min_idx_2].item()

        return elev_1, elev_2, dist_between_pts
    except Exception as e:
        print(f"Error in compute_crossover_error: {e}")
        return None, None, None  # Return sentinel values on error

def compute_crossover_error(stac_item_1, stac_item_2, intersection_geometry, preloaded_layer_1=None, preloaded_layer_2=None):
    """Wrapper that handles delayed objects properly"""
    # These will be delayed objects or None
    return compute_crossover_error_impl(stac_item_1, stac_item_2, intersection_geometry, preloaded_layer_1, preloaded_layer_2)

future_db_layers = {}
future_results = {}
for idx, row in intersections.iterrows():
    stac_item_1 = stac_items_df.loc[row['id_1']].to_dict()
    stac_item_2 = stac_items_df.loc[row['id_2']].to_dict()
    stac_item_1['id'] = row['id_1']
    stac_item_2['id'] = row['id_2']

    # Fetch the layers from the database if available
    # Updated to use opr:segment instead of opr:flight (schema change)
    db_key_1 = (row['collection_1'], row['opr:date_1'], row['opr:segment_1'])
    db_key_2 = (row['collection_2'], row['opr:date_2'], row['opr:segment_2'])
    if db_key_1 not in future_db_layers:
        future_db_layers[db_key_1] = safe_get_layers_db(stac_item_1)
    if db_key_2 not in future_db_layers:
        future_db_layers[db_key_2] = safe_get_layers_db(stac_item_2)

    # Create delayed task but DON'T compute yet
    r = compute_crossover_error(
        stac_item_1, stac_item_2, row.intersection_geometry,
        preloaded_layer_1=future_db_layers.get(db_key_1),
        preloaded_layer_2=future_db_layers.get(db_key_2)
    )
    future_results[idx] = r  # Store the delayed object, not the computed result

# Compute all results in parallel
results = compute(*future_results.values())

# Process the results
for idx, result in zip(future_results.keys(), results):
    try:
        elev_1, elev_2, dist_between_pts = result
        if elev_1 is not None and elev_2 is not None:  # Skip failed computations
            intersections.at[idx, 'wgs84_1'] = elev_1
            intersections.at[idx, 'wgs84_2'] = elev_2
            intersections.at[idx, 'layer_pt_distance'] = dist_between_pts
        else:
            print(f"Skipping intersection {idx} due to computation error")
    except Exception as e:
        print(f"Error processing intersection {idx}: {repr(e)}")

We now have all of our crossovers. Don’t worry about a few errors in the cell above. This is normal as it tries to load layer data from both possible sources (files and database).

intersections['elev_diff'] = np.abs(intersections['wgs84_1'] - intersections['wgs84_2'])
# Set elev_diff to NaN where layer_pt_distance is large
intersections.loc[intersections['layer_pt_distance'] > 100, 'elev_diff'] = np.nan

intersections_success = intersections.dropna().reset_index(drop=True)
intersections_success['intersection_geometry_x'] = intersections_success['intersection_geometry'].apply(lambda geom: geom.x)
intersections_success = intersections_success.sort_values(by='intersection_geometry_x', ascending=False).reset_index(drop=True)
intersections_success['idx'] = intersections_success.index
hover_tooltips = [
    ("Index", "@idx"),
    ("Collection 1", "@collection_1"),
    ("Collection 2", "@collection_2"),
    ("Difference", "@elev_diff{0.00} m"),
]

vlim = intersections_success['elev_diff'].abs().quantile(0.9)

hv_int = intersections_success.hvplot(color='elev_diff', hover_cols=['idx', 'collection_1', 'collection_2', 'elev_diff'], hover_tooltips=hover_tooltips, clim=(0, vlim))
#hv_int = hv_int.opts(scalebar=True) # Can enable if you want - requires hvplot >= 0.12.0
(background_map * region_hv * flight_lines * hv_int).opts(frame_width=600, aspect='equal', active_tools=['pan', 'wheel_zoom'])

This map shows the differences in picked bed elevation at every crossing point where layer data was available. Explore around and see what you notice!

Step 4: Investigate individual crossovers by looking at the radar data¶

If you hover over any crossover point in the map above, you’ll get the index associated with each crossover. Enter one of those indices below to load and plot the corresponding radar data to see what’s happening.

selected_idx = 4

# Get intersection details
intersect = intersections_success.loc[selected_idx]
stac_1 = stac_items_df.loc[intersect['id_1']].to_dict()
stac_2 = stac_items_df.loc[intersect['id_2']].to_dict()

# Load frames
frame_1 = opr.load_frame(stac_1)
frame_1 = xopr.radar_util.add_along_track(frame_1)
frame_1 = xopr.radar_util.interpolate_to_vertical_grid(frame_1, vertical_coordinate='wgs84')
frame_2 = opr.load_frame(stac_2)
frame_2 = xopr.radar_util.add_along_track(frame_2)
frame_2 = xopr.radar_util.interpolate_to_vertical_grid(frame_2, vertical_coordinate='wgs84')

# Project to EPSG:3031 and find closest points to intersection
x_int, y_int = intersect.intersection_geometry.coords[0]
frame_1_proj = xopr.geometry.project_dataset(frame_1, "EPSG:3031")
frame_2_proj = xopr.geometry.project_dataset(frame_2, "EPSG:3031")

# Find indices closest to intersection
dist_1 = np.sqrt((frame_1_proj['x'] - x_int)**2 + (frame_1_proj['y'] - y_int)**2)
dist_2 = np.sqrt((frame_2_proj['x'] - x_int)**2 + (frame_2_proj['y'] - y_int)**2)
idx_1 = dist_1.argmin().item()
idx_2 = dist_2.argmin().item()

print(f"Frame 1: {intersect['id_1']} from {intersect['collection_1']}")
print(f"Frame 2: {intersect['id_2']} from {intersect['collection_2']}")
print(f"Intersection at index {idx_1} (frame 1) and {idx_2} (frame 2)")
print(f"Bed elevation difference: {intersect['elev_diff']:.2f} m")

Frame 1: Data_20221212_01_013 from 2022_Antarctica_BaslerMKB
Frame 2: Data_20231212_02_012 from 2023_Antarctica_BaslerMKB
Intersection at index 264 (frame 1) and 1275 (frame 2)
Bed elevation difference: 481.48 m

def add_layers_to_frame(frame):
    speed_in_ice = scipy.constants.c / np.sqrt(3.15)  # Approximate speed of light in ice (m/s)

    layers = opr.get_layers(frame)
    if len(layers) == 0:
        print(f"No layers found for granule {frame.attrs['granule']}")
        return frame
    
    has_surface = False

    for layer_idx in layers:
        frame[f'layer_{layer_idx}_twtt'] = layers[layer_idx]['twtt'].interp(coords={'slow_time': frame['slow_time']})

        if layer_idx == 1:
            has_surface = True
            frame[f'layer_{layer_idx}_range'] = frame['layer_1_twtt'] * (scipy.constants.c / 2)
            frame[f'layer_{layer_idx}_elevation'] = frame['Elevation'] - frame[f'layer_1_range']
        elif has_surface:
            twtt_from_surface = frame[f'layer_{layer_idx}_twtt'] - frame['layer_1_twtt']
            frame[f'layer_{layer_idx}_range'] = (twtt_from_surface * (speed_in_ice / 2)) + frame['layer_1_range']
            frame[f'layer_{layer_idx}_elevation'] = frame['Elevation'] - frame[f'layer_{layer_idx}_range']

    return frame

# Load layers for both frames
frame_1 = add_layers_to_frame(frame_1)
frame_2 = add_layers_to_frame(frame_2)

clb_min_pct, clb_max_pct = 30, 97

# Plot radargrams in elevation coordinates with layers
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 8))

# Frame 1 radargram in elevation
pwr_1_elev = 10*np.log10(np.abs(frame_1.Data))
vmax_1 = np.percentile(pwr_1_elev, clb_max_pct)
vmin_1 = np.percentile(pwr_1_elev, clb_min_pct)
pwr_1_elev.plot.imshow(x='along_track', y='wgs84', cmap='gray', ax=ax1, vmin=vmin_1, vmax=vmax_1)
ax1.axvline(frame_1.along_track[idx_1].values, color='red', linestyle='--', linewidth=2, label='Crossover')

# Plot layers using elevation data
if 'layer_1_elevation' in frame_1:
    frame_1['layer_1_elevation'].plot(ax=ax1, x='along_track', linestyle=':', color='tab:blue', linewidth=2, label='Surface')
if 'layer_2_elevation' in frame_1:
    frame_1['layer_2_elevation'].plot(ax=ax1, x='along_track', linestyle=':', color='tab:orange', linewidth=2, label='Bed')

ax1.set_title(f"{intersect['collection_1']} - {intersect['id_1']} (Elevation view)")
ax1.set_ylabel('Elevation (m)')
ax1.legend()

# Frame 2 radargram in elevation
pwr_2_elev = 10*np.log10(np.abs(frame_2.Data))
vmax_2 = np.percentile(pwr_2_elev, clb_max_pct)
vmin_2 = np.percentile(pwr_2_elev, clb_min_pct)
pwr_2_elev.plot.imshow(x='along_track', y='wgs84', cmap='gray', ax=ax2, vmin=vmin_2, vmax=vmax_2)
ax2.axvline(frame_2.along_track[idx_2].values, color='red', linestyle='--', linewidth=2, label='Crossover')

# Plot layers using elevation data
if 'layer_1_elevation' in frame_2:
    frame_2['layer_1_elevation'].plot(ax=ax2, x='along_track', linestyle=':', color='tab:blue', linewidth=2, label='Surface')
if 'layer_2_elevation' in frame_2:
    frame_2['layer_2_elevation'].plot(ax=ax2, x='along_track', linestyle=':', color='tab:orange', linewidth=2, label='Bed')

ax2.set_title(f"{intersect['collection_2']} - {intersect['id_2']} (Elevation view)")
ax2.set_xlabel('Along track distance (m)')
ax2.set_ylabel('Elevation (m)')
ax2.legend()

plt.suptitle(f"Radargrams in elevation coordinates - Bed elev diff: {intersect['elev_diff']:.2f} m", fontsize=14)
plt.tight_layout()
plt.show()

# Zoomed elevation plots around crossover
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 8))

window_size = 300  # Number of traces on each side of intersection
elev_window = np.maximum(intersect['elev_diff']*2.5, 100)  # Elevation window in meters around bed

# Frame 1 zoomed
idx_start_1 = max(0, idx_1 - window_size)
idx_end_1 = min(len(frame_1.slow_time), idx_1 + window_size)

pwr_1_zoom = frame_1.Data.isel(slow_time=slice(idx_start_1, idx_end_1))
pwr_1_zoom_db = 10*np.log10(np.abs(pwr_1_zoom))
pwr_1_zoom_db.plot.imshow(x='along_track', y='wgs84', cmap='gray', ax=ax1, vmin=vmin_1, vmax=vmax_1)
ax1.axvline(frame_1.along_track[idx_1].values, color='red', linestyle='--', linewidth=2, label='Crossover')

# Plot layers
if 'layer_1_elevation' in frame_1:
    layer_1_zoom = frame_1['layer_1_elevation'].isel(slow_time=slice(idx_start_1, idx_end_1))
    layer_1_zoom.plot.scatter(ax=ax1, x='along_track', linestyle=':', color='tab:blue', linewidth=2, label='Surface', s=4)
if 'layer_2_elevation' in frame_1:
    layer_2_zoom = frame_1['layer_2_elevation'].isel(slow_time=slice(idx_start_1, idx_end_1))
    layer_2_zoom.plot.scatter(ax=ax1, x='along_track', linestyle=':', color='tab:orange', linewidth=2, label='Bed', s=4)

    bed_elev_1 = frame_1['layer_2_elevation'].isel(slow_time=idx_1).values
    elev_min_1 = bed_elev_1 - elev_window/2
    elev_max_1 = bed_elev_1 + elev_window/2
    ax1.set_ylim(elev_min_1, elev_max_1)
ax1.set_title(f"{intersect['collection_1']} - Zoomed (Bed: {intersect['wgs84_1']:.1f} m)")
ax1.set_ylabel('Elevation (m)')
ax1.legend()

# Frame 2 zoomed
idx_start_2 = max(0, idx_2 - window_size)
idx_end_2 = min(len(frame_2.slow_time), idx_2 + window_size)

pwr_2_zoom = frame_2.Data.isel(slow_time=slice(idx_start_2, idx_end_2))
pwr_2_zoom_db = 10*np.log10(np.abs(pwr_2_zoom))
pwr_2_zoom_db.plot.imshow(x='along_track', y='wgs84', cmap='gray', ax=ax2, vmin=vmin_2, vmax=vmax_2)
ax2.axvline(frame_2.along_track[idx_2].values, color='red', linestyle='--', linewidth=2, label='Crossover')

# Plot layers
if 'layer_1_elevation' in frame_2:
    layer_1_zoom = frame_2['layer_1_elevation'].isel(slow_time=slice(idx_start_2, idx_end_2))
    layer_1_zoom.plot.scatter(ax=ax2, x='along_track', linestyle=':', color='tab:blue', linewidth=2, label='Surface', s=4)
if 'layer_2_elevation' in frame_2:
    layer_2_zoom = frame_2['layer_2_elevation'].isel(slow_time=slice(idx_start_2, idx_end_2))
    layer_2_zoom.plot.scatter(ax=ax2, x='along_track', linestyle=':', color='tab:orange', linewidth=2, label='Bed', s=4)

    bed_elev_2 = frame_2['layer_2_elevation'].isel(slow_time=idx_2).values
    elev_min_2 = bed_elev_2 - elev_window/2
    elev_max_2 = bed_elev_2 + elev_window/2
    
    ax2.set_ylim(elev_min_2, elev_max_2)
ax2.set_title(f"{intersect['collection_2']} - Zoomed (Bed: {intersect['wgs84_2']:.1f} m)")
ax2.set_xlabel('Along track distance (m)')
ax2.set_ylabel('Elevation (m)')
ax2.legend()

plt.suptitle(f"Elevation crossover comparison - Bed difference: {intersect['elev_diff']:.2f} m", fontsize=14)
plt.tight_layout()
plt.show()

Notebooks

Reading parameters