Demonstration of function resample_in_space()¶

This notebook demonstrates the functionality of the resample_in_space function. We show how it can be used for three different purposes:

• Rectification
• Reprojection to another Coordinate Reference System
• Up- and Downsampling

In [1]:

import xarray as xr
from xcube_resampling.gridmapping import GridMapping
from xcube_resampling import resample_in_space
from zarr.storage import ZipStore

import xarray as xr from xcube_resampling.gridmapping import GridMapping from xcube_resampling import resample_in_space from zarr.storage import ZipStore

First we open a Sentinel-3 OLCI Level-1B product.

In [2]:

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

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

CPU times: user 42.7 ms, sys: 2.9 ms, total: 45.6 ms
Wall time: 44.8 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

• 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

Next we can plot the raw data.

In [3]:

%%time
source_ds.rtoa_8.plot(vmax=0.25)

%%time source_ds.rtoa_8.plot(vmax=0.25)

CPU times: user 661 ms, sys: 132 ms, total: 793 ms
Wall time: 734 ms

Out[3]:

<matplotlib.collections.QuadMesh at 0x725a8c7e0d70>

The dataset uses 2D coordinates lon and lat, each defined over two dimensions (x and y) to represent the satellite viewing geometry. These coordinates are expressed in the EPSG:4326 geographic coordinate system (WGS84).

Rectification¶

We can create a GridMapping object from the dataset, which contains all the information required by the rectification algorithm.

In [4]:

%%time
source_gm = GridMapping.from_dataset(source_ds)
source_gm

%%time source_gm = GridMapping.from_dataset(source_ds) source_gm

CPU times: user 55.5 ms, sys: 5.13 ms, total: 60.6 ms
Wall time: 50.2 ms

Out[4]:

class: Coords2DGridMapping

• is_regular: False
• is_j_axis_up: False
• is_lon_360: False
• crs: EPSG:4326
• xy_res: (0.00447132, 0.00255235) estimated
• xy_bbox: (12.6920653, 55.1989038, 20.0071197, 60.6394662)
• ij_bbox: (0, 0, 1189, 1890)
• xy_dim_names: ('x', 'y')
• xy_var_names: ('lon', 'lat')
• size: (1189, 1890)
• tile_size: (512, 512)

We aim to rectify this dataset so that it uses a regular EPSG:4326 CRS with 1D coordinates lon and lat having constant spacing. To do this, we first create a new regular target grid mapping as follows:

In [5]:

%%time
target_gm = source_gm.to_regular(tile_size=512)
target_gm

%%time target_gm = source_gm.to_regular(tile_size=512) target_gm

CPU times: user 1.1 ms, sys: 16 μs, total: 1.12 ms
Wall time: 650 μs

Out[5]:

class: RegularGridMapping

• is_regular: True
• is_j_axis_up: False
• is_lon_360: False
• crs: EPSG:4326
• xy_res: (0.00447132, 0.00255235)
• xy_bbox: (12.6920653, 55.1989038, 20.0071448, 60.640514)
• ij_bbox: (0, 0, 1636, 2132)
• xy_dim_names: ('lon', 'lat')
• xy_var_names: ('lon', 'lat')
• size: (1636, 2132)
• tile_size: (512, 512)

The target grid mapping is created lazily: the large 2D x and y coordinates used for rectification are stored as Dask arrays.

In [6]:

%%time
x_coords, y_coords = target_gm.xy_coords
y_coords

%%time x_coords, y_coords = target_gm.xy_coords y_coords

CPU times: user 3.4 ms, sys: 0 ns, total: 3.4 ms
Wall time: 3.4 ms

Out[6]:

<xarray.DataArray 'xy_coords' (lat: 2132, lon: 1636)> Size: 28MB
dask.array<getitem, shape=(2132, 1636), dtype=float64, chunksize=(512, 512), chunktype=numpy.ndarray>
Dimensions without coordinates: lat, lon

xarray.DataArray

'xy_coords'

• lat: 2132
• lon: 1636

• dask.array<chunksize=(512, 512), meta=np.ndarray>

  Array Chunk Bytes 26.61 MiB 2.00 MiB Shape (2132, 1636) (512, 512) Dask graph 20 chunks in 10 graph layers Data type float64 numpy.ndarray

We can now perform resample_in_space to rectify the dataset.

In [7]:

%%time
target_ds = resample_in_space(source_ds, source_gm=source_gm, target_gm=target_gm, interp_methods="nearest")
target_ds

