Rectification to a Sentinel-3 OLCI Scene¶

In this notebook, we apply a rectification algorithm to a Sentinel-3 OLCI Level-1B scene.

The Sentinel-3 OLCI Level-1B product is provided with two-dimensional geolocation coordinates, meaning that each pixel is associated with its own latitude and longitude values. This results in an irregular (curvilinear) grid rather than a regular, evenly spaced grid.

To enable standard geospatial processing and analysis, we rectify the dataset onto a regular grid using the rectification algorithm.

In [1]:

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]:

%%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 107 ms, sys: 20 ms, total: 127 ms
Wall time: 127 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

xarray.Dataset

• Dimensions:
  • y: 1890
  • x: 1189
• Coordinates: (2)
  • lat
    (y, x)
    float64
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    Array Chunk Bytes 17.14 MiB 2.00 MiB Shape (1890, 1189) (512, 512) Dask graph 12 chunks in 2 graph layers Data type float64 numpy.ndarray

  • lon
    (y, x)
    float64
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    Array Chunk Bytes 17.14 MiB 2.00 MiB Shape (1890, 1189) (512, 512) Dask graph 12 chunks in 2 graph layers Data type float64 numpy.ndarray

• Data variables: (4)
  • quality_flags
    (y, x)
    uint32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    flag_coding_name :

    quality_flags

    flag_descriptions :

    flag_masks :

    [-2147483648, 1073741824, 536870912, 268435456, 134217728, 67108864, 33554432, 16777216, 8388608, 4194304, 2097152, 1048576, 524288, 262144, 131072, 65536, 32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1]

    flag_meanings :

    land coastline fresh_inland_water tidal_region bright straylight_risk invalid cosmetic duplicated sun_glint_risk dubious saturated_Oa01 saturated_Oa02 saturated_Oa03 saturated_Oa04 saturated_Oa05 saturated_Oa06 saturated_Oa07 saturated_Oa08 saturated_Oa09 saturated_Oa10 saturated_Oa11 saturated_Oa12 saturated_Oa13 saturated_Oa14 saturated_Oa15 saturated_Oa16 saturated_Oa17 saturated_Oa18 saturated_Oa19 saturated_Oa20 saturated_Oa21

    long_name :

    Classification and quality flags

    Array Chunk Bytes 8.57 MiB 1.00 MiB Shape (1890, 1189) (512, 512) Dask graph 12 chunks in 2 graph layers Data type uint32 numpy.ndarray

  • rtoa_3
    (y, x)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    2.0

    units :

    1

    wavelength :

    442.5

    Array Chunk Bytes 8.57 MiB 1.00 MiB Shape (1890, 1189) (512, 512) Dask graph 12 chunks in 2 graph layers Data type float32 numpy.ndarray

  • rtoa_6
    (y, x)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    5.0

    units :

    1

    wavelength :

    560.0

    Array Chunk Bytes 8.57 MiB 1.00 MiB Shape (1890, 1189) (512, 512) Dask graph 12 chunks in 2 graph layers Data type float32 numpy.ndarray

  • rtoa_8
    (y, x)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    7.0

    units :

    1

    wavelength :

    665.0

    Array Chunk Bytes 8.57 MiB 1.00 MiB Shape (1890, 1189) (512, 512) Dask graph 12 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (0)
• Attributes: (4)

  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 latitude (lat) and longitude (lon) coordinates are terrain-corrected using a digital elevation model (DEM). This means that the geolocation accounts for topographic variations instead of assuming a smooth ellipsoidal Earth surface.

As a result, the coordinate fields are not smoothly varying across the scene. Instead, they reflect terrain-induced distortions introduced by the DEM. Areas with strong topographic variability therefore show stronger spatial irregularities in the coordinate gradients (highlighted in the right figure below).

In [3]:

%%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 763 ms, sys: 333 ms, total: 1.1 s
Wall time: 1.3 s

Out[3]:

<matplotlib.collections.QuadMesh at 0x7c887a9cd810>

For some applications, however, we prefer a regular grid that is not influenced by terrain-related distortions. To achieve this, we rectify the dataset onto a regular grid using the rectify_dataset method.

As interpolation method, we use nearest. Alternatively, we could specify 0, which is shorthand for nearest-neighbor interpolation.

In [4]:

%%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 787 ms, sys: 36.9 ms, total: 824 ms
Wall time: 718 ms

Out[4]:

<xarray.Dataset> Size: 56MB
Dimensions:        (lon: 1636, lat: 2132)
Coordinates:
    spatial_ref    int64 8B 0
  * lon            (lon) float64 13kB 12.69 12.7 12.7 12.71 ... 20.0 20.0 20.0
  * lat            (lat) float64 17kB 60.64 60.64 60.63 ... 55.21 55.2 55.2
