Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Export your netCDF file to SEAScope

import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
import os
import datetime
import pyproj
import numpy as np
from netCDF4 import Dataset, num2date
import pyproj

from IPython.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

Export your netCDF file to SEAScope

This notebook will give you an example on how to export your own data to SEAScope. In this example, we will open CFOSAT SWIM L2S product and compute the proper ground control point necessary for SEAScope. We will create the Collection and variable before sending data to SEAScope.

Before running this notebook:

  • Open SEAScope to start the communication

Read SWIM L2S spectral data

In the following cell, we will

  • Uncompress it

  • read a netcdf file and retrieve the variable of interest to be send to SEAScope

# Get the data from the ftp (not working on Windows)
file = 'CFO_OP06_SWI_L2S____F_20240217T081133_20240217T094440_1.0.2'
print(f'Check if {file} exists')
if not os.path.exists(file):
    print(f'Check if zip file exists: {file} does not exists')
    if not os.path.exists(f'{file}.zip'):
        print("The zip file will be downloaded, it will take ~5min")
        !wget https://ftp.odl.bzh/odl/events/space_week_2024/CFO_OP06_SWI_L2S____F_20240217T081133_20240217T094440_1.0.2.zip
# Uncompress the directory
if not os.path.exists(file):
    print(f'unzip {file}.zip')    
    import zipfile
    with zipfile.ZipFile('./CFO_OP06_SWI_L2S____F_20240217T081133_20240217T094440_1.0.2.zip', 'r') as zip_ref:
        zip_ref.extractall('./')

In the uncompressed directory, you should have one file per beam angle There are four beams:

  • ‘02’

  • ‘04’

  • ‘06’

  • ‘08’

  • ‘10’

In the following cell choose the beam you want to display in SEAScope in the variable `beam`
# Choose the beam you want to process
beam = '04'
# Read one netcdf file

prod = '20240217T081133_20240217T094440_1.0.2'

l2s_path = 'CFO_OP06_SWI_L2S____F_'+prod+'/CFO_OP06_SWI_L2S'+beam+'__F_'+prod+'.nc'
clat = 30
clon = -22 # Fev 17 2024 storm atlantic



dset = Dataset(l2s_path, 'r')

fluct_specs = dset.variables['fluctuation_spectra'][:]
mod_specs = dset.variables['modulation_spectra'][:]
wave_spec = dset.variables['wave_spectra'][:]
seg_sea_ice_concentration = dset.variables['seg_sea_ice_concentration'][:]
klin = dset.variables['klin'][:]
near_lat = dset.variables['near_lat'][:]
near_lon = dset.variables['near_lon'][:]
far_lat = dset.variables['far_lat'][:]
far_lon = dset.variables['far_lon'][:]
lat = dset.variables['lat'][:]
lon = dset.variables['lon'][:]
time = dset.variables['time'][:]

dr = abs(lon-clon)+abs(lat-clat)
deb = np.argmin(dr)-200
fin = np.argmin(dr)+200

#sigma0 = db2linear(dset.variables['sigma0'][:])
sigma0 = (dset.variables['sigma0'][:])
sigma0_nv = dset.variables['sigma0_nv'][:]
#sigma0_trend = db2linear(dset.variables['sigma0_trend'][:])
# sigma0_trend = (dset.variables['sigma0_trend'][:])
# sigma0 = sigma0 /sigma0_trend

seg_flag = dset.variables['seg_flag'][:]

phi = dset.variables['phi'][:]
phi_geo = dset.variables['phi_geo'][:]

klin = dset.variables['klin'][:]
k = dset.variables['k'][:]
dk = dset.variables['dk'][:]

inc = dset.variables['incidence'][:]
l2s_inc = dset.l2s_angle

u10 = dset.variables['seg_model_u10'][:].mean(axis=1)
v10 = dset.variables['seg_model_v10'][:].mean(axis=1)
wind = np.sqrt(u10**2+v10**2)

mtf = dset.variables['mtf'][:].mean(axis=1)

