基于站点属性合并多份海洋观测NC文件的问题与解决

问题描述 #

尝试合并多份包含不同经纬度、不同深度的物理海洋学NC文件，使用xr.open_mfdataset实现，但合并后绘图仅显示单一结果，存在异常。

初始尝试代码 #

##Combining using concat_dim and nested method
ds = xr.open_mfdataset("33HQ20150809*.nc", concat_dim=['latitude'], combine= "nested")
ds.to_netcdf('geotraces2015_combined.nc')
df = xr.open_dataset("geotraces2015_combined.nc")

##Setting up values. Oxygen values are transposed so it matches same shape as lat and pressure. 
oxygen = df['oxygen'].values.transpose()
##Plotting using colourf
fig = plt.figure()
ax = fig.add_subplot(111)
plt.contourf(oxygen, cmap = 'inferno')
plt.gca().invert_yaxis()
cbar = plt.colorbar(label = 'Oxygen Concentration (umol kg-1')

NC文件结构说明 #

单个文件的xarray数据集结构如下（不同文件的pressure维度长度不同，经纬度、时间也各有差异）：

<xarray.Dataset>
Dimensions:         (pressure: 744, time: 1, latitude: 1, longitude: 1)
Coordinates:
  * pressure        (pressure) float64 0.0 1.0 2.0 3.0 ... 741.0 742.0 743.0
  * time            (time) datetime64[ns] 2015-08-12T18:13:00
  * latitude        (latitude) float32 60.25
  * longitude       (longitude) float32 -179.1
Data variables: (12/19)
    pressure_QC     (pressure) int16 ...
    temperature     (pressure) float64 ...
    temperature_QC  (pressure) int16 ...
    salinity        (pressure) float64 ...
    salinity_QC     (pressure) int16 ...
    oxygen          (pressure) float64 ...
    ...              ...
    CTDNOBS         (pressure) float64 ...
    CTDETIME        (pressure) float64 ...
    woce_date       (time) int32 ...
    woce_time       (time) int16 ...
    station         |S40 ...
    cast            |S40 ...
Attributes:
    EXPOCODE:                   33HQ20150809
    Conventions:                COARDS/WOCE
    WOCE_VERSION:               3.0
...

另一个示例文件结构：

<xarray.Dataset>
Dimensions:         (pressure: 179, time: 1, latitude: 1, longitude: 1)
Coordinates:
  * pressure        (pressure) float64 0.0 1.0 2.0 3.0 ... 176.0 177.0 178.0
  * time            (time) datetime64[ns] 2015-08-18T19:18:00
  * latitude        (latitude) float32 73.99
  * longitude       (longitude) float32 -168.8
Data variables: (12/19)
    pressure_QC     (pressure) int16 ...
    temperature     (pressure) float64 ...
    temperature_QC  (pressure) int16 ...
    salinity        (pressure) float64 ...
    salinity_QC     (pressure) int16 ...
    oxygen          (pressure) float64 ...
    ...              ...
    CTDNOBS         (pressure) float64 ...
    CTDETIME        (pressure) float64 ...
    woce_date       (time) int32 ...
    woce_time       (time) int16 ...
    station         |S40 ...
    cast            |S40 ...
Attributes:
    EXPOCODE:                   33HQ20150809
    Conventions:                COARDS/WOCE
    WOCE_VERSION:               3.0

中间尝试的方法 #

曾尝试用preprocess函数处理文件，进行坐标设置、维度压缩和扩展，但问题未解决：

def preprocess(ds):
    return ds.set_coords('station').squeeze(["latitude", "longitude", "time"]).expand_dims('station')
ds = xr.open_mfdataset('33HQ20150809*.nc', concat_dim='station', combine='nested', preprocess=preprocess)

最终解决方案 #

以station作为唯一合并维度，通过preprocess对每个文件执行维度压缩、坐标设置和扩展操作，最终实现正确合并与绘图：

import pandas as pd
import numpy as np
import os
import netCDF4 
import pathlib
import matplotlib.pyplot as plt

def preprocess(ds):
    return ds.set_coords('station').squeeze(["latitude", "longitude", "time"]).expand_dims('station')

ds = xr.open_mfdataset('filename*.nc', preprocess=preprocess, parallel=True)
ds = ds.sortby('latitude').transpose()

ds.oxygen.plot.contourf(x="latitude", y="pressure")
plt.gca().invert_yaxis()

内容的提问来源于stack exchange，提问作者Maria Cristina Alvarez