Example for Sea Surface Calibration (SSC) and Sea Surface Leveling (SSL)

The Sea Surface Leveling (SSL) has been introduced by Rott et al. 2022 to align scanning lidars in offshore wind farms.

The Sea Surface Calibration (SSC) builds on this work to calibration the elevation offsets of lidar scanners from multi-elevation scans of the sea surface. Further, it corrects the lidar-sea surface range determination.

Lets apply both to an example dataset, that has been uploaded on Zenodo with the contribution Meyer et al. 2026.

from lidalign.SSC import SSC
import numpy as np
import pandas as pd
import xarray as xr 
import plotly.graph_objects as go
import plotly.express as px
from plotly.subplots import make_subplots
from itertools import cycle
from lidalign.SSC import db2linear, linear2db, WaterRangeDetection, GaussianTruncatedPulse

Visualisation of data

Let us visualise our dataset:

fig = px.imshow(ds['cnr'].T, origin = 'lower', color_continuous_scale = 'viridis', zmax = 0, zmin = -28)
fig.update_layout(yaxis_title = 'Radial distance [m]',  xaxis_title = 'Time', title_text = 'Carrier to noise ratio (CNR) as function of the range over time')
fig.show(renderer= 'notebook')

fig = go.Figure()
fig.add_trace(go.Scattergl(x = ds['time'], y = ds['elevation']))
fig.update_layout(yaxis_title = 'Elevation [deg]',  xaxis_title = 'Time', title_text = 'Elevation of the scans over time')
fig.show(renderer = 'notebook')

You can clearly see, that with increased negative elevation angle, the CNR drop to the noise level (~-28dB) happens at shorter ranges. With only one elevation angle, we can only determine the external misalignments. For the internal misalignment (elevation offset), multiple elevations must be considered.

Determination of Lidar-sea surface distance

For each Line of Sight, a distance to the sea surface from the lidar needs to be determined. Some methods exist for this, which we want to visualize here. All of them are based on the work by Rott et al. 2022.

🔔 Important: We highly recommend the use of the LinSig or Convo method, as the other methods lead to faulty lidar-sea surface ranges, which lead to significant errors in the obtained static elevation offset.

Method Name	Reference	Description
Rot22	Rott et al 2022	Inverse sigmoid, initial method for following work
Gra24	Gramitzky et al. 2024	Extension of inverse sigmoid for linear (in dB scale) decay/increas of the CNR due to atmospheric backscatter at longer ranges
Gra24-PV2	Gramitzky et al. 2025	Observed possible error in range determination, removed half of the probe volume length from the obtained lidar-sea surface range.
LinSig	Meyer et al. 2026 (in preparation)	Conversion of the inverse sigmoid into linear scale of CNR and fit of function in linear scale
dBSig	Meyer et al. 2026 (in preparation)	Conversion of the inverse sigmoid into linear scale of CNR and fit of function in dB scale
Convo	Meyer et al. 2026 (in preparation)	Convolution approach. The probe length is here given as fixed input and does not need to be fitted anymore

The LinSig method uses: $$ S_{\rm sig, lin} (r) = \begin{cases} 10^{\left(A \cdot (r - r_{\rm w})\right)/10} \cdot S_{\rm sig}, & r< r_{\rm w} \ S_{\rm sig}, & \rm{else}. \end{cases}, \mathrm{with} \ S_{\rm sig} = \frac{S_{\rm high}}{1 +\exp{(g \cdot (r-r_{\rm w}))}}. $$

For more details on the methods see the description in Meyer et al. 2026.

Conversion between dB and linear scale

Whereas the carrier to noise ratio (CNR) normally is provided in decibel (dB), it is measured in linear scale. Thus, the inverse sigmoid decay should describe the CNR in linear scale. Conversion between both is simple:

$$ \begin{equation} S_{\rm dB} = 10 \log_{10}(S_{\rm lin}), \end{equation} $$

For details of the corresponding methods, we refere to the mentioned references.

Examplary determination of the lidar-sea surface range

We demonstrate all methods for a single Line of Sight from our dataset. Shown are the CNR from the measurements and from the fitted models for the CNR in decibel scale (upper plot) and in linear scale (lower plot).

The obtained ranges from the methods are plotted in the corresponding figures as vertical lines.

dst = ds.sel(time = '2025-08-12 05:09', method = 'nearest')
#dst = ds.sel(time = '2025-08-12 05:05:44.83', method = 'nearest')

pulse = GaussianTruncatedPulse(75)
fig = make_subplots(rows = 2, shared_xaxes=True)
kw = dict(mode = 'markers', marker = dict(color = 'black', size = 10), name = 'Measurement', line= dict(color = 'black'), legendgroup= 'Measurement')
fig.add_trace(go.Scattergl(x = dst.range, y = dst['cnr'], **kw), row = 1, col = 1)
fig.add_trace(go.Scattergl(x = dst.range, y = db2linear(dst['cnr']), showlegend= False, **kw), row = 2, col = 1)
fig.update_layout(yaxis_title = 'CNR [dB]', height = 800)
fig.update_yaxes(title = 'CNR [linear]', col= 1, row = 2)
fig.update_xaxes(title = 'Range [m]', col = 1, row = 2)