Data variables:
    quality_flags  (lat, lon) uint32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_3         (lat, lon) float32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_6         (lat, lon) float32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_8         (lat, lon) float32 14MB 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

xarray.Dataset

• Dimensions:
  • lon: 1636
  • lat: 2132
• Coordinates: (3)
  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    grid_mapping_name :

    latitude_longitude

    array(0)

  • lon
    (lon)
    float64
    12.69 12.7 12.7 ... 20.0 20.0 20.0

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    bounds :

    lon_bnds

    array([12.694301, 12.698772, 12.703244, ..., 19.995967, 20.000438, 20.004909],
          shape=(1636,))

  • lat
    (lat)
    float64
    60.64 60.64 60.63 ... 55.2 55.2

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    bounds :

    lat_bnds

    array([60.639238, 60.636685, 60.634133, ..., 55.205285, 55.202732, 55.20018 ],
          shape=(2132,))

• Data variables: (4)
  • quality_flags
    (lat, lon)
    uint32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    flag_coding_name :

    quality_flags

    flag_descriptions :

    flag_masks :

    [-2147483648, 1073741824, 536870912, 268435456, 134217728, 67108864, 33554432, 16777216, 8388608, 4194304, 2097152, 1048576, 524288, 262144, 131072, 65536, 32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1]

    flag_meanings :

    land coastline fresh_inland_water tidal_region bright straylight_risk invalid cosmetic duplicated sun_glint_risk dubious saturated_Oa01 saturated_Oa02 saturated_Oa03 saturated_Oa04 saturated_Oa05 saturated_Oa06 saturated_Oa07 saturated_Oa08 saturated_Oa09 saturated_Oa10 saturated_Oa11 saturated_Oa12 saturated_Oa13 saturated_Oa14 saturated_Oa15 saturated_Oa16 saturated_Oa17 saturated_Oa18 saturated_Oa19 saturated_Oa20 saturated_Oa21

    long_name :

    Classification and quality flags

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 1.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 71 graph layers Data type uint32 numpy.ndarray

  • rtoa_3
    (lat, lon)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    2.0

    units :

    1

    wavelength :

    442.5

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 1.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 71 graph layers Data type float32 numpy.ndarray

  • rtoa_6
    (lat, lon)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    5.0

    units :

    1

    wavelength :

    560.0

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 1.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 71 graph layers Data type float32 numpy.ndarray

  • rtoa_8
    (lat, lon)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    7.0

    units :

    1

    wavelength :

    665.0

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 1.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 71 graph layers Data type float32 numpy.ndarray

• Indexes: (2)
  • lon
    PandasIndex

    PandasIndex(Index([       12.69430096, 12.698772279987768, 12.703243599975535,
           12.707714919963303,  12.71218623995107, 12.716657559938838,
           12.721128879926605, 12.725600199914373,  12.73007151990214,
           12.734542839889908,
           ...
            19.96466726011009,  19.96913858009786, 19.973609900085627,
           19.978081220073392,  19.98255254006116,  19.98702386004893,
           19.991495180036697, 19.995966500024466,  20.00043782001223,
                  20.00490914],
          dtype='float64', name='lon', length=1636))

  • lat
    PandasIndex

    PandasIndex(Index([      60.639237825,       60.636685475, 60.634133125000005,
                 60.631580775,       60.629028425,       60.626476075,
           60.623923725000004,       60.621371375,       60.618819025,
           60.616266675000006,
           ...
                 55.223151125,       55.220598775,       55.218046425,
           55.215494074999995,       55.212941725,       55.210389375,
           55.207837024999996, 55.205284674999994,       55.202732325,
                 55.200179975],
          dtype='float64', name='lat', length=2132))

• Attributes: (4)

  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]:

%%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 8.43 s, sys: 1.91 s, total: 10.3 s
Wall time: 1.76 s

Out[5]:

<matplotlib.collections.QuadMesh at 0x7c886010e710>

You can also change the chunk size in the output dataset with the keyword argument tile_size.

In [6]:

%%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 380 ms, sys: 29.1 ms, total: 409 ms
Wall time: 349 ms

Out[6]:

<xarray.Dataset> Size: 56MB
Dimensions:        (lon: 1636, lat: 2132)
Coordinates:
    spatial_ref    int64 8B 0
  * lon            (lon) float64 13kB 12.69 12.7 12.7 12.71 ... 20.0 20.0 20.0
  * lat            (lat) float64 17kB 60.64 60.64 60.63 ... 55.21 55.2 55.2
