COG multi-band example

import rioxarray

# from https://openaerialmap.org/
cog_url = (
    "https://oin-hotosm.s3.amazonaws.com/5d7dad0becaf880008a9bc88/0/5d7dad0becaf880008a9bc89.tif"
)

Open a pre-computed lower resolution overview

• overview_level=4 is an argument specific to the rasterio engine that allows opening a pre-computed lower resolution “overview” of the data

dset = rioxarray.open_rasterio(cog_url, masked=True, overview_level=4)

dset

<xarray.DataArray (band: 3, y: 312, x: 688)>
[643968 values with dtype=float32]
Coordinates:
  * band         (band) int64 1 2 3
  * x            (x) float64 -1.333e+07 -1.333e+07 ... -1.333e+07 -1.333e+07
  * y            (y) float64 4.34e+06 4.34e+06 4.34e+06 ... 4.339e+06 4.339e+06
    spatial_ref  int64 0
Attributes:
    acquisitionEndDate:    2019-09-14T17:10:25+00:00
    acquisitionStartDate:  2019-09-14T16:40:25+00:00
    AREA_OR_POINT:         Area
    isCalibrated:          False
    scale_factor:          1.0
    add_offset:            0.0

xarray.DataArray

• band: 3
• y: 312
• x: 688

• ...

  [643968 values with dtype=float32]

• Coordinates: (4)
  • band
    (band)
    int64
    1 2 3

    array([1, 2, 3])

  • x
    (x)
    float64
    -1.333e+07 ... -1.333e+07

    array([-13334018.583529, -13334017.389201, -13334016.194872, ...,
           -13333200.468461, -13333199.274132, -13333198.079804])

  • y
    (y)
    float64
    4.34e+06 4.34e+06 ... 4.339e+06

    array([4339530.145784, 4339528.951456, 4339527.757127, ..., 4339161.098257,
           4339159.903929, 4339158.7096  ])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

    spatial_ref :

    PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

    GeoTransform :

    -13334019.180693429 1.194328566955879 0.0 4339530.742948703 0.0 -1.194328566955879

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index([1, 2, 3], dtype='int64', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([-13334018.583529146,  -13334017.38920058, -13334016.194872012,
           -13334015.000543445, -13334013.806214878, -13334012.611886311,
           -13334011.417557744, -13334010.223229177,  -13334009.02890061,
           -13334007.834572043,
           ...
            -13333208.82876075, -13333207.634432184, -13333206.440103617,
            -13333205.24577505, -13333204.051446483, -13333202.857117916,
           -13333201.662789349, -13333200.468460781, -13333199.274132214,
           -13333198.079803647],
          dtype='float64', name='x', length=688))

  • y
    PandasIndex

    PandasIndex(Index([4339530.145784419, 4339528.951455852, 4339527.757127285,
           4339526.562798718, 4339525.368470151, 4339524.174141584,
           4339522.979813017,  4339521.78548445, 4339520.591155883,
           4339519.396827317,
           ...
           4339169.458557199, 4339168.264228632, 4339167.069900065,
           4339165.875571498, 4339164.681242931, 4339163.486914364,
           4339162.292585797,  4339161.09825723, 4339159.903928663,
           4339158.709600096],
          dtype='float64', name='y', length=312))

• Attributes: (6)

  acquisitionEndDate :

  2019-09-14T17:10:25+00:00

  acquisitionStartDate :

  2019-09-14T16:40:25+00:00

  AREA_OR_POINT :

  Area

  isCalibrated :

  False

  scale_factor :

  1.0

  add_offset :

  0.0

dset.sel(band=1).plot(cmap="Greys_r")

<matplotlib.collections.QuadMesh at 0x7f26bc81e550>

Visualize as RGB image

dset.astype("int").plot.imshow(rgb="band")

/home/runner/micromamba/envs/jupytergis-examples/lib/python3.11/site-packages/xarray/core/duck_array_ops.py:188: RuntimeWarning: invalid value encountered in cast
  return data.astype(dtype, **kwargs)