colors = cycle(px.colors.qualitative.Plotly)
yvals = np.array([dst['cnr'].min(), min(dst['cnr'].max(),10)+3])

for method in ["LinSig", "dBSig", "Gra24", "Rot21","Convo", 'Gra24-PV']:
    WRD = WaterRangeDetection(dst.copy(), pulse= pulse, verbose=2)

    if method == 'Gra24-PV':
        res = WRD.get_water_range_from_cnr(func = 'Gra24', return_fit=True)
        res.r_water -= pulse.FWHM//2
    else:
        res = WRD.get_water_range_from_cnr(func = method, return_fit=True)
    
    kw = dict(name = f"{method}: {res.r_water:.1f}m", legendgroup = method, line = dict(color = next(colors)))
    fig.add_trace(go.Scattergl(x = dst.range, y = res.fit_data['fit_db'], **kw), row =1, col = 1)
    fig.add_trace(go.Scattergl(x = dst.range, y = res.fit_data['fit_lin'],showlegend=False, **kw), row =2, col = 1)

    fig.add_trace(go.Scattergl(x = [res.r_water, res.r_water], y = yvals, mode = 'lines', showlegend=False, **kw), row = 1, col = 1)
    fig.add_trace(go.Scattergl(x = [res.r_water, res.r_water], y = db2linear(yvals), mode = 'lines', showlegend=False, **kw), row = 2, col = 1)
fig.show(renderer = 'notebook')
    

first guess for middle range: 1870
first_guess: [np.int32(1870), np.float64(0.0007379042301291007), np.float64(0.08035261221856173), 0.01, 0.01]
bounds: Bounds(array([ 2.5e+02, -4.0e+01, -2.3e+01,  5.0e-04, -1.5e-02]), array([ 4.96e+03, -2.00e+01, -3.00e+00,  1.00e+00,  3.00e-02]))
First guess 0.0007379042301291007 not in bounds -40.0 - -20.0, clipping
First guess 0.08035261221856173 not in bounds -23.0 - -3.0, clipping
Using fit method: LSQ
first guess for middle range: 1870
first_guess: [np.int32(1870), np.float64(-31.32), np.float64(-10.95), 0.01, 0.01]
bounds: Bounds(array([ 2.5e+02, -4.0e+01, -2.3e+01,  5.0e-04, -1.5e-02]), array([ 4.96e+03, -2.00e+01, -3.00e+00,  1.00e+00,  3.00e-02]))
Using fit method: LSQ
first guess for middle range: 1870
first_guess: [np.int32(1870), np.float64(-31.32), np.float64(-10.95), 0.01, 0.01]
bounds: Bounds(array([ 2.5e+02, -4.0e+01, -2.3e+01,  5.0e-04, -1.5e-02]), array([ 4.96e+03, -2.00e+01, -3.00e+00,  1.00e+00,  3.00e-02]))
Using fit method: LSQ
first guess for middle range: 1870
first_guess: [np.int32(1870), np.float64(-31.32), np.float64(-10.95), 0.01]
bounds: Bounds(array([ 2.5e+02, -4.0e+01, -2.3e+01,  5.0e-04]), array([4960,  -20,   -3,    1]))
Using fit method: LSQ
first guess for middle range: 1870
Had to switch to minimize for convolution fit method, due to fitting problems with curve fit and ODR
first_guess: [np.int32(1870), 0.01, np.float64(-10.95), -32]
bounds: Bounds(array([ 2.5e+02, -1.5e-02, -2.2e+01, -4.0e+01]), array([ 4.96e+03,  3.00e-02,  0.00e+00, -2.50e+01]))
Using fit method: minimize
first guess for middle range: 1870
first_guess: [np.int32(1870), np.float64(-31.32), np.float64(-10.95), 0.01, 0.01]
bounds: Bounds(array([ 2.5e+02, -4.0e+01, -2.3e+01,  5.0e-04, -1.5e-02]), array([ 4.96e+03, -2.00e+01, -3.00e+00,  1.00e+00,  3.00e-02]))
Using fit method: LSQ

We can see the differences in the obtained CNR curves. Whereas the methods (except Rot21) show similar CNR curves, the Linsig, Convo and Gra24-PV yield significantly lower lidar-sea surface ranges.

🔔 Important: Note, that we despite the fact, that the Gra24-PV shows similar results as the LinSig, we recommend using the latter due to its insensitivity to the probe volume definition. See Meyer et al. 2026 for more details.

Parameter explanation

There are a few parameters, that need to be explained for the usage of the scripts. As each lidar system might have different internal noise and measurement capabilities, also the CNR distribution will look different and needs to be adjusted. Also, atmospheric conditions will change the resulting CNR distribution.