partition_label = dset.variables['partition_label'][:]
partition_wavelength = dset.variables['partition_wavelength'][:]
partition_direction = dset.variables['partition_direction'][:]
partition_lon = dset.variables['partition_lon'][:]
partition_lat = dset.variables['partition_lat'][:]
ptime =  dset.variables['partition_time'][:]
partition_hs = dset.variables['partition_hs'][:]
tcs = dset.getncattr('time_coverage_start')
tce = dset.getncattr('time_coverage_end')

# srt = ptime.argsort()

# partition_lon = partition_lon[srt]
# partition_lat = partition_lat[srt]
# partition_direction = partition_direction[srt]
# partition_wavelength = partition_wavelength[srt]
# partition_hs = partition_hs[srt]

drpart = abs(partition_lon-clon)+abs(partition_lat-clat)
npart0 = np.argmin(drpart)-100

dset.close()

fluct_spec = fluct_specs.mean(axis=-1)
mod_spec = mod_specs.mean(axis=-1)

#sli = slice(deb,fin)#lat.shape[0]

sli = slice(np.clip(deb,0,None),fin)

Plot your SWIM data

# Plot your data to check that it is ok
fig = plt.figure(figsize = (20,10))
w = np.where(partition_wavelength>40)
plt.plot(partition_lat[w],partition_wavelength[w],'.')
plt.grid()

Compute Center points of coordinates

The near and far coordinates are given in this L2S file, we need to recompute the center of the measurement to provide to SEAScope the proper point


_nk = 28 # fluc_spec
print(klin[_nk - 1], 2. * np.pi / klin[_nk - 1])
_lat = np.zeros((fluct_spec[sli, :].shape[0], _nk))
_lon = np.zeros((fluct_spec[sli, :].shape[0], _nk))
_lat1 =  np.zeros((fluct_spec[sli, :].shape[0], _nk))
_lat2 =  np.zeros((fluct_spec[sli, :].shape[0], _nk))
_lon1 =  np.zeros((fluct_spec[sli, :].shape[0], _nk))
_lon2 =  np.zeros((fluct_spec[sli, :].shape[0], _nk))

geod= pyproj.Geod(ellps='WGS84')
for i in range(_lat.shape[0]):
    az,back_azimuth,dist = geod.inv(near_lon[sli][i], near_lat[sli][i], far_lon[sli][i], far_lat[sli][i])
    endlon, endlat, backaz = geod.fwd(near_lon[sli][i], near_lat[sli][i], az+90, dist/60)
    dlon = endlon - near_lon[sli][i]
    dlat = endlat - near_lat[sli][i]
    ll = geod.npts(near_lon[sli][i], near_lat[sli][i], far_lon[sli][i], far_lat[sli][i], npts=_nk)
    _lat[i, :] = np.array(ll)[:, 1]
    _lon[i, :] = np.array(ll)[:, 0]
    ll1 = geod.npts(near_lon[sli][i]-dlon, near_lat[sli][i]-dlat, far_lon[sli][i]-dlon, far_lat[sli][i]-dlat, npts=_nk)
    _lat1[i, :] = np.array(ll1)[:, 1]
    _lon1[i, :] = np.array(ll1)[:, 0]    
    ll2 = geod.npts(near_lon[sli][i]+dlon, near_lat[sli][i]+dlat, far_lon[sli][i]+dlon, far_lat[sli][i]+dlat, npts=_nk)
    _lat2[i, :] = np.array(ll2)[:, 1]
    _lon2[i, :] = np.array(ll2)[:, 0]

Remove low frequencies for the fluctation and modulation spectrum

Low frequencies are artefacts and not wave signals