Data variables:
    quality_flags  (lat, lon) uint32 14MB dask.array<chunksize=(2048, 1636), meta=np.ndarray>
    rtoa_3         (lat, lon) float32 14MB dask.array<chunksize=(2048, 1636), meta=np.ndarray>
    rtoa_6         (lat, lon) float32 14MB dask.array<chunksize=(2048, 1636), meta=np.ndarray>
    rtoa_8         (lat, lon) float32 14MB dask.array<chunksize=(2048, 1636), 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

xarray.Dataset

• Dimensions:
  • lon: 1636
  • lat: 2132
• Coordinates: (3)
  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    grid_mapping_name :

    latitude_longitude

    array(0)

  • lon
    (lon)
    float64
    12.69 12.7 12.7 ... 20.0 20.0 20.0

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    bounds :

    lon_bnds

    array([12.694301, 12.698772, 12.703244, ..., 19.995967, 20.000438, 20.004909],
          shape=(1636,))

  • lat
    (lat)
    float64
    60.64 60.64 60.63 ... 55.2 55.2

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    bounds :

    lat_bnds

    array([60.639238, 60.636685, 60.634133, ..., 55.205285, 55.202732, 55.20018 ],
          shape=(2132,))

• Data variables: (4)
  • quality_flags
    (lat, lon)
    uint32
    dask.array<chunksize=(2048, 1636), meta=np.ndarray>

    flag_coding_name :

    quality_flags

    flag_descriptions :

    flag_masks :

    [-2147483648, 1073741824, 536870912, 268435456, 134217728, 67108864, 33554432, 16777216, 8388608, 4194304, 2097152, 1048576, 524288, 262144, 131072, 65536, 32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1]

    flag_meanings :

    land coastline fresh_inland_water tidal_region bright straylight_risk invalid cosmetic duplicated sun_glint_risk dubious saturated_Oa01 saturated_Oa02 saturated_Oa03 saturated_Oa04 saturated_Oa05 saturated_Oa06 saturated_Oa07 saturated_Oa08 saturated_Oa09 saturated_Oa10 saturated_Oa11 saturated_Oa12 saturated_Oa13 saturated_Oa14 saturated_Oa15 saturated_Oa16 saturated_Oa17 saturated_Oa18 saturated_Oa19 saturated_Oa20 saturated_Oa21

    long_name :

    Classification and quality flags

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 12.78 MiB Shape (2132, 1636) (2048, 1636) Dask graph 2 chunks in 17 graph layers Data type uint32 numpy.ndarray

  • rtoa_3
    (lat, lon)
    float32
    dask.array<chunksize=(2048, 1636), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    2.0

    units :

    1

    wavelength :

    442.5

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 12.78 MiB Shape (2132, 1636) (2048, 1636) Dask graph 2 chunks in 17 graph layers Data type float32 numpy.ndarray

  • rtoa_6
    (lat, lon)
    float32
    dask.array<chunksize=(2048, 1636), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    5.0

    units :

    1

    wavelength :

    560.0

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 12.78 MiB Shape (2132, 1636) (2048, 1636) Dask graph 2 chunks in 17 graph layers Data type float32 numpy.ndarray

  • rtoa_8
    (lat, lon)
    float32
    dask.array<chunksize=(2048, 1636), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    7.0

    units :

    1

    wavelength :

    665.0

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 12.78 MiB Shape (2132, 1636) (2048, 1636) Dask graph 2 chunks in 17 graph layers Data type float32 numpy.ndarray

• Indexes: (2)
  • lon
    PandasIndex

    PandasIndex(Index([       12.69430096, 12.698772279987768, 12.703243599975535,
           12.707714919963303,  12.71218623995107, 12.716657559938838,
           12.721128879926605, 12.725600199914373,  12.73007151990214,
           12.734542839889908,
           ...
            19.96466726011009,  19.96913858009786, 19.973609900085627,
           19.978081220073392,  19.98255254006116,  19.98702386004893,
           19.991495180036697, 19.995966500024466,  20.00043782001223,
                  20.00490914],
          dtype='float64', name='lon', length=1636))

  • lat
    PandasIndex

    PandasIndex(Index([      60.639237825,       60.636685475, 60.634133125000005,
                 60.631580775,       60.629028425,       60.626476075,
           60.623923725000004,       60.621371375,       60.618819025,
           60.616266675000006,
           ...
                 55.223151125,       55.220598775,       55.218046425,
           55.215494074999995,       55.212941725,       55.210389375,
           55.207837024999996, 55.205284674999994,       55.202732325,
                 55.200179975],
          dtype='float64', name='lat', length=2132))

• Attributes: (4)

  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]:

%%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.74 s, sys: 740 ms, total: 2.48 s
Wall time: 1.77 s

Out[7]:

<matplotlib.collections.QuadMesh at 0x7c88605d0f50>

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

