Regrid SWOT data to regular grid

[nilsleh]

Dear SWOT-community,

I downloaded SWOT L3 data from Aviso for a research project. I am new to SWOT data and these sort of data sources in general, so apologies in advance.

swot = xr.open_dataset("SWOT_L3_LR_SSH_Basic_001_584_20230811T012726_20230811T021852_v1.0.nc")

However, I would like to match the SWOT data with other data products that come in a regular grid like (time, lat, lon) format, like so:

I have done some searching with respect to possible tools like:

• xESMF
• pygeo-pyinterp
• HarmonyPy
• and likely others

But was wondering if you have a recommended workflow or possibly even a suggestion on how to do this? The SWOT data wouldn't need to be in the same resolution as the other data sources, but should be regridded to be in lat/lon format. Thank you for any hints, suggestions, and possible solutions. Cheers.

[iuryt]

This is an old post so you probably have found a workaround for that. An easy way to do that is to interpolate your gridded data to SWOT, since you have the coordinate longitude, latitude defined as dimensions.

gridded.interp(longitude=swot.longitude, latitude=swot.latitude)

I hope that helps.

[nilsleh]

@iuryt Thanks for your reply, actually I have still not found a good solution for this. I have tried something like:

import xarray as xr
import xesmf as xe
import numpy as np

def regrid_dataset_to_target_grid(source_dataset, target_latitudes, target_longitudes, method='bilinear'):
    """
    Regrids a source dataset to a target grid using the specified method.

    Parameters:
    - source_dataset (xarray.Dataset): The source dataset to be regridded.
    - target_latitudes (numpy.ndarray): Array of target latitudes for the new grid.
    - target_longitudes (numpy.ndarray): Array of target longitudes for the new grid.
    - method (str): The regridding method to use. Defaults to 'bilinear'.

    Returns:
    - xarray.Dataset: The regridded dataset with the target grid.
    """

    # Define the target grid
    target_grid = xr.Dataset({
        'lat': (['lat'], target_latitudes),
        'lon': (['lon'], target_longitudes)
    })

    # Create the regridder object
    regridder = xe.Regridder(source_dataset, target_grid, method)

    # Fix for numpy 2 while waiting for xESMF 0.8.7
    regridder.shape_in = tuple(map(int, regridder.shape_in))
    regridder.shape_out = tuple(map(int, regridder.shape_out))

    # Perform the regridding
    regridded_ds = regridder(source_dataset)

    # Copy the 'time' variable from the source dataset to the regridded dataset
    if 'time' in source_dataset:
        regridded_ds['time'] = source_dataset['time']

    return regridded_ds

swot = xr.open_dataset("SWOT_L2_LR_SSH_Basic_001_149_20230726T122755_20230726T131924_PGC0_01.nc")


lats = swot.latitude.values
lons = swot.longitude.values

grid_lats = np.linspace(lats.min(), lats.max(), 100)
grid_lons = np.linspace(lons.min(), lons.max(), 100)

regridded_dataset = regrid_dataset_to_target_grid(swot, grid_lats, grid_lons)

With Level 2 or Level 3 data downloaded from Earthdata or Aviso.

However, I think what I am looking for a more nuanced solution. Basically for a smaller lat/lon extent I am looking to extract a patch that respects these bounds and has dimension lat/lon and time.

So with the same data input i have also tried the following:

def extract_swot_patch(ds, lat_bounds: tuple[float], lon_bounds: tuple[float], resolution=0.1):
    """Extract a gridded patch of SWOT SSH data within given boundaries.
    
    Args:
        ds: xarray Dataset with SWOT data
        lat_bounds: tuple of (lat_min, lat_max), overrides center_lat/patch_size
        lon_bounds: tuple of (lon_min, lon_max), overrides center_lon/patch_size
        center_lat/lon: patch center coordinates (used if bounds not provided)
        patch_size: size in degrees (used if bounds not provided)
        resolution: grid spacing in degrees
    """
    # Determine boundaries
    lat_min, lat_max = lat_bounds
    lon_min, lon_max = lon_bounds
        
    # Create target grid
    # n_points = int((lat_max - lat_min)/resolution)
    n_points=512
    lat_grid = np.linspace(lat_min, lat_max, n_points)
    lon_grid = np.linspace(lon_min, lon_max, n_points)
    
    # Get SWOT data as 1D arrays 
    lats = ds.latitude.values.ravel()
    lons = ds.longitude.values.ravel()
    ssh = ds.ssha_karin.values.ravel()
    
    # Filter points within region
    mask = ((lats >= lat_min) & (lats <= lat_max) &
            (lons >= lon_min) & (lons <= lon_max))
    
    if not np.any(mask):
        return None
    L_x = 960e3  # size of domain
    L_y = 960e3  # size of domain
    # Grid the data using binned statistics (similar to bin_ssh)
    ssh_grid, _, _, _ = stats.binned_statistic_2d(
        lats[mask], lons[mask], ssh[mask],
        statistic='mean',
        bins=[n_points, n_points], 
        range=[[lat_min, lat_max], [lon_min, lon_max]]
        # range=[[-L_x / 2, L_x / 2], [-L_y / 2, L_y / 2]]
    )
    
    # Create output dataset
    patch_ds = xr.Dataset(
        data_vars={
            'ssha_karin': (['latitude', 'longitude'], ssh_grid),
        },
        coords={
            'latitude': lat_grid,
            'longitude': lon_grid
        }
    )
    
    # Copy attributes
    patch_ds.attrs = ds.attrs.copy()
    patch_ds.attrs.update({
        'lat_bounds': (lat_min, lat_max),
        'lon_bounds': (lon_min, lon_max),
        'resolution': resolution
    })
    
    return patch_ds

which is based on a script I found here.

So to summarize I am looking for a way to take Level 2 or 3 data that comes in the above outlined dimension of num_lines, num_pixels and for a particular lat/lon time of interest obtain a patch that contains interpolated/regridded SWOT data, such that I could combine it with other lat/lon time indexed data such as SST.