GridMapping instance

In [1]:

import dask.array as da
import numpy as np
import pyproj as pp
import xarray as xr

import dask.array as da import numpy as np import pyproj as pp import xarray as xr

In [3]:

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

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

Out[3]:

<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

In [4]:

crs_4326 = pp.crs.CRS(4326)
crs_4326

crs_4326 = pp.crs.CRS(4326) crs_4326

Out[4]:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

In [5]:

crs_4326.is_geographic

crs_4326.is_geographic

Out[5]:

True

In [6]:

crs_4326.is_geocentric

crs_4326.is_geocentric

Out[6]:

False

In [7]:

crs_32633 = pp.crs.CRS(32633)
crs_32633

crs_32633 = pp.crs.CRS(32633) crs_32633

Out[7]:

<Projected CRS: EPSG:32633>
Name: WGS 84 / UTM zone 33N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: 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.
- bounds: (12.0, 0.0, 18.0, 84.0)
Coordinate Operation:
- name: UTM zone 33N
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

In [8]:

t = pp.transformer.Transformer.from_crs(crs_4326, crs_32633)
t

t = pp.transformer.Transformer.from_crs(crs_4326, crs_32633) t

Out[8]:

<Conversion Transformer: pipeline>
Description: UTM zone 33N
Area of Use:
- name: Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore.
- bounds: (12.0, 0.0, 18.0, 84.0)

In [9]:

lon = ds.lon.values
lat = ds.lat.values

lon = ds.lon.values lat = ds.lat.values

In [10]:

%%time
x, y = t.transform(lat, lon)

%%time x, y = t.transform(lat, lon)

CPU times: user 454 ms, sys: 10.1 ms, total: 464 ms
Wall time: 463 ms

In [11]:

x, y

x, y

Out[11]:

(array([[464353.17350437, 464753.99068697, 464753.99068697, ...,
         778366.98327271, 778663.97989779, 778960.82095417],
        [464300.61922255, 464697.86336861, 464697.86336861, ...,
         778311.52553234, 778608.48730308, 778905.40483205],
        [464244.50782782, 464643.5297553 , 464643.5297553 , ...,
         778256.11946402, 778553.04638088, 778849.92904976],
        ...,
        [355582.31759961, 356009.5675282 , 356009.5675282 , ...,
         668917.33742756, 669212.77356855, 669508.22326567],
        [355522.80213548, 355950.07579013, 355950.07579013, ...,
         668857.14267177, 669152.60061632, 669448.00001141],
        [355463.34306198, 355890.64413473, 355890.64413473, ...,
         668796.87484134, 669092.350345  , 669387.77154332]],
       shape=(1890, 1189)),
 array([[6722665.0879521 , 6722587.51099095, 6722587.51099095, ...,
         6662383.08675777, 6662325.42149084, 6662267.76727586],
        [6722380.12777498, 6722303.24876914, 6722303.24876914, ...,
         6662095.72015417, 6662038.04563426, 6661980.39060315],
        [6722095.87204537, 6722018.63563135, 6722018.63563135, ...,
         6661808.24764388, 6661750.56387706, 6661692.89959194],
        ...,
        [6186558.12992815, 6186468.08028835, 6186468.08028835, ...,
         6120997.19266252, 6120934.81417697, 6120872.34396918],
        [6186274.18860283, 6186184.02228292, 6186184.02228292, ...,
         6120710.71380218, 6120648.21975818, 6120585.85404474],
        [6185990.24598573, 6185900.07425244, 6185900.07425244, ...,
         6120424.23357397, 6120361.73520792, 6120299.25393304]],
       shape=(1890, 1189)))

In [12]:

coords_4326 = da.stack([ds.lat.data, ds.lon.data]).rechunk((2, None, None))
coords_4326

coords_4326 = da.stack([ds.lat.data, ds.lon.data]).rechunk((2, None, None)) coords_4326

Out[12]:

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

In [13]:

def transform(block, t=None):
    x1, y1 = block
    x2, y2 = t.transform(x1, y1)
    return np.stack([x2, y2])

def transform(block, t=None): x1, y1 = block x2, y2 = t.transform(x1, y1) return np.stack([x2, y2])

In [14]:

coords_32633 = da.apply_gufunc(
    transform, "()->()", coords_4326, output_dtypes=np.float64, t=t
)
coords_32633

coords_32633 = da.apply_gufunc( transform, "()->()", coords_4326, output_dtypes=np.float64, t=t ) coords_32633

Out[14]:

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

In [15]:

%%time
x, y = coords_32633.compute()

%%time x, y = coords_32633.compute()

CPU times: user 955 ms, sys: 107 ms, total: 1.06 s
Wall time: 343 ms

In [16]:

x, y

x, y

Out[16]:

