Clarification about the right use of DataArray and Dataset

[MaxDragonheart]

I've a lot of scenes from Landsat 8 that I manage as DataArray. I put 9 bands from a scene into one DataArray and I set band_name; filename is the band's name:

band_dataarray = band_dataarray.squeeze()
band_dataarray = band_dataarray.drop('band')
band_dataarray = band_dataarray.assign_coords(band_name=filename)
band_dataarray = band_dataarray.expand_dims(dim='band_name')

Below the output from one band:

<xarray.DataArray (band_name: 1, y: 367, x: 727)>
array([[[65535, 65535, 65535, ..., 65535, 65535, 65535],
        [65535, 65535, 65535, ..., 65535, 65535, 65535],
        [65535, 65535, 65535, ..., 65535, 65535, 65535],
        ...,
        [65535, 65535, 65535, ..., 65535, 65535, 65535],
        [65535, 65535, 65535, ..., 65535, 65535, 65535],
        [65535, 65535, 65535, ..., 65535, 65535, 65535]]], dtype=uint16)
Coordinates:
  * x            (x) float64 4.361e+05 4.361e+05 ... 4.579e+05 4.579e+05
  * y            (y) float64 4.525e+06 4.525e+06 ... 4.514e+06 4.514e+06
    spatial_ref  int64 0
  * band_name    (band_name) <U2 'B1'
Attributes:
    STATISTICS_MAXIMUM:        2586
    STATISTICS_MEAN:           481.89105443705
    STATISTICS_MINIMUM:        261
    STATISTICS_STDDEV:         121.53471756427
    STATISTICS_VALID_PERCENT:  36.71
    _FillValue:                65535.0
    scale_factor:              0.0001
    add_offset:                0.0

All bands are stored into a list(band_list) and I use concat to put it all together:

scene_bands = xr.concat(band_list, dim='band_name')

Now I add the sensing date as another dimension:

scene_bands = scene_bands.assign_coords(time=sensing_datetime)
scene_bands = scene_bands.expand_dims(dim='time')

Below my output at this point:

<xarray.DataArray (time: 1, band_name: 9, y: 367, x: 727)>
array([[[[ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         ...,
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.]],

        [[-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         ...,
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.]],

        [[-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         ...,
...
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.]],

        [[ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         ...,
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.]],

        [[ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         ...,
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.]]]],
      dtype=float32)
Coordinates:
  * x            (x) float64 4.361e+05 4.361e+05 ... 4.579e+05 4.579e+05
  * y            (y) float64 4.525e+06 4.525e+06 ... 4.514e+06 4.514e+06
    spatial_ref  int64 0
  * band_name    (band_name) <U3 'B1' 'B10' 'B11' 'B2' 'B3' 'B4' 'B5' 'B6' 'B7'
  * time         (time) datetime64[ns] 2022-06-28
Attributes:
    STATISTICS_MAXIMUM:        2586
    STATISTICS_MEAN:           481.89105443705
    STATISTICS_MINIMUM:        261
    STATISTICS_STDDEV:         121.53471756427
    STATISTICS_VALID_PERCENT:  36.71
    _FillValue:                65535.0
    scale_factor:              0.0001
    add_offset:                0.0

As I said I've a lot of scenes, so I put all the above data into a list and using again concat I obtain my final output:

time_series_scenes = xr.concat(scenes, dim='time')


<xarray.DataArray (time: 113, band_name: 9, y: 367, x: 727)>
array([[[[ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         ...,
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.]],

        [[-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         ...,
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.]],

        [[-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         [-32768., -32768., -32768., ..., -32768., -32768., -32768.],
         ...,
...
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.]],

        [[ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         ...,
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.]],

        [[ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         ...,
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.],
         [ 65535.,  65535.,  65535., ...,  65535.,  65535.,  65535.]]]],
      dtype=float32)
Coordinates:
  * x            (x) float64 4.361e+05 4.361e+05 ... 4.579e+05 4.579e+05
  * y            (y) float64 4.525e+06 4.525e+06 ... 4.514e+06 4.514e+06
    spatial_ref  int64 0
  * band_name    (band_name) <U3 'B1' 'B10' 'B11' 'B2' 'B3' 'B4' 'B5' 'B6' 'B7'
  * time         (time) datetime64[ns] 2021-09-12 2021-09-28 ... 2018-12-25
Attributes:
    _FillValue:    65535.0
    scale_factor:  0.0001
    add_offset:    0.0

With this output I'm be able to compute many spectral indexes, for examples NDVI.

My question: is this a right approach?
I think is better to use a Dataset for every scenes instead of a DataArray. I'm a little bit confused...
What is your approach in this case?

[scottyhq]

@MaxDragonheart it's a matter of preference to use the DataArray or DataSet representation with bands as variables. Xarray allows you to move back and forth with to_dataset() and to_array() methods.

Just be careful that the CRS (spatial_ref) of your DataArrays match as Xarray will not automatically reproject data in different CRS when constructing a cube. You might check out this example notebook for a DataSet representation https://planetarycomputer.microsoft.com/dataset/landsat-c2-l2#Example-Notebook:

There are so many ways to get Landsat data these days, where are you obtaining it from? Libraries that build on Xarray like stackstac or odc-stac make constructing DataSets easy and will even handle reprojection if the Landsat data has accompanying STAC metadata (like the planetary computer example above).

If you don't have the STAC metadata, your approach works fine. But if the CRS of scenes is different you need to manually deal with matching the grids of different data arrays, which the rioxarray extension can do

[MaxDragonheart]

Thank you for the answer @scottyhq :)
I've putted all inside one Dataset with which I've easely computed NDVI and other indexs useful for me.

Unfortunately I haven't a STAC catalogue, my data comes from earthexplorer but thanks to you I found out Miscrosoft Planetary Computer Program. I send the registration form and I hope that they accept my submission :)