使用Planetary Computer调用OTCI数据时遇xarray读取报错求助

错误信息 #

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

原问题代码 #

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. 

解决方案 #

错误根源在于Sentinel-3 OLCI L2的NetCDF资产并非常规单波段栅格，StackSTAC的默认栅格处理逻辑无法适配其结构。改用直接通过xarray读取签名后的NetCDF文件的方式，可解决该问题：

修正后的代码 #

import pystac_client
import planetary_computer
import numpy as np
import xarray as xr
import rasterio
import warnings

# 忽略地理参考警告（NetCDF的地理信息存储方式与常规栅格不同）
warnings.filterwarnings("ignore", category=rasterio.errors.NotGeoreferencedWarning)

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)
min_lon, min_lat, max_lon, max_lat = bbox

# 初始化STAC客户端
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")
items = list(search.items())

# 遍历每个item，提取OTCI数据并裁剪到感兴趣区域
otci_list = []
for item in items:
    # 获取签名后的OTCI资产URL
    otci_url = item.assets["otci"].href
    # 打开NetCDF文件
    ds = xr.open_dataset(otci_url)
    # 裁剪到目标区域
    ds_cropped = ds.where(
        (ds.lon >= min_lon) & (ds.lon <= max_lon) & 
        (ds.lat >= min_lat) & (ds.lat <= max_lat),
        drop=True
    )
    # 添加时间维度
    ds_cropped = ds_cropped.assign_coords(time=item.datetime)
    otci_list.append(ds_cropped["otci"])

# 合并所有时间步的数据
otci_timeseries = xr.concat(otci_list, dim="time")

# 提取指定像素的时间序列并转为NumPy数组
# 注意：这里需要根据实际裁剪后的坐标选择索引，示例用[0,0]
pixel_timeseries = otci_timeseries[:, 0, 0].values
print(pixel_timeseries)

关键说明 #

1. 直接读取NetCDF：跳过StackSTAC，用xarray直接打开签名后的OTCI资产URL，适配Sentinel-3的NetCDF结构。
2. 区域裁剪：利用xarray的where方法根据经纬度范围裁剪数据，保留目标区域。
3. 时间维度合并：为每个时间步的数据集添加time坐标，再通过xr.concat合并为完整的时间序列。
4. 忽略无关警告：NetCDF的地理参考方式与常规栅格不同，相关警告不影响数据处理，可直接忽略。

内容的提问来源于stack exchange，提问作者Dávid D.Kovács