Hello,

I've been trying to use Cubo to download cubes from the Sentinel 5p collection available on the planetary computer:

https://planetarycomputer.microsoft.com/dataset/sentinel-5p-l2-netcdf

The download process works quite well by simply adapting the code snippet you provide in the documentation. But whenever I then try to apply the to_numpy() method to the resulting DataArray i get the rasterio ValueError "Can't read closed raster file". This doesn't happen when attempting to do the same with the sentinel-2 l2a collection.

Here's the code snippet I'm using:

import cubo
import xarray as xr

Not working

da = cubo.create(
lat=49.80,
lon=9.42,
collection="sentinel-5p-l2-netcdf",
bands=["no2"],
start_date="2021-01-01",
end_date="2021-01-01",
)

da.to_numpy()

Thank you in advance and best regards!

I'd like to be able to use something like EPSG:4326 to select from my xarray. Reprojection support would be a valuable feature to me. I can understand how the focus of this library should probably just be the creation of Xarray cubes.

The saved coord cubo:distance_from_center prevents the data to be reprojected.

data_sentinel2_subset = data_sentinel2.isel(time=0)
data_sentinel2_subset = data_sentinel2_subset.drop_vars(['cubo:distance_from_center']) # needed line
data_sentinel2_subset_reprojected = data_sentinel2_subset.rio.reproject('EPSG:3035')

Error:

ValueError: IndexVariable objects must be 1-dimensional

If you do not have a project registered in gcloud, it is better to remove this and leave it outside, I think the initialization should depend on the user

Is there a way to get a batched cubo?

da3 = cubo.create(
lat=45.30493249820944,
lon=132.6933519142134,
collection="sentinel-2-l2a",
bands=["B02","B03","B04"],
start_date="2020-01-01",
end_date="2021-01-01",
edge_size=64,
resolution=10,
query={"eo:cloud_cover": {"lt": 10}}
)
da3
(da3.sel(band=["B04","B03","B02"])/2000).clip(0,1).plot.imshow(col="time",col_wrap = 5)

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

ValueError Traceback (most recent call last)
Cell In[20], line 1
----> 1 (da3.sel(band=["B04","B03","B02"])/5000).clip(0,1).plot.imshow(col="time",col_wrap = 5)

File /srv/conda/envs/notebook/lib/python3.11/site-packages/xarray/plot/accessor.py:430, in DataArrayPlotAccessor.imshow(self, *args, **kwargs)
428 @functools.wraps(dataarray_plot.imshow)
429 def imshow(self, *args, **kwargs) -> AxesImage:
--> 430 return dataarray_plot.imshow(self._da, *args, **kwargs)

File /srv/conda/envs/notebook/lib/python3.11/site-packages/xarray/plot/dataarray_plot.py:1495, in _plot2d..newplotfunc(failed resolving arguments)
1493 # Need the decorated plotting function
1494 allargs["plotfunc"] = globals()[plotfunc.name]
-> 1495 return _easy_facetgrid(darray, kind="dataarray", **allargs)
1497 if darray.ndim == 0 or darray.size == 0:
1498 # TypeError to be consistent with pandas
1499 raise TypeError("No numeric data to plot.")

File /srv/conda/envs/notebook/lib/python3.11/site-packages/xarray/plot/facetgrid.py:1042, in _easy_facetgrid(data, plotfunc, kind, x, y, row, col, col_wrap, sharex, sharey, aspect, size, subplot_kws, ax, figsize, **kwargs)
1039 sharex = False
1040 sharey = False
-> 1042 g = FacetGrid(
1043 data=data,
1044 col=col,
1045 row=row,
1046 col_wrap=col_wrap,
1047 sharex=sharex,
1048 sharey=sharey,
1049 figsize=figsize,
1050 aspect=aspect,
1051 size=size,
1052 subplot_kws=subplot_kws,
1053 )
1055 if kind == "line":
1056 return g.map_dataarray_line(plotfunc, x, y, **kwargs)

File /srv/conda/envs/notebook/lib/python3.11/site-packages/xarray/plot/facetgrid.py:175, in FacetGrid.init(self, data, col, row, col_wrap, sharex, sharey, figsize, aspect, size, subplot_kws)
173 rep_row = row is not None and not data[row].to_index().is_unique
174 if rep_col or rep_row:
--> 175 raise ValueError(
176 "Coordinates used for faceting cannot "
177 "contain repeated (nonunique) values."
178 )
180 # single_group is the grouping variable, if there is exactly one
181 single_group: bool | Hashable

ValueError: Coordinates used for faceting cannot contain repeated (nonunique) values.

Running the example on the README.md produces the following error

AssertionError: `out_epsg` not found. out_epsg='EPSG:32618'

The code that produces this error:

import cubo
import xarray as xr  

da = cubo.create(
    lat=4.31, # Central latitude of the cube
    lon=-76.2, # Central longitude of the cube
    collection="sentinel-2-l2a", # Name of the STAC collection
    bands=["B02","B03","B04"], # Bands to retrieve
    start_date="2021-06-01", # Start date of the cube
    end_date="2021-06-10", # End date of the cube
    edge_size=64, # Edge size of the cube (px)
    resolution=10, # Pixel size of the cube (m)
)

Thank you for developing this! Given that it is one function and the main dependency is stackstac, integrating it into stackstac would mean that I can just mamba install stackstac instead of additionally pip installing cubo (and potentially mixing pypi and conda dependencies). Happy to work on this! It would make it possible for us to then use cubo in https://carpentries-incubator.github.io/geospatial-python/ without adding another dependency.

Hello Dave,

I hope this message finds you well. I recently started working with Cubo to develop a pipeline that generates cloudless mosaics from Sentinel-2 satellite data. While working on this project, I encountered a few issues that I believe would be useful to share with you. Additionally, I am available to assist you in further developing this package if you are interested.

In higher latitudes, using latitude and longitude as the center of patches can cause issues because the points are not on the same line. In other cases, the generated patches at the borders contain NaN values.

Kind regards,
Behzad

Hello,
Trying the basic example:

da = cubo.create(
    lat=50,
    lon=10,
    collection="sentinel-2-l2a",
    bands=["B02","B03","B04"],
    start_date="2021-06-01",
    end_date="2021-06-10",
    edge_size=32,
    resolution=10,
)

I have this error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/bdeneu/miniconda3/envs/glc/lib/python3.9/site-packages/cubo/cubo.py", line 111, in create
    cube = stackstac.stack(
  File "/home/bdeneu/miniconda3/envs/glc/lib/python3.9/site-packages/stackstac/stack.py", line 287, in stack
    asset_table, spec, asset_ids, plain_items = prepare_items(
  File "/home/bdeneu/miniconda3/envs/glc/lib/python3.9/site-packages/stackstac/prepare.py", line 166, in prepare_items
    assert isinstance(out_epsg, int), f"`out_epsg` not found. {out_epsg=}"
AssertionError: `out_epsg` not found. out_epsg='EPSG:32632'

It seems that stackstac out_epsg expect an int and received a string !
My version of python is 3.9.16 and stackstac 0.4.4.
Best regards

The extent of a cubo.create operation is defined using a centerpoint coordinate (lat, lon), and a radius (edge_size).
This is a feature request for a GEOJSON based load. Where either a GEOJSON bounding box or polygon is used to define the area extent of the load.