(array([[464353.17350437, 464753.99068697, 464753.99068697, ...,
         778366.98327271, 778663.97989779, 778960.82095417],
        [464300.61922255, 464697.86336861, 464697.86336861, ...,
         778311.52553234, 778608.48730308, 778905.40483205],
        [464244.50782782, 464643.5297553 , 464643.5297553 , ...,
         778256.11946402, 778553.04638088, 778849.92904976],
        ...,
        [355582.31759961, 356009.5675282 , 356009.5675282 , ...,
         668917.33742756, 669212.77356855, 669508.22326567],
        [355522.80213548, 355950.07579013, 355950.07579013, ...,
         668857.14267177, 669152.60061632, 669448.00001141],
        [355463.34306198, 355890.64413473, 355890.64413473, ...,
         668796.87484134, 669092.350345  , 669387.77154332]],
       shape=(1890, 1189)),
 array([[6722665.0879521 , 6722587.51099095, 6722587.51099095, ...,
         6662383.08675777, 6662325.42149084, 6662267.76727586],
        [6722380.12777498, 6722303.24876914, 6722303.24876914, ...,
         6662095.72015417, 6662038.04563426, 6661980.39060315],
        [6722095.87204537, 6722018.63563135, 6722018.63563135, ...,
         6661808.24764388, 6661750.56387706, 6661692.89959194],
        ...,
        [6186558.12992815, 6186468.08028835, 6186468.08028835, ...,
         6120997.19266252, 6120934.81417697, 6120872.34396918],
        [6186274.18860283, 6186184.02228292, 6186184.02228292, ...,
         6120710.71380218, 6120648.21975818, 6120585.85404474],
        [6185990.24598573, 6185900.07425244, 6185900.07425244, ...,
         6120424.23357397, 6120361.73520792, 6120299.25393304]],
       shape=(1890, 1189)))

In [17]:

def transform2(x1, y1, t=None):
    x2, y2 = t.transform(x1, y1)
    return x2, y2

def transform2(x1, y1, t=None): x2, y2 = t.transform(x1, y1) return x2, y2

In [18]:

x, y = da.apply_gufunc(
    transform2,
    "(),()->(),()",
    ds.lat.data,
    ds.lon.data,
    output_dtypes=(np.float64, np.float64),
    t=t,
)
x, y

x, y = da.apply_gufunc( transform2, "(),()->(),()", ds.lat.data, ds.lon.data, output_dtypes=(np.float64, np.float64), t=t, ) x, y

Out[18]:

(dask.array<transpose, shape=(1890, 1189), dtype=float64, chunksize=(512, 512), chunktype=numpy.ndarray>,
 dask.array<transpose, shape=(1890, 1189), dtype=float64, chunksize=(512, 512), chunktype=numpy.ndarray>)

In [19]:

%%time
xx = x.compute()
yy = y.compute()

%%time xx = x.compute() yy = y.compute()

CPU times: user 1.69 s, sys: 64.8 ms, total: 1.76 s
Wall time: 570 ms

In [20]:

xx, yy

xx, yy

Out[20]:

(array([[464353.17350437, 464753.99068697, 464753.99068697, ...,
         778366.98327271, 778663.97989779, 778960.82095417],
        [464300.61922255, 464697.86336861, 464697.86336861, ...,
         778311.52553234, 778608.48730308, 778905.40483205],
        [464244.50782782, 464643.5297553 , 464643.5297553 , ...,
         778256.11946402, 778553.04638088, 778849.92904976],
        ...,
        [355582.31759961, 356009.5675282 , 356009.5675282 , ...,
         668917.33742756, 669212.77356855, 669508.22326567],
        [355522.80213548, 355950.07579013, 355950.07579013, ...,
         668857.14267177, 669152.60061632, 669448.00001141],
        [355463.34306198, 355890.64413473, 355890.64413473, ...,
         668796.87484134, 669092.350345  , 669387.77154332]],
       shape=(1890, 1189)),
 array([[6722665.0879521 , 6722587.51099095, 6722587.51099095, ...,
         6662383.08675777, 6662325.42149084, 6662267.76727586],
        [6722380.12777498, 6722303.24876914, 6722303.24876914, ...,
         6662095.72015417, 6662038.04563426, 6661980.39060315],
        [6722095.87204537, 6722018.63563135, 6722018.63563135, ...,
         6661808.24764388, 6661750.56387706, 6661692.89959194],
        ...,
        [6186558.12992815, 6186468.08028835, 6186468.08028835, ...,
         6120997.19266252, 6120934.81417697, 6120872.34396918],
        [6186274.18860283, 6186184.02228292, 6186184.02228292, ...,
         6120710.71380218, 6120648.21975818, 6120585.85404474],
        [6185990.24598573, 6185900.07425244, 6185900.07425244, ...,
         6120424.23357397, 6120361.73520792, 6120299.25393304]],
       shape=(1890, 1189)))

In [21]:

pp.crs.CRS.from_string("urn:ogc:def:crs:OGC:1.3:CRS84")

pp.crs.CRS.from_string("urn:ogc:def:crs:OGC:1.3:CRS84")

Out[21]:

<Geographic 2D CRS: OGC:CRS84>
Name: WGS 84 (CRS84)
Axis Info [ellipsoidal]:
- Lon[east]: Geodetic longitude (degree)
- Lat[north]: Geodetic latitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

In [22]:

pp.crs.CRS.from_user_input("epsg:4326")

pp.crs.CRS.from_user_input("epsg:4326")

Out[22]:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich