xarray Can't read closed raster file ValueError by using Planetary Computer

51 Views Asked by At

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 :)

0

There are 0 best solutions below