/home/runner/micromamba/envs/jupytergis-examples/lib/python3.11/site-packages/xarray/plot/utils.py:253: RuntimeWarning: overflow encountered in scalar absolute
  vlim = max(abs(vmin - center), abs(vmax - center))
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

<matplotlib.image.AxesImage at 0x7f26be72d810>

Open with Xarray

• We need to open and read a data array

• we can pass an argument to open a pre-computed lower resolution overview

import xarray as xr

dset = xr.open_dataarray(cog_url, engine='rasterio', open_kwargs={"overview_level": 2})
dset

<xarray.DataArray 'band_data' (band: 3, y: 1248, x: 2752)>
[10303488 values with dtype=float32]
Coordinates:
  * band         (band) int64 1 2 3
  * x            (x) float64 -1.333e+07 -1.333e+07 ... -1.333e+07 -1.333e+07
  * y            (y) float64 4.34e+06 4.34e+06 4.34e+06 ... 4.339e+06 4.339e+06
    spatial_ref  int64 ...
Attributes:
    acquisitionEndDate:    2019-09-14T17:10:25+00:00
    acquisitionStartDate:  2019-09-14T16:40:25+00:00
    AREA_OR_POINT:         Area
    isCalibrated:          False

xarray.DataArray

'band_data'

• band: 3
• y: 1248
• x: 2752

• ...

  [10303488 values with dtype=float32]

• Coordinates: (4)
  • band
    (band)
    int64
    1 2 3

    array([1, 2, 3])

  • x
    (x)
    float64
    -1.333e+07 ... -1.333e+07

    array([-13334019.031402, -13334018.73282 , -13334018.434238, ...,
           -13333198.229095, -13333197.930513, -13333197.63193 ])

  • y
    (y)
    float64
    4.34e+06 4.34e+06 ... 4.339e+06

    array([4339530.593658, 4339530.295075, 4339529.996493, ..., 4339158.858891,
           4339158.560309, 4339158.261727])

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

    spatial_ref :

    PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

    GeoTransform :

    -13334019.180693429 0.29858214173896974 0.0 4339530.742948703 0.0 -0.29858214173896974

    [1 values with dtype=int64]

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index([1, 2, 3], dtype='int64', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([-13334019.031402359, -13334018.732820217, -13334018.434238076,
           -13334018.135655934, -13334017.837073792,  -13334017.53849165,
            -13334017.23990951, -13334016.941327367, -13334016.642745225,
           -13334016.344163083,
           ...
           -13333200.319169711,  -13333200.02058757, -13333199.722005427,
           -13333199.423423285, -13333199.124841144, -13333198.826259002,
            -13333198.52767686, -13333198.229094718, -13333197.930512577,
           -13333197.631930435],
          dtype='float64', name='x', length=2752))

  • y
    PandasIndex

    PandasIndex(Index([ 4339530.593657631,   4339530.29507549,  4339529.996493348,
            4339529.697911207,  4339529.399329064,  4339529.100746923,
            4339528.802164781,   4339528.50358264,  4339528.205000497,
            4339527.906418356,
           ...
            4339160.948966159,  4339160.650384017,  4339160.351801875,
            4339160.053219734,  4339159.754637592,   4339159.45605545,
            4339159.157473308,  4339158.858891167, 4339158.5603090245,
            4339158.261726883],
          dtype='float64', name='y', length=1248))

• Attributes: (4)

  acquisitionEndDate :

  2019-09-14T17:10:25+00:00

  acquisitionStartDate :

  2019-09-14T16:40:25+00:00

  AREA_OR_POINT :

  Area

  isCalibrated :

  False

dset.astype("int").plot.imshow(rgb="band")

/home/runner/micromamba/envs/jupytergis-examples/lib/python3.11/site-packages/xarray/core/duck_array_ops.py:188: RuntimeWarning: invalid value encountered in cast
  return data.astype(dtype, **kwargs)
/home/runner/micromamba/envs/jupytergis-examples/lib/python3.11/site-packages/xarray/plot/utils.py:253: RuntimeWarning: overflow encountered in scalar absolute
  vlim = max(abs(vmin - center), abs(vmax - center))
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

<matplotlib.image.AxesImage at 0x7f26bc745810>