Sentinel3 OLCI (chl) Python 上 netcdf 文件的平均值

Question

我在尝试获得 Sentinel 3 图像的月平均值时遇到了一些麻烦……真的。 Python，Matlab，我们两个人都陷入了这个问题。

主要原因在于这些图像的信息不在单个 netcdf 文件中，而是整齐地与坐标和产品一起放置。相反，它们都位于一天文件夹中的单独文件中 different .nc files with different information each, about one single satellite image。据我了解，SNAP 使用 xmlxs 文件来处理所有这些单独的 .nc 文件。

现在，我认为尝试合并和创建/编辑 .nc 文件以创建一个新的每日 .nc 文件是一个好主意，其中包括叶绿素、坐标，还不如添加时间。稍后，我将合并这些新的，以便能够使用 xarray 计算每月平均值。至少那是我的想法，但我不能做第一部分。这可能是一个明显的解决方案，但这是我尝试过的，使用 xarray 模块

import os
import numpy as np
import xarray as xr
import netCDF4
from netCDF4 import Dataset
  
nc_folder = df_try.iloc[0] #folder where the image files are

#open dataset in xarray
nc_chl = xr.open_dataset(str(nc_folder['path']) + '/' + 'chl_nn.nc') #path to chlorophyll file 
nc_chl

n_coord =xr.open_dataset(str(nc_folder['path'])+ '/'+ 'geo_coordinates.nc') #path to coordinates file

n_time = xr.open_dataset(str(nc_folder['path'])+ '/' + 'time_coordinates.nc') #path to time file

ds_grid = [[nc_chl], [n_coord], [n_time]]

combined = xr.combine_nested(ds_grid, concat_dim=[None, None]) 
combined #dataset with all but not recognizing coordinates

ds = combined.rename({'latitude': 'lat', 'longitude': 'lon', 'time_stamp' : 'time'}).set_coords(['lon', 'lat', 'time']) #dataset recognizing coordinates as coordinates
ds 

它给出了一个数据集

维度：列 4865 行：4091

3 个坐标（纬度、经度和时间）和 chl 变量。

现在，它不会保存到 netcdf4（我尝试过，但出现错误），但我也在想是否有人知道另一种平均方法？我有三年的图像（从 2017 年开始到 2019 年结束），我需要以不同的方式（每月、季节性......）进行平均。我目前的主要问题是叶绿素值与地理坐标是分开的，因此仅直接使用叶绿素文件是行不通的，只会弄得一团糟。

有什么建议吗？

Answer 1

这里有两个选项：

使用 xarray

​​>

在 xarray 中，您可以将它们添加为 坐标。这有点棘手，因为geo_coordinates.nc 文件中的坐标也是多维的。

可能的解决方案如下：

import netCDF4 
import xarray as xr
import matplotlib.pyplot as plt

# paths
root = r'C:<your_path>\S3B_OL_2_WFR____20201015.SEN3\chl_nn.nc'        #set path to chl file
coor = r'C:<your_path>\S3B_OL_2_WFR____20201015.SEN3\geo_coordinates.nc' #set path to the coordinates file

# loading xarray datasets
ds = xr.open_dataset(root)
olci_geo_coords = xr.open_dataset(coor)

# extracting coordinates
lat = olci_geo_coords.latitude.data
lon = olci_geo_coords.longitude.data

# assign coordinates to the chl dataset (needs to refer to both the dimensions of our dataset)
ds = ds.assign_coords({"lon":(["rows","columns"], lon), "lat":(["rows","columns"], lat)})

# clip the image (add your own coordinates)
area_of_interest = ds.where((10 < ds.lon) & (ds.lon < 12) & (58 < ds.lat) & (ds.lat < 59), drop=True)

# simple plot with coordinates as axis
plt.figure(figsize=(15,15))
area_of_interest["CHL_NN"].plot(x="lon",y="lat")

更简单的方法是将它们作为变量添加到新数据集中：

# path to the folder
root = r'C:<your_path>\S3B_OL_2_WFR____20201015.SEN3\*.nc'        #set path to chl file

# create a dataset by combining nc files (coordinates will become variables)
ds = xr.open_mfdataset(root,combine = 'by_coords')

但是在这种情况下，当您绘制图像或对其进行剪辑时，您不能直接使用坐标。

使用 snappy

​​>

在 python 中，snappy 包可用并且基于 SNAP 工具箱（在 JAVA 上实现）。检查：https://senbox.atlassian.net/wiki/spaces/SNAP/pages/19300362/How+to+use+the+SNAP+API+from+Python

安装后（不幸的是，snappy 仅支持 python 2.7、3.3 或 3.4），您可以直接在 python 上使用可用的 SNAP 函数来聚合您的卫星图像并创建周/月平均值。然后，您无需合并 lon、lat netcdf 文件，因为您将在 xfdumanifest.xml 上工作，SNAP 会处理这些。

这是一个例子。它也执行聚合（根据两个 chl nc 文件计算平均值）：

from snappy import ProductIO, WKTReader
from snappy import jpy
from snappy import GPF
from snappy import HashMap

# setting the aggregator method
aggregator_average_config = snappy.jpy.get_type('org.esa.snap.binning.aggregators.AggregatorAverage$Config')
agg_avg_chl = aggregator_average_config('CHL_NN')

# creating the hashmap to store the parameters
HashMap = snappy.jpy.get_type('java.util.HashMap')
parameters = HashMap()

#creating the aggregator array
aggregators = snappy.jpy.array('org.esa.snap.binning.aggregators.AggregatorAverage$Config', 1)
#adding my aggregators in the list
aggregators[0] = agg_avg_chl

# set parameters
# output directory 
dir_out = 'level-3_py_dynamic.dim'
parameters.put('outputFile', dir_out)

# number of rows (directly linked with resolution)
parameters.put('numRows', 66792) # to have about 300 meters spatial resolution

# aggregators list
parameters.put('aggregators', aggregators)

# Region to clip the aggregation on
wkt="POLYGON ((8.923302175377243 59.55648108694149, 13.488748662344074 59.11388968719029,12.480488185001589 56.690625338725155, 8.212366327767503 57.12425256476263,8.923302175377243 59.55648108694149))"
geom = WKTReader().read(wkt)
parameters.put('region', geom)

# Source product path 
path_15 = r"C:<your_path>\S3B_OL_2_WFR____20201015.SEN3\xfdumanifest.xml"
path_16 = r"C:\<your_path>\S3B_OL_2_WFR____20201016.SEN3\xfdumanifest.xml"
path = path_15 + "," + path_16
parameters.put('sourceProductPaths', path)

#result = snappy.GPF.createProduct('Binning', parameters, (source_p1, source_p2))

# create results
result = snappy.GPF.createProduct('Binning', parameters) #to be used with product paths specified in the parameters hashmap

print("results stored in: {0}".format(dir_out) )

我很新，对这个话题很感兴趣，很高兴听到你/其他的解决方案！