Rectification to a Sentinel-3 OLCI scene¶

In this notebook, we apply the rectification algorithm to a Setninel-3 OLCI Level-1b scene. The Sentinel-3 OLCI Level-1b product comes with 2D coordinates – irregular grids where each pixel has its own pair of x/y orlat/lon coordinates. We will there rectify the dataset to a regular grid.

In [1]:

Copied!





import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from xcube_resampling.rectify import rectify_dataset
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from xcube_resampling.rectify import rectify_dataset

In [2]:

Copied!

%%time
source_ds = xr.open_zarr("./inputdata/S3-OLCI-L2A.zarr.zip", consolidated=False)
source_ds
%%time
source_ds = xr.open_zarr("./inputdata/S3-OLCI-L2A.zarr.zip", consolidated=False)
source_ds

CPU times: user 116 ms, sys: 15.9 ms, total: 132 ms
Wall time: 151 ms

Out[2]:

<xarray.Dataset> Size: 72MB
Dimensions:        (y: 1890, x: 1189)
Coordinates:
    lat            (y, x) float64 18MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    lon            (y, x) float64 18MB dask.array<chunksize=(512, 512), meta=np.ndarray>
Dimensions without coordinates: y, x
Data variables:
    quality_flags  (y, x) uint32 9MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_3         (y, x) float32 9MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_6         (y, x) float32 9MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_8         (y, x) float32 9MB dask.array<chunksize=(512, 512), meta=np.ndarray>
Attributes:
    Conventions:   CF-1.4
    product_type:  C2RCC_OLCI
    start_date:    04-JUL-2018 09:21:55.677316
    stop_date:     04-JUL-2018 09:23:18.811790

The provided coordinates, latitude (lat) and longitude (lon), are terrain-corrected using a digital elevation model (DEM) that approximates the Earth's actual fractal surface. As a result, the gradient vectors do not vary monotonically across the scene. Instead, they reflect the surface roughness captured by the DEM. Areas of high terrain variability are highlighted in the right figure below.

In [3]:

Copied!





%%time
lon_grad = source_ds.lat.differentiate("y").rename("Vertical Gradient")

fig, ax = plt.subplots(1, 2, figsize=(14, 5))
source_ds.rtoa_8.plot(ax=ax[0], cmap="gray", vmax=0.2)
lon_grad.plot(ax=ax[1], cmap="viridis", robust=True)
%%time
lon_grad = source_ds.lat.differentiate("y").rename("Vertical Gradient")

fig, ax = plt.subplots(1, 2, figsize=(14, 5))
source_ds.rtoa_8.plot(ax=ax[0], cmap="gray", vmax=0.2)
lon_grad.plot(ax=ax[1], cmap="viridis", robust=True)

CPU times: user 893 ms, sys: 198 ms, total: 1.09 s
Wall time: 979 ms

Out[3]:

<matplotlib.collections.QuadMesh at 0x760b10ebc910>

No description has been provided for this image

Now, we can rectify the dataset. Note that this can be done with the resample_in_space method.

In [4]:

Copied!

%%time
target_ds = rectify_dataset(source_ds, interp_methods="nearest")
target_ds
%%time
target_ds = rectify_dataset(source_ds, interp_methods="nearest")
target_ds

CPU times: user 2.42 s, sys: 165 ms, total: 2.59 s
Wall time: 1.38 s

Out[4]:

<xarray.Dataset> Size: 102MB
Dimensions:        (lon: 2926, lat: 2177)
Coordinates:
    spatial_ref    int64 8B 0
  * lon            (lon) float64 23kB 12.69 12.7 12.7 12.7 ... 20.0 20.0 20.01
  * lat            (lat) float64 17kB 60.64 60.64 60.64 ... 55.21 55.2 55.2
Data variables:
    quality_flags  (lat, lon) uint32 25MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_3         (lat, lon) float32 25MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_6         (lat, lon) float32 25MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_8         (lat, lon) float32 25MB dask.array<chunksize=(512, 512), meta=np.ndarray>
Attributes:
    Conventions:   CF-1.4
    product_type:  C2RCC_OLCI
    start_date:    04-JUL-2018 09:21:55.677316
    stop_date:     04-JUL-2018 09:23:18.811790

In [5]:

Copied!

%%time
target_ds.rtoa_8.plot(cmap="gray", vmax=0.2)
%%time
target_ds.rtoa_8.plot(cmap="gray", vmax=0.2)