%%time target_ds = resample_in_space(source_ds, source_gm=source_gm, target_gm=target_gm, interp_methods="nearest") target_ds

CPU times: user 860 ms, sys: 55.5 ms, total: 915 ms
Wall time: 793 ms

Out[7]:

<xarray.Dataset> Size: 56MB
Dimensions:        (lon: 1636, lat: 2132)
Coordinates:
  * 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
    spatial_ref    int64 8B 0
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)
  • lon
    (lon)
    float64
    12.69 12.7 12.7 ... 20.0 20.0 20.0

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

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

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

  • 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)

• 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

• 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 [8]:

%%time
target_ds.rtoa_8.plot(vmax=0.25)

%%time target_ds.rtoa_8.plot(vmax=0.25)

CPU times: user 12.4 s, sys: 3.08 s, total: 15.4 s
Wall time: 2.03 s

Out[8]:

<matplotlib.collections.QuadMesh at 0x725a74f0a850>

Reprojection¶

We can also switch to a different coordinate reference system, for example UTM zone 33N, also known as EPSG:32633.

Since we need the bounding box and the target resolution in UTM coordinates, we first transform the source grid mapping. This produces an irregular grid mapping with 2D UTM coordinates.

In [9]:

%%time
temp_utm_target_gm = source_gm.transform("EPSG:32633")
temp_utm_target_gm

%%time temp_utm_target_gm = source_gm.transform("EPSG:32633") temp_utm_target_gm

CPU times: user 338 ms, sys: 6.62 ms, total: 345 ms
Wall time: 176 ms

Out[9]:

class: Coords2DGridMapping

• is_regular: False
• is_j_axis_up: False
• is_lon_360: True
• crs: EPSG:32633
• xy_res: (264.534, 285.506) estimated
• xy_bbox: (355331, 6120157, 779093, 6722808)
• ij_bbox: (0, 0, 1189, 1890)
• xy_dim_names: ('x', 'y')
• xy_var_names: ('transformed_x', 'transformed_y')
• size: (1189, 1890)
• tile_size: (512, 512)

In the next cell, we convert the irregular grid into a regular grid:

In [10]:

%%time
utm_target_gm = temp_utm_target_gm.to_regular()
utm_target_gm

%%time utm_target_gm = temp_utm_target_gm.to_regular() utm_target_gm

CPU times: user 92 μs, sys: 0 ns, total: 92 μs
Wall time: 95.1 μs

Out[10]:

class: RegularGridMapping

• is_regular: True
• is_j_axis_up: False
• is_lon_360: False
• crs: EPSG:32633
• xy_res: (264.534, 285.506)
• xy_bbox: (355330.733, 6120157.247, 779114.267, 6722860.753)
• ij_bbox: (0, 0, 1602, 2111)
• xy_dim_names: ('x', 'y')
• xy_var_names: ('x', 'y')
• size: (1602, 2111)
• tile_size: (512, 512)

In [11]:

%%time
x_coords, y_coords = utm_target_gm.xy_coords
y_coords

%%time x_coords, y_coords = utm_target_gm.xy_coords y_coords

CPU times: user 3.1 ms, sys: 23 μs, total: 3.12 ms
Wall time: 3.13 ms

Out[11]:

<xarray.DataArray 'xy_coords' (y: 2111, x: 1602)> Size: 27MB
dask.array<getitem, shape=(2111, 1602), dtype=float64, chunksize=(512, 512), chunktype=numpy.ndarray>
Dimensions without coordinates: y, x

xarray.DataArray

'xy_coords'

• y: 2111
• x: 1602

• dask.array<chunksize=(512, 512), meta=np.ndarray>

  Array Chunk Bytes 25.80 MiB 2.00 MiB Shape (2111, 1602) (512, 512) Dask graph 20 chunks in 10 graph layers Data type float64 numpy.ndarray

Next, we apply resample_in_space to resample the dataset onto the newly defined target grid.

In [12]:

%%time
transformed_target_ds = resample_in_space(
    source_ds, source_gm=source_gm, target_gm=utm_target_gm, interp_methods=0
)
transformed_target_ds

%%time transformed_target_ds = resample_in_space( source_ds, source_gm=source_gm, target_gm=utm_target_gm, interp_methods=0 ) transformed_target_ds

CPU times: user 2.41 s, sys: 82.7 ms, total: 2.49 s
Wall time: 1.14 s

