来自特定地理位置的 CMEMS netcdf 文件的温度和盐度

Question

我想获取特定地理位置的海面（仅顶层）的温度['thetao']和盐度['so']。

我在此网站上找到了如何执行此操作的指南。

import netCDF4 as nc
import numpy as np

fn = "\\...\...\Downloads\global-analysis-forecast-phy-001-024_1647367066622.nc"  # path to netcdf file

ds = nc.Dataset(fn)  # read as netcdf dataset

print(ds)
print(ds.variables.keys()) # get all variable names

temp = ds.variables['thetao']
sal = ds.variables['so']
lat,lon = ds.variables['latitude'], ds.variables['longitude']

# extract lat/lon values (in degrees) to numpy arrays
latvals = lat[:]; lonvals = lon[:]

# a function to find the index of the point closest pt
# (in squared distance) to give lat/lon value.
def getclosest_ij(lats,lons,latpt,lonpt):
  # find squared distance of every point on grid
  dist_sq = (lats-latpt)**2 + (lons-lonpt)**2
  # 1D index of minimum dist_sq element
  minindex_flattened = dist_sq.argmin()
  # Get 2D index for latvals and lonvals arrays from 1D index
  return np.unravel_index(minindex_flattened, lats.shape)

iy_min, ix_min = getclosest_ij(latvals, lonvals, 50., -140)
print(iy_min)
print(ix_min)


# Read values out of the netCDF file for temperature and salinity
print('%7.4f %s' % (temp[0,0,iy_min,ix_min], temp.units))
print('%7.4f %s' % (sal[0,0,iy_min,ix_min], sal.units))

我正在使用的 nc 文件的一些详细信息：

dimensions(sizes): time(1), depth(1), latitude(2041), longitude(4320)
 variables(dimensions): float32 depth(depth), float32 latitude(latitude), int16 thetao(time, depth, latitude, longitude), float32 time(time), int16 so(time, depth, latitude, longitude), float32 longitude(longitude)
    groups:
dict_keys(['depth', 'latitude', 'thetao', 'time', 'so', 'longitude'])

我收到此错误：

dist_sq = (lats-latpt)**2 + (lons-lonpt)**2
ValueError: operands could not be broadcast together with shapes (2041,) (4320,)

我怀疑形状/数组存在问题。在网站的示例中（上面的链接），纬度和经度具有 (x,y)，但是此 NC 文件仅具有纬度 (2041,) 和经度 (4320,)。

我该如何解决这个问题？

Answer 1

这是因为纬度和经度是大小不同的向量...

如果使用 WGS84 或度作为单位，我通常会这样做：

lonm,latm = np.meshgrid(lons,lats)
dmat = (np.cos(latm*np.pi/180.0)*(lonm-lonpt)*60.*1852)**2+((latm-latpt)*60.*1852)**2

现在您可以找到最近的点：

kkd = np.where(dmat==np.nanmin(dmat))