{ "cells": [ { "cell_type": "markdown", "id": "93d1b56e", "metadata": {}, "source": [ "# Example for Sea Surface Calibration (SSC) and Sea Surface Leveling (SSL)\n", "\n", "The **Sea Surface Leveling (SSL)** has been introduced by [Rott et al. 2022](https://wes.copernicus.org/articles/7/283/2022/) to align scanning lidars in offshore wind farms. \n", "\n", "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. \n" ] }, { "cell_type": "markdown", "id": "f54381c6", "metadata": {}, "source": [ "Lets apply both to an example dataset, that has been [uploaded on Zenodo](https://zenodo.org/records/18698332) with the contribution Meyer et al. 2026. " ] }, { "cell_type": "code", "execution_count": null, "id": "a0d7a569", "metadata": {}, "outputs": [], "source": [ "from lidalign.SSC import SSC\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr \n", "import plotly.graph_objects as go\n", "import plotly.express as px\n", "from plotly.subplots import make_subplots\n", "from itertools import cycle\n", "from lidalign.SSC import db2linear, linear2db, WaterRangeDetection, GaussianTruncatedPulse\n" ] }, { "cell_type": "markdown", "id": "712c199b", "metadata": {}, "source": [ "### Example Dataset\n", "We have provided one multi-elevation scan for the Sea Surface Calibration, that has been obtained offshore by a scanning wind lidar (WindCubeScan 400s) at three negative elevation angles. For other applications, see the scripts of the validation campaign in Heligoland in this repository." ] }, { "cell_type": "code", "execution_count": 5, "id": "a9eeefb0", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 4MB\n",
       "Dimensions:    (time: 1740, range: 315)\n",
       "Coordinates:\n",
       "  * range      (range) int32 1kB 250 265 280 295 310 ... 4915 4930 4945 4960\n",
       "  * time       (time) datetime64[ns] 14kB 2025-08-12T05:04:03.036000 ... 2025...\n",
       "Data variables:\n",
       "    azimuth    (time) float64 14kB ...\n",
       "    elevation  (time) float64 14kB ...\n",
       "    cnr        (time, range) float64 4MB ...
" ], "text/plain": [ " Size: 4MB\n", "Dimensions: (time: 1740, range: 315)\n", "Coordinates:\n", " * range (range) int32 1kB 250 265 280 295 310 ... 4915 4930 4945 4960\n", " * time (time) datetime64[ns] 14kB 2025-08-12T05:04:03.036000 ... 2025...\n", "Data variables:\n", " azimuth (time) float64 14kB ...\n", " elevation (time) float64 14kB ...\n", " cnr (time, range) float64 4MB ..." ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds = xr.open_dataset('/home/paul/Downloads/Offshore_SSC_Scan.nc') # you will need to adjust the path here.\n", "ds" ] }, { "cell_type": "markdown", "id": "5a0ce770", "metadata": {}, "source": [ "### Visualisation of data\n", "Let us visualise our dataset:" ] }, { "cell_type": "code", "execution_count": 7, "id": "9d6ae860", "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = px.imshow(ds['cnr'].T, origin = 'lower', color_continuous_scale = 'viridis', zmax = 0, zmin = -28)\n", "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')\n", "fig.show(renderer= 'notebook')\n", "\n", "fig = go.Figure()\n", "fig.add_trace(go.Scattergl(x = ds['time'], y = ds['elevation']))\n", "fig.update_layout(yaxis_title = 'Elevation [deg]', xaxis_title = 'Time', title_text = 'Elevation of the scans over time')\n", "fig.show(renderer = 'notebook')" ] }, { "cell_type": "markdown", "id": "bf32b141", "metadata": {}, "source": [ "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. " ] }, { "cell_type": "markdown", "id": "fea5b56b", "metadata": {}, "source": [ "## Determination of Lidar-sea surface distance \n", "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](https://wes.copernicus.org/articles/7/283/2022/wes-7-283-2022.html).\n", "
\n", " 🔔 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. \n", "
\n", "\n", "\n", "|Method Name|Reference|Description|\n", "|---|---|---|\n", "|Rot22|[Rott et al 2022](https://wes.copernicus.org/articles/7/283/2022/wes-7-283-2022.html)| Inverse sigmoid, initial method for following work|\n", "|Gra24|[Gramitzky et al. 2024](https://iopscience.iop.org/article/10.1088/1742-6596/2767/4/042014)| Extension of inverse sigmoid for linear (in dB scale) decay/increas of the CNR due to atmospheric backscatter at longer ranges|\n", "|Gra24-PV2|[Gramitzky et al. 2025](https://wes.copernicus.org/articles/11/861/2026/wes-11-861-2026.html)|Observed possible error in range determination, removed half of the probe volume length from the obtained lidar-sea surface range. |\n", "|LinSig|Meyer et al. 2026 (in preparation)|Conversion of the inverse sigmoid into linear scale of CNR and fit of function in linear scale| \n", "|dBSig|Meyer et al. 2026 (in preparation)|Conversion of the inverse sigmoid into linear scale of CNR and fit of function in **dB** scale| \n", "|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| \n", "\n", "The *LinSig* method uses:\n", "$$\n", "S_{\\rm sig, lin} (r) = \n", "\\begin{cases}\n", " 10^{\\left(A \\cdot (r - r_{\\rm w})\\right)/10} \\cdot \n", " S_{\\rm sig}, & r< r_{\\rm w} \\\\\n", " S_{\\rm sig}, & \\rm{else}.\n", "\\end{cases}, \\mathrm{with} \\ S_{\\rm sig} = \\frac{S_{\\rm high}}{1 +\\exp{(g \\cdot (r-r_{\\rm w}))}}.\n", "$$\n", "\n", "For more details on the methods see the description in Meyer et al. 2026.\n", "### Conversion between dB and linear scale\n", "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:\n", "\n", "$$\n", " \\begin{equation}\n", " S_{\\rm dB} = 10 \\log_{10}(S_{\\rm lin}), \n", "\\end{equation}\n", "$$\n", "\n", "\n", "For details of the corresponding methods, we refere to the mentioned references." ] }, { "cell_type": "markdown", "id": "0648a3a6", "metadata": {}, "source": [ "### Examplary determination of the lidar-sea surface range\n", "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). \n", "\n", "The obtained ranges from the methods are plotted in the corresponding figures as vertical lines. " ] }, { "cell_type": "code", "execution_count": null, "id": "dd4df005", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "first guess for middle range: 1870\n", "first_guess: [np.int32(1870), np.float64(0.0007379042301291007), np.float64(0.08035261221856173), 0.01, 0.01]\n", "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]))\n", "First guess 0.0007379042301291007 not in bounds -40.0 - -20.0, clipping\n", "First guess 0.08035261221856173 not in bounds -23.0 - -3.0, clipping\n", "Using fit method: LSQ\n", "first guess for middle range: 1870\n", "first_guess: [np.int32(1870), np.float64(-31.32), np.float64(-10.95), 0.01, 0.01]\n", "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]))\n", "Using fit method: LSQ\n", "first guess for middle range: 1870\n", "first_guess: [np.int32(1870), np.float64(-31.32), np.float64(-10.95), 0.01, 0.01]\n", "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]))\n", "Using fit method: LSQ\n", "first guess for middle range: 1870\n", "first_guess: [np.int32(1870), np.float64(-31.32), np.float64(-10.95), 0.01]\n", "bounds: Bounds(array([ 2.5e+02, -4.0e+01, -2.3e+01, 5.0e-04]), array([4960, -20, -3, 1]))\n", "Using fit method: LSQ\n", "first guess for middle range: 1870\n", "Had to switch to minimize for convolution fit method, due to fitting problems with curve fit and ODR\n", "first_guess: [np.int32(1870), 0.01, np.float64(-10.95), -32]\n", "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]))\n", "Using fit method: minimize\n", "first guess for middle range: 1870\n", "first_guess: [np.int32(1870), np.float64(-31.32), np.float64(-10.95), 0.01, 0.01]\n", "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]))\n", "Using fit method: LSQ\n" ] }, { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dst = ds.sel(time = '2025-08-12 05:09', method = 'nearest')\n", "#dst = ds.sel(time = '2025-08-12 05:05:44.83', method = 'nearest')\n", "\n", "pulse = GaussianTruncatedPulse(75)\n", "fig = make_subplots(rows = 2, shared_xaxes=True)\n", "kw = dict(mode = 'markers', marker = dict(color = 'black', size = 10), name = 'Measurement', line= dict(color = 'black'), legendgroup= 'Measurement')\n", "fig.add_trace(go.Scattergl(x = dst.range, y = dst['cnr'], **kw), row = 1, col = 1)\n", "fig.add_trace(go.Scattergl(x = dst.range, y = db2linear(dst['cnr']), showlegend= False, **kw), row = 2, col = 1)\n", "fig.update_layout(yaxis_title = 'CNR [dB]', height = 800)\n", "fig.update_yaxes(title = 'CNR [linear]', col= 1, row = 2)\n", "fig.update_xaxes(title = 'Range [m]', col = 1, row = 2)\n", "\n", "\n", "colors = cycle(px.colors.qualitative.Plotly)\n", "yvals = np.array([dst['cnr'].min(), min(dst['cnr'].max(),10)+3])\n", "\n", "for method in [\"LinSig\", \"dBSig\", \"Gra24\", \"Rot21\",\"Convo\", 'Gra24-PV']:\n", " WRD = WaterRangeDetection(dst.copy(), pulse= pulse, verbose=2)\n", "\n", " if method == 'Gra24-PV':\n", " res = WRD.get_water_range_from_cnr(func = 'Gra24', return_fit=True)\n", " res.r_water -= pulse.FWHM//2\n", " else:\n", " res = WRD.get_water_range_from_cnr(func = method, return_fit=True)\n", " \n", " kw = dict(name = f\"{method}: {res.r_water:.1f}m\", legendgroup = method, line = dict(color = next(colors)))\n", " fig.add_trace(go.Scattergl(x = dst.range, y = res.fit_data['fit_db'], **kw), row =1, col = 1)\n", " fig.add_trace(go.Scattergl(x = dst.range, y = res.fit_data['fit_lin'],showlegend=False, **kw), row =2, col = 1)\n", "\n", " fig.add_trace(go.Scattergl(x = [res.r_water, res.r_water], y = yvals, mode = 'lines', showlegend=False, **kw), row = 1, col = 1)\n", " fig.add_trace(go.Scattergl(x = [res.r_water, res.r_water], y = db2linear(yvals), mode = 'lines', showlegend=False, **kw), row = 2, col = 1)\n", "fig.show(renderer = 'notebook')\n", " \n" ] }, { "cell_type": "markdown", "id": "434a7a91", "metadata": {}, "source": [ "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. " ] }, { "cell_type": "markdown", "id": "e0df664b", "metadata": {}, "source": [ "
\n", " 🔔 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. \n", "
\n", "\n", "\n", "### Parameter explanation\n", "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. \n", "\n", "The most important parameters are:\n", "|Parameter|Recommendation|Explanation|\n", "|---|---|---|\n", "|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)|\n", "|cnr_hard_target|0 [dB]| CNR of LOS must be lower than this value, otherwise it is considered a hard target|\n", "|use_bounds|True| Use given bounds for fit (improves fit)|\n", "|func| See above | see above|\n", "|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!|\n", "|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)|\n", "|_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|\n", "\n", "
\n", " 🔔 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*.\n", "
\n", "\n", "\n", "\n", "There are some more parameters, see the code for explanation. " ] }, { "cell_type": "markdown", "id": "e8da9231", "metadata": {}, "source": [ "## For the whole scan\n", "Lets perform the detection for the whole scan (and make some adjustments through the kw arguments)\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "f5f14693", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Obtaining lidar-water range: 0%| | 0/1740 [00:00\n", " window.PlotlyConfig = {MathJaxConfig: 'local'};\n", " if (window.MathJax && window.MathJax.Hub && window.MathJax.Hub.Config) {window.MathJax.Hub.Config({SVG: {font: \"STIX-Web\"}});}\n", " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = px.imshow(ds['cnr'].T, origin = 'lower', color_continuous_scale = 'viridis', zmax = 0, zmin = -28)\n", "\n", "fig.add_trace(go.Scattergl(x = dist_ds.time, y = dist_ds['water_range'], name = 'Obtained Lidar-sea surface range', mode = 'markers', showlegend=True))\n", "fig.update_layout(legend=dict(\n", " yanchor=\"top\",\n", " y=0.99,\n", " xanchor=\"left\",\n", " x=0.05\n", "))\n", "fig.show(renderer= 'notebook')" ] }, { "cell_type": "markdown", "id": "6e6af5b0", "metadata": {}, "source": [ "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. " ] }, { "cell_type": "markdown", "id": "ccfab8bf", "metadata": {}, "source": [ "## Sea Surface Leveling and Calibration\n", "SSC and SSL have different purposes:\n", "- **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. \n", "- **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. \n", "\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 19, "id": "4aaa640d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Pitch: 0.012°\n", "Roll: -0.015°\n", "This means: \n", " \t- the lidar west side (both legs) must go up by 0.15 rotations of the leg screws\n", " \t- the lidar north side (both legs) must go down by 0.12 rotations of the leg screws\n" ] } ], "source": [ "res = SSC.get_misalignment(dist_ds, print_help = True)" ] }, { "cell_type": "markdown", "id": "695e1b01", "metadata": {}, "source": [ "Lets plot the final results and see, how well the model fits the actual obtained lidar-sea surface ranges:" ] }, { "cell_type": "code", "execution_count": 20, "id": "8db193cf", "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "disti = SSC.rotated_water_range(dist_ds[\"elevation\"], dist_ds[\"azimuth\"], *res.x)\n", "\n", "fig = go.Figure()\n", "fig.add_trace(go.Scattergl(x = dist_ds['azimuth'],y = dist_ds['water_range'], name= 'Measurements', mode = 'markers', marker = dict(color = 'black', size = 10)))\n", "fig.add_trace(go.Scattergl(x = dist_ds['azimuth'],y = disti, name= 'Model', mode = 'markers', marker = dict(color = 'red')))\n", "\n", "fig.add_annotation(\n", " x=0.02,\n", " y=0.8,\n", " xref=\"paper\",\n", " yref=\"paper\",\n", " align=\"left\",\n", " showarrow=False,\n", " bgcolor=\"rgba(255,255,255,0.8)\",\n", " bordercolor=\"black\",\n", " borderwidth=1,\n", " borderpad=6,\n", " font=dict(size=12),\n", " text=(\n", " \"SSC Fit:
\"\n", " f\"Pitch: {res.result_dict['pitch']:.3f}°
\"\n", " f\"Roll: {res.result_dict['roll']:.3f}°
\"\n", " f\"Δφ: {res.result_dict['ele_offset']:.3f}\"\n", " )\n", ")\n", "\n", "fig.update_layout(yaxis_title = 'Lidar-sea surface range [m]', xaxis_title = 'Lidar Azimuth [deg]')\n", "fig.show(renderer = 'notebook')" ] }, { "cell_type": "markdown", "id": "45f711e7", "metadata": {}, "source": [ "## Application to RHI data:\n", "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. " ] } ], "metadata": { "kernelspec": { "display_name": ".general_venv", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }