Python如何反转netCDF数组解决高程数据上下颠倒问题

问题说明 #

从netCDF格式高程文件中提取CSV坐标点对应高程值时，发现文件存储的高程数据南北方向上下颠倒，直接提取的结果位置完全错位。此前尝试将CSV读取的纬度值乘以-1的方案会篡改原始坐标，不符合最终输出的坐标准确性要求。

原有问题代码 #

import numpy as np
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Path, PathPatch

csv_data = np.loadtxt('CSV with target coordinates',skiprows=1,delimiter=',')
num_el = csv_data[:,0]
lat = csv_data[:,1]
lon = csv_data[:,2]
value = csv_data[:,3]

data = Dataset("elevation Data",'r')
lon_range = data.variables['x_range'][:]
lat_range = data.variables['y_range'][:]
topo_range = data.variables['z_range'][:]
spacing = data.variables['spacing'][:]
dimension = data.variables['dimension'][:]
z = data.variables['z'][:]
lon_num =  dimension[0]
lat_num =  dimension[1]

etopo_lon = np.linspace(lon_range[0],lon_range[1],dimension[0])
etopo_lat = np.linspace(lat_range[0],lat_range[1],dimension[1])
topo = np.reshape(z, (lat_num, lon_num))

height = np.empty_like(num_el)
desired_lat_idx = np.empty_like(num_el)
desired_lon_idx = np.empty_like(num_el)
for i in range(len(num_el)): 
    tmp_lat = np.abs(etopo_lat - lat[i]).argmin()
    tmp_lon = np.abs(etopo_lon - lon[i]).argmin()
    desired_lat_idx[i] = tmp_lat
    desired_lon_idx[i] = tmp_lon
    height[i] = topo[tmp_lat,tmp_lon]
height[height<-10]=0 

print(len(desired_lat_idx))
print(len(desired_lon_idx))
print(len(height))


dfl= pd.DataFrame({
    'Latitude' : lat.reshape(-1),
    'Longitude': lon.reshape(-1),
    'Altitude': height.reshape(-1)
});
print(dfl)

# 要求此处纬度值必须和原始CSV一致，不能修改
df =dfl
lat=np.array(df['Latitude'])
lon=np.array(df['Longitude'])
val=np.array(df['Altitude'])

m = Basemap(projection='robin', lon_0=0, lat_0=0, resolution='l',area_thresh=1000)
m.drawcoastlines(color = 'black')
x,y = m(lon,lat)
colormesh= m.contourf(x,y,val,100, tri=True,  cmap = 'terrain')
plt.colorbar(location='bottom',pad=0.04,fraction=0.06)
plt.show()

无效方案 #

直接反转原始CSV纬度值的方案会导致坐标准确性失效，不可用：

lat = csv_data[:,1]
lat= lat*(-1)

根因分析 #

绝大多数公开地形netCDF文件（比如ETOPO系列）的z数组存储顺序为纬度从北向南排列（即第一行对应最北边界，最后一行对应最南边界），但代码中用np.linspace生成的etopo_lat是按纬度值从小到大（从南向北）排列，两者顺序不匹配，导致按纬度索引取高程时取到了南北对称位置的错误值。

修正方法 #

不需要修改任何原始CSV的坐标值，只需要在reshape得到topo数组后，沿纬度轴（第0轴）上下翻转数组，让topo的行顺序和etopo_lat的纬度顺序一一对应即可。
找到代码中这一行：

topo = np.reshape(z, (lat_num, lon_num))

在它后面新增一行翻转操作：

topo = np.flipud(topo)

其余代码完全不需要改动，后续的索引匹配、DataFrame生成、绘图逻辑都可以正常运行，最终输出的纬度字段完全保留原始CSV的准确值，提取的高程值也会和坐标正确匹配。

验证提示：可以先选取已知高程的点位（如青藏高原点位高程应在4000m以上、沿海平原点位高程在0-100m区间）做校验，确认提取结果正确。

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