The most important parameters are:

Parameter	Recommendation	Explanation
min_cnr	-22 [dB]	CNR of LOS must reach at least this value, to be considered valid (if beam is blocked, CNR is always < than this value)
cnr_hard_target	0 [dB]	CNR of LOS must be lower than this value, otherwise it is considered a hard target
use_bounds	True	Use given bounds for fit (improves fit)
func	See above	see above
growth_rate_bounds	[0.0005, 1]	Fit boundaries for the growth rate. This parameter actually depends on the probe volume and should be adjusted accordingly. Larger probe volumes (~100m) give smaller values, whereas shorte probe volumes result in higher values. The lower bound might need to be adjusted according to your prove volume!
cnr_noise_cut	-40 dB	Noise level for CNR cut. If the cnr goes below this level, all CNR readings afterwards are thrown away, as they are considered as noise. Has an impact on the sigmoid functions and should be used for LinSig (-28 e.g., in the case of a WindCube400S)
_first_guess_roll_window	10	Depending on the radial resolution of your input CNR data, you should adjust this value, which is used to find the steepest CNR drop (of a rolling window with size _first_guess_roll_window) as first guess for r_water

🔔 Important: Due to the nature of the fit functions with multiple parameters, the first guess for the range is of high importance for a good fit. The first guess can also be specified through the parameter *dist_guess*.

There are some more parameters, see the code for explanation.

For the whole scan

Lets perform the detection for the whole scan (and make some adjustments through the kw arguments)

dist_ds = SSC(ds, verbose = 1).get_all_water_ranges(func = 'LinSig',growth_rate_bounds = [0.01, 1], min_cnr = -20, cnr_hard_target = -3).distance_ds

fig = px.imshow(ds['cnr'].T, origin = 'lower', color_continuous_scale = 'viridis', zmax = 0, zmin = -28)

fig.add_trace(go.Scattergl(x = dist_ds.time, y = dist_ds['water_range'], name = 'Obtained Lidar-sea surface range', mode = 'markers', showlegend=True))
fig.update_layout(legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.05
))
fig.show(renderer= 'notebook')

As we can see, the obtained ranges follow nicely the CNR cut off. You should always somehow make sure, they align. Otherwise, your following steps might result in erroneous parameters.

Sea Surface Leveling and Calibration

SSC and SSL have different purposes:

SSL is used to align your scanning lidar offshore. Its evaluation and the application of the scans is quick and can be done during the installation while staff is on site to correct possible external misalignment.
SSC is used to calibrate the elevation offset of the lidar scanner. This process takes a bit more time and can be done fully remotely, as the correction must be applied in the software/stearing of the device.

Lets perform SSL and SSC with the obtained lidar-sea surface ranges, the print_help keyword translates the misalignments into actual rotations required on the lidar legs (in case of a WindCubeScan) that would be required to align the lidar properly.

res = SSC.get_misalignment(dist_ds, print_help = True)

Pitch: 0.012°
Roll: -0.015°
This means: 
 	- the lidar west side (both legs) must go up by 0.15 rotations of the leg screws
 	- the lidar north side (both legs) must go down by 0.12 rotations of the leg screws

Lets plot the final results and see, how well the model fits the actual obtained lidar-sea surface ranges:

disti = SSC.rotated_water_range(dist_ds["elevation"], dist_ds["azimuth"], *res.x)

fig = go.Figure()
fig.add_trace(go.Scattergl(x = dist_ds['azimuth'],y = dist_ds['water_range'], name= 'Measurements', mode = 'markers', marker = dict(color = 'black', size = 10)))
fig.add_trace(go.Scattergl(x = dist_ds['azimuth'],y = disti, name= 'Model', mode = 'markers', marker = dict(color = 'red')))

fig.add_annotation(
    x=0.02,
    y=0.8,
    xref="paper",
    yref="paper",
    align="left",
    showarrow=False,
    bgcolor="rgba(255,255,255,0.8)",
    bordercolor="black",
    borderwidth=1,
    borderpad=6,
    font=dict(size=12),
    text=(
        "<b>SSC Fit: </b> <br>"
        f"<b>Pitch:</b> {res.result_dict['pitch']:.3f}°<br>"
        f"<b>Roll:</b> {res.result_dict['roll']:.3f}°<br>"
        f"<b>Δφ:</b> {res.result_dict['ele_offset']:.3f}"
    )
)

fig.update_layout(yaxis_title = 'Lidar-sea surface range [m]', xaxis_title = 'Lidar Azimuth [deg]')
fig.show(renderer = 'notebook')

Application to RHI data:

The method is not limited to the application to multi-elevation PPI data, but can also be applied to RHI measurements. If multiple RHIs at different azimuths are measured, the same code as above can be applied. If only one RHI is available, the tilt in that direction must be known exactly in advance an then could be used to limit the variable in the bounds keyword.