Resample altimetry swath data to a grid

[nilsleh]

I recently started working with swath altimetry data, namely Level 2 SWOT data which comes in the .nc files of the following format:

However, for actual analysis I would like to regrid/resample that collection of observations to a regular grid of (time, lat, lon) such that given spatial extents for specific timestamps can easily be compared with other available data products of form like:

I am therefore wondering, whether pyresample would be capable of doing such a task? I have uploaded a few files of the observational data to google drive if that is of any help.

As a reference, I also tried inquiring on the Earthdata Forum from which the L2 SWOT product is downloaded.

Thanks in advance for any help and suggestions!

[djhoese]

Yes, pyresample should be able to handle that. You would have to resample each time-step separately I think, but it may depend on the interface you use. This might almost be easiest with Satpy which although it can't read SWOT data directly with a satpy "reader", you could make a Satpy Scene object, add a DataArray for each variable you'll resample (with an added "area" attr to represent the SwathDefinition of the data), then call new_scn = scn.resample(output_lonlat_grid). Then new_scn will be a container of the resampled DataArrays.

Or there are the existing pyresample interfaces, but it's been a while since I've used them myself and the other pyresample maintainers are currently on holiday so sorry I can't give a more detailed answer.

[nilsleh]

@djhoese sorry to ping you, but I was wondering whether you or your colleagues are available to take a look at this again. Thanks a lot in advance.

[djhoese]

Sorry. Life caused me to lose track this browser tab. A couple things to fix:

1. Don't do .values on the lon/lats. SwathDefinition for Satpy should have DataArrays.
2. Your grid_lons/lats could be a SwathDefinition but they'd have to be 2D arrays. But because they are a grid you could make it an AreaDefinition with a +proj=longlat projection. You'd just have to do the math to find the extents (the outer edges of the edge pixels, not the centers of those pixels). You'd then pass this AreaDefinition to resample.

I'm not sure why the input is the same as the output. Satpy should have screamed at you for passing an xr.Dataset as the target grid of resampling (which is not supported). It needs to be an AreaDefinition or in some cases a SwathDefinition might work.

[nilsleh]

No worries, I must admit I am still stuck due to my lack of knowledge of projections, satpy, pyresample, proj4 etc.

from satpy import Scene
from pyresample.geometry import AreaDefinition
from pyresample.utils import check_and_wrap

# single along track with dimensions num_pixels, num_line
swot_data = xr.open_dataset("/mnt/SSD2/nils/datasets/swot/downloaded_data/SWOT_L2_LR_SSH_Basic_010_210_20240201T113157_20240201T122325_PIC0_01.nc")

projection = "EPSG:4326"

area_extent = (-90, 90, -180, 180)

area_def = AreaDefinition(
    area_id="swot_area",
    description="SWOT", 
    proj_id="swot_proj",
    projection=projection,
    width=1000,
    height=1000,
    area_extent=area_extent
)

swot_scene = Scene()
swot_scene["ssh_karin"] = swot_data["ssh_karin"]

new_scene = swot_scene.resample(area_def)
new_scene.values()

Essentially, what I am looking for is that from all the individual files that contain a single along track to generate a "composite" similar to this:

I am sorry for the lack or precision with these questions, so very appreciative for any further remarks.

[mraspaud]

If the aim is to generate this plot, I would recommend looking at eg. cartopy (this might work https://stackoverflow.com/a/61053504).
Or are you using pyresample for performance or running non-interactively?

[nilsleh]

The ultimate goal would be to have such along track data in a xarray dataset that has dimensions time / lat / lon where the time index would be one composite over the revisiting period essentially. I would like to use the gridded L2 SWOT data in machine learning applications, where I need to extract patches of a given spatial lat/lon extent across different data sources. Thus having all data sources defined on a grid would be very convenient to make data indexing more easy for the downstream application.