python - 为什么纬度和经度在netcdf文件中是二维数组？

Question

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

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

所以我从 netCDf 文件中提取了纬度和经度坐标数据，但我期望它们是一维数组，而不是两个坐标都具有 1450x900 形状的二维数组。

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

假设我们有 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个)呢？希望你明白我的困惑。

是否有可能以某种方式映射这些二维数组并将其“降级”为 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)

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

如果需要，这里是文件的元数据:

Dimensions:                       (bnds: 2, projection_x_coordinate: 900, projection_y_coordinate: 1450, time: 31)
Coordinates:
  * time                          (time) datetime64[ns] 2018-12-01T12:00:00 ....
  * projection_y_coordinate       (projection_y_coordinate) float64 -1.995e+0...
  * projection_x_coordinate       (projection_x_coordinate) float64 -1.995e+0...
    latitude                      (projection_y_coordinate, projection_x_coordinate) float64 ...
    longitude                     (projection_y_coordinate, projection_x_coordinate) float64 ...
Dimensions without coordinates: bnds
Data variables:
    tasmin                        (time, projection_y_coordinate, projection_x_coordinate) float64 ...
    transverse_mercator           int32 ...
    time_bnds                     (time, bnds) datetime64[ns] ...
    projection_y_coordinate_bnds  (projection_y_coordinate, bnds) float64 ...
    projection_x_coordinate_bnds  (projection_x_coordinate, bnds) float64 ...
Attributes:
    comment:        Daily resolution gridded climate observations
    creation_date:  2019-08-21T21:26:02
    frequency:      day
    institution:    Met Office
    references:     doi: 10.1002/joc.1161
    short_name:     daily_mintemp
    source:         HadUK-Grid_v1.0.1.0
    title:          Gridded surface climate observations data for the UK
    version:        v20190808
    Conventions:    CF-1.5

Answer 1

您的数据符合 Climate and Forecast conventions 的 1.5 版.

描述此版本约定的文档是 here ，尽管相关部分在公约的许多版本中基本没有变化。

参见第 5.2 节:

5.2. Two-Dimensional Latitude, Longitude, Coordinate Variables

The latitude and longitude coordinates of a horizontal grid that was not defined as a Cartesian product of latitude and longitude axes, can sometimes be represented using two-dimensional coordinate variables. These variables are identified as coordinates by use of the coordinates attribute.

看起来您正在使用 HadOBS 1km 分辨率网格化每日最低温度，尤其是此文件:

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 ;

你还会看到坐标变量projection_x_coordinate和 projection_y_coordinate单位为米。

所讨论的网格是使用数字网格引用的英国军械测量局网格。
参见例如这个 description操作系统网格(来自维基百科)。

如果您希望在规则的经纬度网格上表达数据，则需要进行某种类型的插值。我看到您正在使用 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]])