I would like to use the OTCI dataset from the Planetary Computer Sentinel-3 Land dataset, by accessning it with the "stacstack" api, convert it to an xarray DataArray and do some further calculations. However, when I try to access the "otci" dataset as an array, the error I get is the following:
usr/local/lib/python3.10/dist-packages/stackstac/rio_reader.py:327: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
ds = SelfCleaningDatasetReader(self.url, sharing=False)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-4-5e459c4b60f6> in <cell line: 10>()
8 otci= xr_stack.squeeze() # shape of dataarray: (time, band,y,x) now converted to (time,y,x) by using .squeeze() and removing the "band" dimension
9 timeseries = otci[:,1,1] # Now we access all the timestamps, and try to get the 1x1 pixel
---> 10 timeseries.values #we want to convert to numpy to plot it later on, but here it fails. The error which it gave me: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.ds = SelfCleaningDatasetReader(self.url, sharing=False)
11
12
10 frames
/usr/local/lib/python3.10/dist-packages/stackstac/rio_reader.py in _open(self)
338 ds.close()
339 raise RuntimeError(
--> 340 f"Assets must have exactly 1 band, but file {self.url!r} has {ds.count}. "
341 "We can't currently handle multi-band rasters (each band has to be "
342 "a separate STAC asset), so you'll need to exclude this asset from your analysis."
rasterio/_base.pyx in rasterio._base.DatasetBase.count.__get__()
ValueError: Can't read closed raster file
I attach the example code where I get this issue. I use Colab, but same error was happening on other IDEs too.
import pystac_client
import stackstac
import matplotlib.pyplot as plt
import warnings
import planetary_computer
import numpy as np
import xarray as xr
import fsspec
import rioxarray as rxr
import rasterio
import odc
area_of_interest = {
"type": "Polygon",
"coordinates":
[[[10.488476562499986, 48.1828494522882],
[10.488476562499986, 47.785761881918084],
[11.191601562499986, 47.785761881918084],
[11.191601562499986, 48.1828494522882]]]
}
bbox = rasterio.features.bounds(area_of_interest)
sentinel_search_url = "https://planetarycomputer.microsoft.com/api/stac/v1"
catalog = pystac_client.Client.open(sentinel_search_url,
modifier=planetary_computer.sign_inplace,
)
search = catalog.search(
bbox=bbox,
collections=["sentinel-3-olci-lfr-l2-netcdf"],
datetime="2020-05-01/2020-06-01")
item = next(search.items())
search_collection = search.item_collection()
xr_stack = stackstac.stack(search_collection,
assets=["otci"],
resolution=0.01, # this is in degrees
epsg = 4326,
bounds = bbox,
)
otci= xr_stack.squeeze() # shape of dataarray: (time, band,y,x) now converted to (time,y,x) by using .squeeze() and removing the "band" dimension
timeseries = otci[:,1,1] # Now we access all the timestamps, and try to get the 1x1 pixel
timeseries.values #we want to convert to numpy to plot it later on, but here it fails.
I would like to have the data as a numpy array, or any format that works for furhter mathematical analysis. Thanks and have a great day :)