为什么netcdf文件中的纬度和经度是二维数组?

wol*_*wol 2 python latitude-longitude coordinates netcdf python-xarray

我有 netCDF 文件,其中包含某个位置的温度数据。数据形状为 1450x900。

我正在我的应用程序中创建搜索功能,以查找具有纬度、经度值的温度数据。

因此,我从 netCDf 文件中提取了纬度和经度坐标数据,但我期望它们是一维数组,而是得到两个坐标形状均为 1450x900 的二维数组。

所以我的问题是:为什么它们是二维数组,而不是 1450 个纬度值和 900 个经度值?1450 纬度值和 900 经度值不是描述整个网格吗?

假设我们有 4x5 的正方形,用于定位网格最右边和最底部点的索引将为 [4, 5]。所以我的 x 索引为 [1, 2, 3, 4],y 索引为:[1, 2, 3, 4, 5]。总共 9 个索引足以定位该网格(由 20 个单元格组成)上的任何点。那么为什么 netcdf 文件中的 lat (x) 和 lon (y) 坐标分别包含 20 个索引(总共 40 个),而不是分别包含 4 个和 5 个索引(总共 9 个)?希望你能明白我困惑的地方。

是否有可能以某种方式映射这些 2D 数组并“降级”到 1450 个纬度值和 900 个经度值?还是像现在这样就可以了?我如何使用这些价值观来实现我的意图?我需要压缩经纬度数组吗?

这是形状:

>>> DS = xarray.open_dataset('file.nc')
>>> DS.tasmin.shape
    (31, 1450, 900)
>>> DS.projection_x_coordinate.shape
    (900,)
>>> DS.projection_y_coordinate.shape
    (1450,)
>>> DS.latitude.shape
    (1450, 900)
>>> DS.longitude.shape
    (1450, 900)
Run Code Online (Sandbox Code Playgroud)

考虑projection_x_coordinate东距projection_y_coordinate/北距值不是纬度/经度

如果需要的话,这是文件的元数据:

尺寸:(bnds:2,投影_x_坐标:900,投影_y_坐标:1450,时间:31)
坐标:
  * 时间(时间)datetime64[ns] 2018-12-01T12:00:00 ....
  * 投影 y 坐标 (投影 y 坐标) float64 -1.995e+0...
  * 投影 x 坐标 (投影 x 坐标) float64 -1.995e+0...
    纬度(投影 y 坐标、投影 x 坐标) float64 ...
    经度(投影 y 坐标、投影 x 坐标) float64 ...
无坐标尺寸:bnds
数据变量:
    tasmin(时间,投影_y_坐标,投影_x_坐标)float64 ...
    横向墨卡托 int32 ...
    time_bnds (时间, bnds) datetime64[ns] ...
    投影 y 坐标 bnds (投影 y 坐标, bnds) float64 ...
    投影_x_坐标_bnds(投影_x_坐标,bnds)float64 ...
属性:
    评论:每日分辨率网格化气候观测
    创建日期:2019-08-21T21:26:02
    频率:天
    机构:英国气象局
    参考文献:doi:10.1002/joc.1161
    短名称:daily_mintemp
    来源:HadUK-Grid_v1.0.1.0
    标题:英国网格化表面气候观测数据
    版本:v20190808
    约定:CF-1.5

ala*_*ani 7

您的数据符合 1.5 版气候和预报公约

描述此版本约定的文档位于此处,尽管相关部分在许多版本的约定中基本没有变化。

参见第 5.2 节:

5.2. 二维纬度、经度、坐标变量

未定义为纬度和经度轴的笛卡尔积的水平网格的纬度和经度坐标有时可以使用二维坐标变量来表示。使用坐标属性将这些变量标识为坐标。

看起来您正在使用 HadOBS 1 公里分辨率网格化每日最低温度,特别是此文件:

http://dap.ceda.ac.uk/thredds/fileServer/badc/ukmo-hadobs/data/insitu/MOHC/HadOBS/HadUK-Grid/v1.0.1.0/1km/tasmin/day/v20190808/tasmin_hadukgrid_uk_1km_day_20181201- 20181231.nc(警告:>300MB 下载)

正如其所述,数据位于横向墨卡托网格上。

如果您查看输出,ncdump -h <filename>您还会看到以下通过transverse_mercator虚拟变量的属性表示的网格描述:

        int transverse_mercator ;
                transverse_mercator:grid_mapping_name = "transverse_mercator" ;
                transverse_mercator:longitude_of_prime_meridian = 0. ;
                transverse_mercator:semi_major_axis = 6377563.396 ;
                transverse_mercator:semi_minor_axis = 6356256.909 ;
                transverse_mercator:longitude_of_central_meridian = -2. ;
                transverse_mercator:latitude_of_projection_origin = 49. ;
                transverse_mercator:false_easting = 400000. ;
                transverse_mercator:false_northing = -100000. ;
                transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ;
Run Code Online (Sandbox Code Playgroud)

您还将看到坐标变量projection_x_coordinateprojection_y_coordinate的单位为米。

所讨论的网格是使用数字网格参考的英国地形测量局网格。
例如,请参阅操作系统网格的描述(来自维基百科)。

如果您希望在常规经纬度网格上表达数据,那么您将需要进行某种类型的插值。我看到你正在使用 xarray。您可以将其与pyresample进行插值结合起来。这是一个例子:

import xarray as xr
import numpy as np
from pyresample.geometry import SwathDefinition
from pyresample.kd_tree import resample_nearest, resample_gauss

ds = xr.open_dataset("tasmin_hadukgrid_uk_1km_day_20181201-20181231.nc")

# Define a target grid. For sake of example, here is one with just 
# 3 longitudes and 4 latitudes.
lons = np.array([-2.1, -2., -1.9])
lats = np.array([51.7, 51.8, 51.9, 52.0])

# The target grid is regular (1-d lon, lat coordinates) but we will need
# a 2d version (similar to the input grid), so use numpy.meshgrid to produce this.
lon2d, lat2d = np.meshgrid(lons, lats)

origin_grid = SwathDefinition(lons=ds.longitude, lats=ds.latitude)
target_grid = SwathDefinition(lons=lon2d, lats=lat2d)

# get a numpy array for the first timestep
data = ds.tasmin[0].to_masked_array()

# nearest neighbour interpolation example
# Note that radius_of_influence has units metres

interpolated = resample_nearest(origin_grid, data, target_grid, radius_of_influence=1000)

# GIVES:
#      array([[5.12490065, 5.02715332, 5.36414835],
#             [5.08337723, 4.96372838, 5.00862833],
#             [6.47538931, 5.53855722, 5.11511239],
#             [6.46571817, 6.17949381, 5.87357538]])


# gaussian weighted interpolation example
# Note that radius_of_influence and sigmas both have units metres

interpolated = resample_gauss(origin_grid, data, target_grid, radius_of_influence=1000, sigmas=1000)

# GIVES:
#      array([[5.20432465, 5.07436805, 5.39693221],
#             [5.09069187, 4.8565934 , 5.08191639],
#             [6.4505963 , 5.44018209, 5.13774416],
#             [6.47345359, 6.2386732 , 5.62121948]])
Run Code Online (Sandbox Code Playgroud)