Sentinel3 OLCI (chl) Среднее значение файлов netcdf на Python

У меня проблемы с попыткой получить среднемесячное значение при включенных изображениях Sentinel 3 ... Все, правда. Python, Matlab, мы два человека, которые застряли в этой проблеме.

Основная причина связана с тем, что информация об этих изображениях не находится в одном файле netcdf, аккуратно помещенном с координатами и продуктами. Вместо этого все они находятся в отдельных файлах внутри однодневной папки в виде разных файлов .nc с разной информацией каждый , об одном спутниковом снимке. 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 года), которые мне нужно было бы усреднить разными способами (ежемесячно, сезонно ...). Моя основная текущая проблема заключается в том, что значения хлорофилла отделены от географических координат, поэтому прямое использование только файлов хлорофилла не должно работать и просто приведет к беспорядку.

Какие-либо предложения?


person Ananas    schedule 06.10.2020    source источник
comment
Лучшим решением было бы добавить координаты к файлам хлорофилла. Хотя без дополнительной информации я не могу предложить, как   -  person Robert Wilson    schedule 06.10.2020
comment
ну, я возился с xarray с comb_nested (согласно коду), пытаясь объединить и хлорофилл, и координаты в один файл .nc, сохраняя с помощью dataset.to_netcdf. Мне удалось его сохранить, но он выдает такое предупреждение? «SerializationWarning: сохранение переменной lat с данными с плавающей запятой как целочисленного типа dtype без какого-либо _FillValue для использования для NaN». Файл, похоже, работает для ежедневных изображений   -  person Ananas    schedule 07.10.2020
comment
cdo может быть проще: code.mpimet.mpg.de/projects/cdo/wiki/   -  person Robert Wilson    schedule 07.10.2020


Ответы (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')

Но в этом случае, когда вы рисуете изображение или обрезаете его, вы не можете напрямую использовать координаты.

Использование мгновенного

В 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), вы можете использовать доступную функцию SNAP непосредственно на python для агрегирования спутниковых изображений и создания средних значений за неделю / месяц. Тогда вам не нужно объединять файл 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) )

Я новичок и заинтересован в этой теме и буду рад услышать ваши / другие решения!

person cmotta    schedule 03.12.2020