def remove_low_freq(spec, ikstart):
    nphi, nk = spec.shape
    ikstop = np.zeros(nphi, dtype='int32') + ikstart
    for ik in range(ikstart, nk - 2):
        ind = np.where(ikstop == ik)[0]
        if ind.size == 0:
            break
        lowf = spec[ind, ik] > spec[ind, ik + 1]
        if lowf.min() == False:
            indnl = ind[~lowf]
            spec[indnl, 0:ik + 1] = 0
        if lowf.max() == True:
            indl = ind[lowf]
            ikstop[indl] = ik + 1
    return spec


# Remove low frequencies artefacts
ikstart = np.where(klin <= 2. * np.pi / 1000.)[0][-1]
flsp = remove_low_freq(fluct_spec[deb:deb+_lat.shape[0], :], ikstart)
modsp = remove_low_freq(mod_spec[deb:deb+_lat.shape[0], :], ikstart)

Export your data to SEAScope

Your collection is named CFOSAT SWIM spectra v1.0

Your variable is named beam + deg modulation spec swell v1 (beam is the beam value angle you have chosen in the beginning)

# import necessary SEAScope librairy
from SEAScope.lib.utils import create_collection, init_ids
from SEAScope.lib.utils import create_granule
import SEAScope.upload
from SEAScope.lib.utils import set_field

from SEAScope.lib.utils import create_variable

# IP address and port used to reach the SEAScope standalone application
host = '127.0.0.1'
port = 11155
# Create collection
collection_id, collection = create_collection('CFOSAT SWIM spectra v1.0')

with SEAScope.upload.connect(host, port) as link:
    # Send create collection job
    SEAScope.upload.collection(link, collection)
# Upload granules in variable  deg modulation spec swell v1  
with SEAScope.upload.connect(host, port) as link:   
    for i in range(_lat.shape[0]):
        if (i%50)==0: print(i)
        gcps = []
        gcps1 = []
        gcps2 = []
        count = 0
        for j in range(_nk):
            gcp = {'lon': _lon[i,j], 'lat': _lat[i,j], 'i': i, 'j': j}
            gcp1 = {'lon': _lon1[i,j], 'lat': _lat1[i,j], 'i': i, 'j': j}
            gcp2 = {'lon': _lon2[i,j], 'lat': _lat2[i,j], 'i': i, 'j': j}
            if _lat[i,j] != 0:
                gcps.append(gcp)
                gcps1.append(gcp1)
                gcps2.append(gcp2)
                count=count+1

        # Create granule
        start_dt = datetime.datetime.strptime(tcs, "%Y-%m-%dT%H:%M:%SZ")
        stop_dt = datetime.datetime.strptime(tce, "%Y-%m-%dT%H:%M:%SZ")
        granule_id, granule = create_granule(collection_id, gcps, start_dt, stop_dt)
        granule_id1, granule1 = create_granule(collection_id, gcps1, start_dt, stop_dt)
        granule_id2, granule2 = create_granule(collection_id, gcps2, start_dt, stop_dt)

        # Set variables
        # Enter a name for each variable field
#         field_name_time = beam+' deg fluc spec swell v1'
        field_name_time = f'{beam} deg modulation spec swell v1'
#         set_field(granule, field_name_time, fluct_spec[deb+i, 0:_nk])
#         SEAScope.upload.granule(link, granule)
#         set_field(granule1, field_name_time, fluct_spec[deb+i, 0:_nk])
#         SEAScope.upload.granule(link, granule1)
#         set_field(granule2, field_name_time, fluct_spec[deb+i, 0:_nk])
#         SEAScope.upload.granule(link, granule2)
        set_field(granule, field_name_time, modsp[i, 0:_nk])
        SEAScope.upload.granule(link, granule)
        set_field(granule1, field_name_time, modsp[i, 0:_nk])
        SEAScope.upload.granule(link, granule1)
        set_field(granule2, field_name_time, modsp[i, 0:_nk])
        SEAScope.upload.granule(link, granule2)
#         field_name_time = beam+' deg wave spec swell'
# #         field_name_time = beam+' deg partition index'
#         set_field(granule, field_name_time, wave_spec[deb+i, 0:_nk])
#         set_field(granule, field_name_time, partition_label[deb+i, 0:_nk])
        
        # Send create granule job    
