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

Question

我想得到特定地理位置的海面(仅顶层(的温度['thetao']和盐度['so']。

我在这个网站上找到了如何做到这一点的指导。

import netCDF4 as nc
import numpy as np
fn = "\......Downloadsglobal-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,)

我怀疑形状/阵列有问题。在网站的例子中(上面的链接(，Lat和Lon有一个(x，y(，但是这个NC文件只有Latitude(2041，(和Longitude(4320，(。

我该如何解决这个问题？

Answer 1

这是因为lats和lons是大小不同的向量。。。

如果使用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))