Out[12]:

<xarray.Dataset> Size: 54MB
Dimensions:        (x: 1602, y: 2111)
Coordinates:
  * x              (x) float64 13kB 3.555e+05 3.557e+05 ... 7.787e+05 7.79e+05
  * y              (y) float64 17kB 6.723e+06 6.722e+06 ... 6.121e+06 6.12e+06
    spatial_ref    int64 8B 0
Data variables:
    quality_flags  (y, x) uint32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_3         (y, x) float32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_6         (y, x) float32 14MB dask.array<chunksize=(512, 512), meta=np.ndarray>
    rtoa_8         (y, x) 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:
  • x: 1602
  • y: 2111
• Coordinates: (3)
  • x
    (x)
    float64
    3.555e+05 3.557e+05 ... 7.79e+05

    array([355463.      , 355727.534041, 355992.068082, ..., 778452.931918,
           778717.465959, 778982.      ], shape=(1602,))

  • y
    (y)
    float64
    6.723e+06 6.722e+06 ... 6.12e+06

    array([6722718.      , 6722432.493839, 6722146.987678, ..., 6120871.012322,
           6120585.506161, 6120300.      ], shape=(2111,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCRS["WGS 84 / UTM zone 33N",BASEGEOGCRS["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]],ID["EPSG",4326]],CONVERSION["UTM zone 33N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",15,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Navigation and medium accuracy spatial referencing."],AREA["Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State."],BBOX[0,12,84,18]],ID["EPSG",32633]]

    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

    projected_crs_name :

    WGS 84 / UTM zone 33N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    15.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    array(0)

• 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

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 12.90 MiB 1.00 MiB Shape (2111, 1602) (512, 512) Dask graph 20 chunks in 71 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

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 12.90 MiB 1.00 MiB Shape (2111, 1602) (512, 512) Dask graph 20 chunks in 70 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

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 12.90 MiB 1.00 MiB Shape (2111, 1602) (512, 512) Dask graph 20 chunks in 70 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

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 12.90 MiB 1.00 MiB Shape (2111, 1602) (512, 512) Dask graph 20 chunks in 70 graph layers Data type float32 numpy.ndarray

• 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 [13]:

%%time
transformed_target_ds.rtoa_8.plot(vmax=0.25)

%%time transformed_target_ds.rtoa_8.plot(vmax=0.25)

CPU times: user 8.91 s, sys: 957 ms, total: 9.87 s
Wall time: 2.22 s

Out[13]:

<matplotlib.collections.QuadMesh at 0x725a6385f610>

The next examples show how to we can apply up- and down-sampling.

Downsampling¶

In [14]:

%%time
downsampled_target_gm = target_gm.scale(0.25, tile_size=128)
downsampled_target_gm

%%time downsampled_target_gm = target_gm.scale(0.25, tile_size=128) downsampled_target_gm

CPU times: user 99 μs, sys: 14 μs, total: 113 μs
Wall time: 116 μs

Out[14]:

class: RegularGridMapping

• is_regular: True
• is_j_axis_up: False
• is_lon_360: False
• crs: EPSG:4326
• xy_res: (0.01788528, 0.0102094)
• xy_bbox: (12.6920653, 55.1989038, 20.0071448, 60.640514)
• ij_bbox: (0, 0, 409, 533)
• xy_dim_names: ('lon', 'lat')
• xy_var_names: ('lon', 'lat')
• size: (409, 533)
• tile_size: (128, 128)

In [15]:

%%time
x_coords, y_coords = downsampled_target_gm.xy_coords
y_coords

%%time x_coords, y_coords = downsampled_target_gm.xy_coords y_coords

CPU times: user 3.27 ms, sys: 41 μs, total: 3.31 ms
Wall time: 3.33 ms

Out[15]:

<xarray.DataArray 'xy_coords' (lat: 533, lon: 409)> Size: 2MB
dask.array<getitem, shape=(533, 409), dtype=float64, chunksize=(128, 128), chunktype=numpy.ndarray>
Dimensions without coordinates: lat, lon

xarray.DataArray

'xy_coords'

• lat: 533
• lon: 409

• dask.array<chunksize=(128, 128), meta=np.ndarray>

  Array Chunk Bytes 1.66 MiB 128.00 kiB Shape (533, 409) (128, 128) Dask graph 20 chunks in 10 graph layers Data type float64 numpy.ndarray

In [16]:

%%time
downsampled_ds = resample_in_space(
    target_ds, source_gm=target_gm, target_gm=downsampled_target_gm
)
downsampled_ds

%%time downsampled_ds = resample_in_space( target_ds, source_gm=target_gm, target_gm=downsampled_target_gm ) downsampled_ds

CPU times: user 200 ms, sys: 4.24 ms, total: 204 ms
Wall time: 199 ms

Out[16]:

<xarray.Dataset> Size: 3MB
Dimensions:        (lat: 533, lon: 409)
Coordinates:
  * lat            (lat) float64 4kB 60.64 60.63 60.61 60.6 ... 55.22 55.21 55.2
  * lon            (lon) float64 3kB 12.7 12.72 12.74 12.75 ... 19.96 19.98 20.0
    spatial_ref    int64 8B 0
Data variables:
    quality_flags  (lat, lon) uint32 872kB dask.array<chunksize=(128, 128), meta=np.ndarray>
    rtoa_3         (lat, lon) float32 872kB dask.array<chunksize=(128, 128), meta=np.ndarray>
    rtoa_6         (lat, lon) float32 872kB dask.array<chunksize=(128, 128), meta=np.ndarray>
    rtoa_8         (lat, lon) float32 872kB dask.array<chunksize=(128, 128), 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:
  • lat: 533
  • lon: 409
• Coordinates: (3)
  • lat
    (lat)
    float64
    60.64 60.63 60.61 ... 55.21 55.2

    array([60.635409, 60.6252  , 60.61499 , ..., 55.224427, 55.214218, 55.204009],
          shape=(533,))

  • lon
    (lon)
    float64
    12.7 12.72 12.74 ... 19.98 20.0

    array([12.701008, 12.718893, 12.736778, ..., 19.962432, 19.980317, 19.998202],
          shape=(409,))

  • 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)