CPU times: user 2.33 s, sys: 384 ms, total: 2.72 s
Wall time: 1.51 s

Out[5]:

<matplotlib.collections.QuadMesh at 0x760b104c3ed0>

You can also change the chunk size in the output dataset with the keyword argument tile_size. Note that interp_methods=0 is identical to interp_methods="nearest".

In [6]:

Copied!

%%time
target_ds = rectify_dataset(source_ds, interp_methods=0, tile_size=2048)
target_ds
%%time
target_ds = rectify_dataset(source_ds, interp_methods=0, tile_size=2048)
target_ds

CPU times: user 2.17 s, sys: 123 ms, total: 2.29 s
Wall time: 1.4 s

Out[6]:

<xarray.Dataset> Size: 102MB
Dimensions:        (lon: 2926, lat: 2177)
Coordinates:
    spatial_ref    int64 8B 0
  * lon            (lon) float64 23kB 12.69 12.7 12.7 12.7 ... 20.0 20.0 20.01
  * lat            (lat) float64 17kB 60.64 60.64 60.64 ... 55.21 55.2 55.2
Data variables:
    quality_flags  (lat, lon) uint32 25MB dask.array<chunksize=(2048, 2048), meta=np.ndarray>
    rtoa_3         (lat, lon) float32 25MB dask.array<chunksize=(2048, 2048), meta=np.ndarray>
    rtoa_6         (lat, lon) float32 25MB dask.array<chunksize=(2048, 2048), meta=np.ndarray>
    rtoa_8         (lat, lon) float32 25MB dask.array<chunksize=(2048, 2048), meta=np.ndarray>
Attributes:
    Conventions:   CF-1.4
    product_type:  C2RCC_OLCI
    start_date:    04-JUL-2018 09:21:55.677316
    stop_date:     04-JUL-2018 09:23:18.811790

In [7]:

Copied!

%%time
target_ds.rtoa_8.plot(cmap="gray", vmax=0.2)
%%time
target_ds.rtoa_8.plot(cmap="gray", vmax=0.2)

CPU times: user 1.5 s, sys: 167 ms, total: 1.67 s
Wall time: 1.31 s

Out[7]:

<matplotlib.collections.QuadMesh at 0x760b1009e710>

We can also apply a different interpolation method, as shown in the next cell.

In [8]:

Copied!

%%time
target_ds_bilin = rectify_dataset(source_ds, interp_methods="bilinear")
target_ds_bilin
%%time
target_ds_bilin = rectify_dataset(source_ds, interp_methods="bilinear")
target_ds_bilin

CPU times: user 2.16 s, sys: 126 ms, total: 2.29 s
Wall time: 1.29 s

Out[8]:

<xarray.Dataset> Size: 102MB
Dimensions:        (lon: 2926, lat: 2177)
Coordinates:
    spatial_ref    int64 8B 0
  * lon            (lon) float64 23kB 12.69 12.7 12.7 12.7 ... 20.0 20.0 20.01
  * lat            (lat) float64 17kB 60.64 60.64 60.64 ... 55.21 55.2 55.2
Data variables:
    quality_flags  (lat, lon) uint32 25MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_3         (lat, lon) float32 25MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_6         (lat, lon) float32 25MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_8         (lat, lon) float32 25MB dask.array<chunksize=(512, 512), meta=np.ndarray>
Attributes:
    Conventions:   CF-1.4
    product_type:  C2RCC_OLCI
    start_date:    04-JUL-2018 09:21:55.677316
    stop_date:     04-JUL-2018 09:23:18.811790

In [9]:

Copied!

%%time
target_ds_bilin.rtoa_8.plot(cmap="gray", vmax=0.2)
%%time
target_ds_bilin.rtoa_8.plot(cmap="gray", vmax=0.2)

CPU times: user 3.04 s, sys: 177 ms, total: 3.22 s
Wall time: 1.48 s

Out[9]:

<matplotlib.collections.QuadMesh at 0x760b084547d0>

Next, we visualize the difference between bilinear and nearest-neighbor interpolation.

In [10]:

Copied!

(target_ds_bilin.rtoa_8 - target_ds.rtoa_8).plot(vmin=-0.02, vmax=+0.02, cmap="RdBu")
(target_ds_bilin.rtoa_8 - target_ds.rtoa_8).plot(vmin=-0.02, vmax=+0.02, cmap="RdBu")

Out[10]:

<matplotlib.collections.QuadMesh at 0x760b10212850>