Python空间相关性计算报错：维度0数组大小不匹配致拼接错误

数据集信息 #

海冰数据 #

• 维度：LON(360)、LAT(173)、TIME(444)、bnds(2)
• 坐标：
  • LON：-179.5 ~ 179.5（步长1）
  • LAT：-82.5 ~ 89.5（步长1）
  • TIME：1985-01-01 ~ 2021-12-01
• 变量：SIC(TIME, LAT, LON)

碳数据 #

• 维度：time(444)、latitude(173)、longitude(360)
• 坐标：
  • time：1985-01-15 ~ 2021-12-15
  • latitude：-82.5 ~ 89.5（步长1）
  • longitude：0.5 ~ 359.5（步长1）
• 变量：spco2(time, latitude, longitude)等

现有代码 #

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from scipy.interpolate import griddata

# Load sea ice data
sea_ice_data = xr.open_dataset(path_ice)
sea_ice_sic = sea_ice_data['SIC']  

# Load carbon data
carbon_data = xr.open_dataset(path_carbon)
carbon_spco2 = carbon_data['spco2']  

# Convert carbon longitudes to -180 to +180 range
carbon_spco2['longitude'] = (carbon_spco2['longitude'] + 180) % 360 - 180

# Flatten latitude and longitude coordinates
carbon_coords = np.column_stack((carbon_spco2['latitude'].values.flatten(), carbon_spco2['longitude'].values.flatten()))

# Calculate spatial correlation map
correlation_map = np.empty((len(sea_ice_sic['LAT']), len(sea_ice_sic['LON'])))
for lat_idx, lat in enumerate(sea_ice_sic['LAT']):
    for lon_idx, lon in enumerate(sea_ice_sic['LON']):
        sic_values = sea_ice_sic.sel(LAT=lat, LON=lon, method='nearest').values

        # Interpolate carbon data to sea ice grid
        spco2_values = griddata(
            carbon_coords,
            carbon_spco2.values.flatten(),
            (lat, lon),
            method='nearest'
        )
        
        correlation_map[lat_idx, lon_idx] = np.corrcoef(sic_values, spco2_values)[0, 1]

# Create a Cartopy projection
projection = ccrs.PlateCarree()

# Plot the spatial correlation map using Cartopy
plt.figure(figsize=(12, 8))
ax = plt.axes(projection=projection)
ax.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())
ax.coastlines()

# Plot the correlation map as an image
plt.imshow(correlation_map, cmap='RdBu_r', vmin=-1, vmax=1, extent=(-180, 180, -90, 90), origin='upper', transform=ccrs.PlateCarree())

# Add a colorbar
cbar = plt.colorbar(label='Correlation Coefficient', orientation='vertical', shrink=0.7)
cbar.ax.tick_params(labelsize=10)

plt.title('Spatial Correlation between Sea Ice and Carbon')
plt.show()

报错信息 #

ValueError                                Traceback (most recent call last)
<ipython-input-8-1d668216ff21> in <cell line: 20>()
     18 
     19 # Flatten latitude and longitude coordinates
---> 20 carbon_coords = np.column_stack((carbon_spco2['latitude'].values.flatten(), carbon_spco2['longitude'].values.flatten()))
     21 
     22 # Calculate spatial correlation map

2 frames
/usr/local/lib/python3.10/dist-packages/numpy/core/overrides.py in column_stack(*args, **kwargs)

/usr/local/lib/python3.10/dist-packages/numpy/lib/shape_base.py in column_stack(tup)
    654             arr = array(arr, copy=False, subok=True, ndmin=2).T
    655         arrays.append(arr)
---> 656     return _nx.concatenate(arrays, 1)
    657 
    658 

/usr/local/lib/python3.10/dist-packages/numpy/core/overrides.py in concatenate(*args, **kwargs)

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 173 and the array at index 1 has size 360

问题分析与解决方案 #

错误原因 #

carbon_spco2['latitude']是173个元素的一维数组，longitude是360个元素的一维数组，直接flatten后拼接会导致维度不匹配。我们需要的是每个网格点的(lat, lon)坐标对（共173×360=62280个），而非单独维度的坐标。

修正步骤 #

1. 正确转换并排序碳数据经度：转换后经度会变成-179.5到179.5，但顺序会混乱，需要重新排序保证经度从西到东递增。
2. 对齐时间维度：海冰数据时间是每月1日，碳数据是每月15日，用线性插值对齐时间。
3. 用xarray内置重网格工具替代循环：xarray.DataArray.interp可以直接将碳数据插值到海冰的网格上，比手动循环+griddata高效且简洁。
4. 计算空间相关性：对每个网格点的时间序列计算相关系数，同时排除NaN值避免报错。

修正后代码 #

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# 加载数据
sea_ice_data = xr.open_dataset(path_ice)
sic = sea_ice_data['SIC']
carbon_data = xr.open_dataset(path_carbon)
spco2 = carbon_data['spco2']

# 1. 转换碳数据经度到-180~180范围并排序
spco2 = spco2.assign_coords(longitude=((spco2.longitude + 180) % 360) - 180)
spco2 = spco2.sortby('longitude')

# 2. 对齐时间维度：将碳数据插值到海冰数据的时间点
spco2_aligned = spco2.interp(time=sic.TIME)

# 3. 将碳数据插值到海冰的空间网格
spco2_regridded = spco2_aligned.interp(LAT=spco2_aligned.latitude, LON=spco2_aligned.longitude).rename(
    {'latitude': 'LAT', 'longitude': 'LON'}
)

# 4. 计算每个网格点的时间序列相关系数
def corr_coeff(x, y):
    mask = ~np.isnan(x) & ~np.isnan(y)
    if np.sum(mask) < 2:
        return np.nan
    return np.corrcoef(x[mask], y[mask])[0, 1]

correlation_map = xr.apply_ufunc(
    corr_coeff, sic, spco2_regridded,
    input_core_dims=[['TIME'], ['time']],
    vectorize=True,
    output_dtypes=[float]
)

# 绘图
plt.figure(figsize=(12, 8))
projection = ccrs.PlateCarree()
ax = plt.axes(projection=projection)
ax.set_extent([-180, 180, -90, 90], crs=projection)
ax.add_feature(cfeature.COASTLINE)

correlation_map.plot(
    ax=ax,
    cmap='RdBu_r',
    vmin=-1,
    vmax=1,
    transform=projection,
    cbar_kwargs={'label': 'Correlation Coefficient', 'shrink': 0.7}
)

plt.title('Spatial Correlation between Sea Ice Concentration (SIC) and Surface pCO2 (spco2)')
plt.show()

额外说明 #

• 若数据集较大，interp可能较慢，可以考虑用xesmf库进行更高效的重网格化（适合规则/不规则网格转换）。
• 时间对齐时，若不想插值，也可以用sel(method='nearest')选择最近的时间点，但插值能更准确匹配月尺度数据。

内容的提问来源于stack exchange，提问作者asheef_ik