# Create a variable for each parameter
variable_name_time = field_name_time
var1 = create_variable(collection, variable_name_time, [field_name_time], dims=1)
    
# Rendering configuration for fluctuation spectrum
rcfg1 = var1['rendering']
rcfg1['min'] = 0
rcfg1['max'] = 0.2 #for fluct spec
rcfg1['colormap'] = 'jet'
rcfg1['color'] = [255,255,255]
rcfg1['opacity'] = 1
rcfg1['zindex'] = .923

# Upload to SEAScope
with SEAScope.upload.connect(host, port) as link:
    # Send create variable jobs
    SEAScope.upload.variable(link, var1)
    # Send create rendering config jobs
    SEAScope.upload.rendering_config(link, rcfg1)
# Upload granules in wavenumber ticks to show the ticks
with SEAScope.upload.connect(host, port) as link:   
    for i in range(_lat.shape[0]):
        nt = [3, 6, 13, 26]
        for j in range(4):
            gcps = []
            az,back_azimuth,dist = geod.inv(_lon[i,nt[j]], _lat[i,nt[j]], _lon[i,nt[j]+1], _lat[i,nt[j]+1])
            endlon, endlat, backaz = geod.fwd(_lon[i,nt[j]], _lat[i,nt[j]], az+90, dist/2)
            dlon = endlon - _lon[i,nt[j]]
            dlat = endlat - _lat[i,nt[j]]
            gcp1 = {'lon': _lon[i,nt[j]]-dlon, 'lat': _lat[i,nt[j]]-dlat, 'i': i, 'j': 0}
            gcp2 = {'lon': _lon[i,nt[j]]+dlon, 'lat': _lat[i,nt[j]]+dlat, 'i': i, 'j': 1}
            if _lat[i,j] != 0:
                gcps.append(gcp1)
                gcps.append(gcp2)

        # Create granule
            start_dt = datetime.datetime.strptime(tcs, "%Y-%m-%dT%H:%M:%SZ")
            stop_dt = datetime.datetime.strptime(tce, "%Y-%m-%dT%H:%M:%SZ")
            granule_id, granule = create_granule(collection_id, gcps, start_dt, stop_dt)

        # Set variables
        # Enter a name for each variable field
            field_name_time = f'{beam} wavenumber ticks'
            set_field(granule, field_name_time,np.ones(2))
#         field_name_time = beam+' deg wave spec swell'
# #         field_name_time = beam+' deg partition index'
#         set_field(granule, field_name_time, wave_spec[deb+i, 0:_nk])
#         set_field(granule, field_name_time, partition_label[deb+i, 0:_nk])
        
        # Send create granule job
            SEAScope.upload.granule(link, granule)
    
# Create a variable for each parameter
variable_name_time = field_name_time
var1 = create_variable(collection, variable_name_time, [field_name_time], dims=1)

# Upload to SEAScope
with SEAScope.upload.connect(host, port) as link:
    # Send create variable jobs
    SEAScope.upload.variable(link, var1)
    
# Rendering configuration for hs
rcfg1 = var1['rendering']
rcfg1['min'] = 1
rcfg1['max'] = 1 #for fluct spec
rcfg1['colormap'] = ''
rcfg1['color'] = [255,255,255]
rcfg1['opacity'] = 0.5
rcfg1['zindex'] = .953

# Upload to SEAScope
with SEAScope.upload.connect(host, port) as link:
    # Send create rendering config jobs
    SEAScope.upload.rendering_config(link, rcfg1)

Go back to SEAScope and look for your data

Your collection is named CFOSAT SWIM spectral v1.0 Your variables are named

  • beam + deg modulation spec swell v1

  • beam + wavenumber ticks (beam is the beam value angle you have chosen in the beginning)

You can start over at the cell where you define your beam angle and choose an other one to process

OVL Manual
Converting data with the IDF converter
OVL Manual
Computing Sentinel-1 radial velocities