• Data variables: (4)
  • quality_flags
    (lat, lon)
    uint32
    dask.array<chunksize=(128, 128), 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 851.55 kiB 64.00 kiB Shape (533, 409) (128, 128) Dask graph 20 chunks in 96 graph layers Data type uint32 numpy.ndarray

  • rtoa_3
    (lat, lon)
    float32
    dask.array<chunksize=(128, 128), 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 851.55 kiB 64.00 kiB Shape (533, 409) (128, 128) Dask graph 20 chunks in 299 graph layers Data type float32 numpy.ndarray

  • rtoa_6
    (lat, lon)
    float32
    dask.array<chunksize=(128, 128), 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 851.55 kiB 64.00 kiB Shape (533, 409) (128, 128) Dask graph 20 chunks in 299 graph layers Data type float32 numpy.ndarray

  • rtoa_8
    (lat, lon)
    float32
    dask.array<chunksize=(128, 128), 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 851.55 kiB 64.00 kiB Shape (533, 409) (128, 128) Dask graph 20 chunks in 299 graph layers Data type float32 numpy.ndarray

• 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 [17]:

%%time
downsampled_ds.rtoa_6.plot(vmax=0.25)

%%time downsampled_ds.rtoa_6.plot(vmax=0.25)

/home/konstantin/micromamba/envs/xcube-multistore/lib/python3.13/site-packages/xcube_resampling/coarsen.py:103: RuntimeWarning: Mean of empty slice
  return nan_reducer(block, axis)

CPU times: user 11.6 s, sys: 1.68 s, total: 13.3 s
Wall time: 2.19 s

Out[17]:

<matplotlib.collections.QuadMesh at 0x725a63369f90>

Upsampling¶

In [18]:

%%time
upsampled_target_gm = target_gm.scale(4, tile_size=1024)
upsampled_target_gm

%%time upsampled_target_gm = target_gm.scale(4, tile_size=1024) upsampled_target_gm

CPU times: user 101 μs, sys: 15 μs, total: 116 μs
Wall time: 118 μs

Out[18]:

class: RegularGridMapping

• is_regular: True
• is_j_axis_up: False
• is_lon_360: False
• crs: EPSG:4326
• xy_res: (0.00111783, 0.0006380875)
• xy_bbox: (12.6920653, 55.1989038, 20.0071448, 60.640514)
• ij_bbox: (0, 0, 6544, 8528)
• xy_dim_names: ('lon', 'lat')
• xy_var_names: ('lon', 'lat')
• size: (6544, 8528)
• tile_size: (1024, 1024)

In [19]:

