为什么纬度和经度在netcdf文件中是二维数组？
Why latitudes and longitudes are two dimensional arrays in netcdf file?
我有 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 是 easting/northing 值而不是 lat/longs
如果需要，这里是文件的元数据：
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
我自己想出了一个答案。
出现的二维经纬度数组用于定义某个位置的“网格”。
换句话说，如果我们zip经纬度值在地图上投影，我们会在某个位置得到“弯曲的网格”（换句话说就是考虑地球曲率），然后用于创建位置的网格参考。
希望所有感兴趣的人都清楚。
您的数据符合 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 分辨率网格化每日最低温度，尤其是此文件：
如其所述，数据位于横向墨卡托网格上。
如果您查看 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的单位都是米
相关网格是使用数字网格引用的英国军械测量局网格。
例如，参见 OS 网格的 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]])
我有 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 是 easting/northing 值而不是 lat/longs
如果需要，这里是文件的元数据：
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
我自己想出了一个答案。
出现的二维经纬度数组用于定义某个位置的“网格”。
换句话说，如果我们zip经纬度值在地图上投影，我们会在某个位置得到“弯曲的网格”（换句话说就是考虑地球曲率），然后用于创建位置的网格参考。
希望所有感兴趣的人都清楚。
您的数据符合 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 分辨率网格化每日最低温度，尤其是此文件：
如其所述，数据位于横向墨卡托网格上。
如果您查看 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的单位都是米
相关网格是使用数字网格引用的英国军械测量局网格。
例如，参见 OS 网格的 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]])