Sentinel2 Scene Classification Map (SCM) and Quality Indicators (QI) band stacking

[goriliukasbuxton]

Is it posible to stack 'eo:bands' with data quality bands into one dask stack?
Would like to do cloud masking by these quality bands

Thank you,

[scottyhq]

Hi @goriliukasbuxton , short answer is yes, but right now you have to combine these different assets yourself due to the different ground-sample-distance (10m versus 20m postings). Currently working on ways to make this easier (#75),

but for now have a look at this example https://nbviewer.jupyter.org/gist/scottyhq/3bda794762139729ace8db457dac0248 which

1. interpolates SCL to match higher resolution data
2. adds SCL as a coordinate mask in another xarray dataset (e.g. b4 or visual S2 asset)
3. use xarray.where() function to mask pixels based on value of SCL

[goriliukasbuxton]

@scottyhq , very nice!
trying to stack all the items:

def scl_item_to_dataset(item):
    da = catalog[item.id].SCL.to_dask()
    da['band'] = ['scl'] 
    da = da.expand_dims(time=[pd.to_datetime(item.datetime)])
    ds = da.to_dataset(dim='band')
    return ds
    

scl_list = []
for i,item in gf.iterrows():
    #print(item)
    stac_scl = client.submit(scl_item_to_dataset,item)

    scl_list.append(stac_scl)
results_scl = client.gather(scl_list)
results_scl

######################
[<xarray.Dataset>
Dimensions: (time: 1, x: 5490, y: 5490)
Coordinates:

• time (time) datetime64[ns] 2020-08-28T16:03:30
• y (y) float64 4e+06 4e+06 4e+06 4e+06 ... 3.89e+06 3.89e+06 3.89e+06
• x (x) float64 2e+05 2e+05 2e+05 ... 3.097e+05 3.098e+05 3.098e+05
  Data variables:
  scl (time, y, x) uint8 dask.array<chunksize=(1, 5490, 5490), meta=np.ndarray>
  Attributes:
  transform: (20.0, 0.0, 199980.0, 0.0, -20.0, 4000020.0)
  crs: +init=epsg:32618
  res: (20.0, 20.0)
  is_tiled: 1
  nodatavals: (0.0,)
  scales: (1.0,)
  offsets: (0.0,)
  AREA_OR_POINT: Area
  OVR_RESAMPLING_ALG: MODE,

also tryinfg to szueeze the stack but its not working:

DS_SCL = xr.concat(results_scl, dim='time')
print('Dataset size: [Gb]', DS_SCL.nbytes/1e9)

DS_SCL2 = DS_SCL.interp(y=DS['y'], x=DS['x'], # from vidual Dataset below
                            #kwargs={"fill_value": 12}, #np.nan by default could use 0 or 12
                            kwargs={'fill_value':'extrapolate'},
                            method='nearest').squeeze('band', drop=False)
DS_SCL2

#######################

KeyError Traceback (most recent call last)
in
5 #kwargs={"fill_value": 12}, #np.nan by default could use 0 or 12
6 kwargs={'fill_value':'extrapolate'},
----> 7 method='nearest').squeeze('band', drop=False)
8 DS_SCL2

~\Anaconda3\envs\rasterio_test\lib\site-packages\xarray\core\common.py in squeeze(self, dim, drop, axis)
365 If drop=True, drop squeezed coordinates instead of making them
366 scalar.
--> 367 axis : None or int or iterable of int, optional
368 Like dim, but positional.
369

~\Anaconda3\envs\rasterio_test\lib\site-packages\xarray\core\common.py in get_squeeze_dims(xarray_obj, dim, axis)
320 dim = list(dim)
321 elif dim is not None:
--> 322 dim = [dim]
323 else:
324 assert axis is not None

~\Anaconda3\envs\rasterio_test\lib\site-packages\xarray\core\common.py in (.0)
320 dim = list(dim)
321 elif dim is not None:
--> 322 dim = [dim]
323 else:
324 assert axis is not None

~\Anaconda3\envs\rasterio_test\lib\site-packages\xarray\core\utils.py in getitem(self, key)
424
425 class Frozen(Mapping[K, V]):
--> 426 """Wrapper around an object implementing the mapping interface to make it
427 immutable. If you really want to modify the mapping, the mutable version is
428 saved under the mapping attribute.

~\Anaconda3\envs\rasterio_test\lib\site-packages\xarray\core\utils.py in getitem(self, key)
455
456 class HybridMappingProxy(Mapping[K, V]):
--> 457 """Implements the Mapping interface. Uses the wrapped mapping for item lookup
458 and a separate wrapped keys collection for iteration.
459

KeyError: 'band'

[RichardScottOZ]

No band in the above, only scl?