%%time
x_coords, y_coords = upsampled_target_gm.xy_coords
y_coords

%%time x_coords, y_coords = upsampled_target_gm.xy_coords y_coords

CPU times: user 4.71 ms, sys: 0 ns, total: 4.71 ms
Wall time: 4.69 ms

Out[19]:

<xarray.DataArray 'xy_coords' (lat: 8528, lon: 6544)> Size: 446MB
dask.array<getitem, shape=(8528, 6544), dtype=float64, chunksize=(1024, 1024), chunktype=numpy.ndarray>
Dimensions without coordinates: lat, lon

xarray.DataArray

'xy_coords'

• lat: 8528
• lon: 6544

• dask.array<chunksize=(1024, 1024), meta=np.ndarray>

  Array Chunk Bytes 425.78 MiB 8.00 MiB Shape (8528, 6544) (1024, 1024) Dask graph 63 chunks in 10 graph layers Data type float64 numpy.ndarray

In [20]:

%%time
upsampled_target_ds = resample_in_space(
    target_ds, source_gm=target_gm, target_gm=upsampled_target_gm
)
upsampled_target_ds

%%time upsampled_target_ds = resample_in_space( target_ds, source_gm=target_gm, target_gm=upsampled_target_gm ) upsampled_target_ds

CPU times: user 91.5 ms, sys: 1.57 ms, total: 93.1 ms
Wall time: 88.6 ms

Out[20]:

<xarray.Dataset> Size: 893MB
Dimensions:        (lat: 8528, lon: 6544)
Coordinates:
  * lat            (lat) float64 68kB 60.64 60.64 60.64 60.64 ... 55.2 55.2 55.2
  * lon            (lon) float64 52kB 12.69 12.69 12.69 ... 20.0 20.01 20.01
    spatial_ref    int64 8B 0
Data variables:
    quality_flags  (lat, lon) uint32 223MB dask.array<chunksize=(1024, 1024), meta=np.ndarray>
    rtoa_3         (lat, lon) float32 223MB dask.array<chunksize=(1024, 1024), meta=np.ndarray>
    rtoa_6         (lat, lon) float32 223MB dask.array<chunksize=(1024, 1024), meta=np.ndarray>
    rtoa_8         (lat, lon) float32 223MB dask.array<chunksize=(1024, 1024), 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:
  • lat: 8528
  • lon: 6544
• Coordinates: (3)
  • lat
    (lat)
    float64
    60.64 60.64 60.64 ... 55.2 55.2

    array([60.640195, 60.639557, 60.638919, ..., 55.200499, 55.199861, 55.199223],
          shape=(8528,))

  • lon
    (lon)
    float64
    12.69 12.69 12.69 ... 20.01 20.01

    array([12.692624, 12.693742, 12.69486 , ..., 20.00435 , 20.005468, 20.006586],
          shape=(6544,))

  • 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)

• Data variables: (4)
  • quality_flags
    (lat, lon)
    uint32
    dask.array<chunksize=(1024, 1024), 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 212.89 MiB 4.00 MiB Shape (8528, 6544) (1024, 1024) Dask graph 63 chunks in 139 graph layers Data type uint32 numpy.ndarray

  • rtoa_3
    (lat, lon)
    float32
    dask.array<chunksize=(1024, 1024), 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 212.89 MiB 4.00 MiB Shape (8528, 6544) (1024, 1024) Dask graph 63 chunks in 139 graph layers Data type float32 numpy.ndarray

  • rtoa_6
    (lat, lon)
    float32
    dask.array<chunksize=(1024, 1024), 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 212.89 MiB 4.00 MiB Shape (8528, 6544) (1024, 1024) Dask graph 63 chunks in 139 graph layers Data type float32 numpy.ndarray

  • rtoa_8
    (lat, lon)
    float32
    dask.array<chunksize=(1024, 1024), 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 212.89 MiB 4.00 MiB Shape (8528, 6544) (1024, 1024) Dask graph 63 chunks in 139 graph layers Data type float32 numpy.ndarray

• 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 [21]:

%%time
upsampled_target_ds.rtoa_6[::4, ::4].plot(vmax=0.25)

%%time upsampled_target_ds.rtoa_6[::4, ::4].plot(vmax=0.25)

CPU times: user 17.9 s, sys: 1.62 s, total: 19.6 s
Wall time: 2.13 s

Out[21]:

<matplotlib.collections.QuadMesh at 0x725a6330b890>