In [8]:

%%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 842 ms, sys: 26.9 ms, total: 869 ms
Wall time: 761 ms

Out[8]:

<xarray.Dataset> Size: 56MB
Dimensions:        (lon: 1636, lat: 2132)
Coordinates:
    spatial_ref    int64 8B 0
  * lon            (lon) float64 13kB 12.69 12.7 12.7 12.71 ... 20.0 20.0 20.0
  * lat            (lat) float64 17kB 60.64 60.64 60.63 ... 55.21 55.2 55.2
Data variables:
    quality_flags  (lat, lon) uint32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_3         (lat, lon) float32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_6         (lat, lon) float32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_8         (lat, lon) float32 14MB 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

xarray.Dataset

• Dimensions:
  • lon: 1636
  • lat: 2132
• Coordinates: (3)
  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    grid_mapping_name :

    latitude_longitude

    array(0)

  • lon
    (lon)
    float64
    12.69 12.7 12.7 ... 20.0 20.0 20.0

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    bounds :

    lon_bnds

    array([12.694301, 12.698772, 12.703244, ..., 19.995967, 20.000438, 20.004909],
          shape=(1636,))

  • lat
    (lat)
    float64
    60.64 60.64 60.63 ... 55.2 55.2

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    bounds :

    lat_bnds

    array([60.639238, 60.636685, 60.634133, ..., 55.205285, 55.202732, 55.20018 ],
          shape=(2132,))

• Data variables: (4)
  • quality_flags
    (lat, lon)
    uint32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    flag_coding_name :

    quality_flags

    flag_descriptions :

    flag_masks :

    [-2147483648, 1073741824, 536870912, 268435456, 134217728, 67108864, 33554432, 16777216, 8388608, 4194304, 2097152, 1048576, 524288, 262144, 131072, 65536, 32768, 16384, 8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1]

    flag_meanings :

    land coastline fresh_inland_water tidal_region bright straylight_risk invalid cosmetic duplicated sun_glint_risk dubious saturated_Oa01 saturated_Oa02 saturated_Oa03 saturated_Oa04 saturated_Oa05 saturated_Oa06 saturated_Oa07 saturated_Oa08 saturated_Oa09 saturated_Oa10 saturated_Oa11 saturated_Oa12 saturated_Oa13 saturated_Oa14 saturated_Oa15 saturated_Oa16 saturated_Oa17 saturated_Oa18 saturated_Oa19 saturated_Oa20 saturated_Oa21

    long_name :

    Classification and quality flags

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 1.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 71 graph layers Data type uint32 numpy.ndarray

  • rtoa_3
    (lat, lon)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    2.0

    units :

    1

    wavelength :

    442.5

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 1.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 71 graph layers Data type float32 numpy.ndarray

  • rtoa_6
    (lat, lon)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    5.0

    units :

    1

    wavelength :

    560.0

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 1.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 71 graph layers Data type float32 numpy.ndarray

  • rtoa_8
    (lat, lon)
    float32
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    bandwidth :

    10.0

    long_name :

    Top-of-atmosphere reflectance

    spectral_band_index :

    7.0

    units :

    1

    wavelength :

    665.0

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 13.31 MiB 1.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 71 graph layers Data type float32 numpy.ndarray

• Indexes: (2)
  • lon
    PandasIndex

    PandasIndex(Index([       12.69430096, 12.698772279987768, 12.703243599975535,
           12.707714919963303,  12.71218623995107, 12.716657559938838,
           12.721128879926605, 12.725600199914373,  12.73007151990214,
           12.734542839889908,
           ...
            19.96466726011009,  19.96913858009786, 19.973609900085627,
           19.978081220073392,  19.98255254006116,  19.98702386004893,
           19.991495180036697, 19.995966500024466,  20.00043782001223,
                  20.00490914],
          dtype='float64', name='lon', length=1636))

  • lat
    PandasIndex

    PandasIndex(Index([      60.639237825,       60.636685475, 60.634133125000005,
                 60.631580775,       60.629028425,       60.626476075,
           60.623923725000004,       60.621371375,       60.618819025,
           60.616266675000006,
           ...
                 55.223151125,       55.220598775,       55.218046425,
           55.215494074999995,       55.212941725,       55.210389375,
           55.207837024999996, 55.205284674999994,       55.202732325,
                 55.200179975],
          dtype='float64', name='lat', length=2132))

• Attributes: (4)

  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]:

%%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 8.42 s, sys: 1.42 s, total: 9.83 s
Wall time: 1.38 s

Out[9]:

<matplotlib.collections.QuadMesh at 0x7c8831b811d0>

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

In [10]:

(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 0x7c8831a27110>