如何将 MetPy 的横截面分析应用于具有二维纬度和经度的数据集？
How to apply MetPy's cross section analysis to a data set with 2-dimensional lat and lon?
我有 4 维数据（时间、深度、y 和 x），但纬度和经度都是二维数组。 y 和 x 只是索引，所以只是从 0、1...end 等开始的整数。
非常类似于MetPy提供的示例数据集：
https://unidata.github.io/MetPy/latest/examples/cross_section.html
不幸的是，这不是最可重现的，因为它非常特定于数据本身。但是我在横截面部分遇到了麻烦。我可以根据metpy解析数据，但是取横截面的时候报错：
data = ds.metpy.parse_cf().squeeze()
print(data)
<xarray.Dataset>
Dimensions: (y: 1347, x: 1379, deptht: 46, axis_nbounds: 2, time_counter: 12)
Coordinates:
nav_lat (y, x) float32 20.92 20.92 20.92 ... 68.49 68.49 68.48
nav_lon (y, x) float32 -78.95 -78.9 -78.85 ... -3.614 -3.546
* deptht (deptht) float32 3.047 9.454 ... 5.625e+03 5.875e+03
time_centered (time_counter) datetime64[ns] 1993-01-16T12:00:00 ....
* time_counter (time_counter) datetime64[ns] 1993-01-16T12:00:00 ....
Dimensions without coordinates: y, x, axis_nbounds
Data variables: (12/17)
deptht_bounds (deptht, axis_nbounds, y, x) float32 0.0 0.0 ... nan
time_centered_bounds (time_counter, axis_nbounds, y, x) datetime64[ns] 1...
time_counter_bounds (time_counter, axis_nbounds, y, x) datetime64[ns] 1...
votemper (time_counter, deptht, y, x) float32 26.63 ... nan
vosaline (time_counter, deptht, y, x) float32 35.88 ... nan
sosstsst (time_counter, y, x) float32 26.63 26.58 ... nan nan
... ...
sohefldo (time_counter, y, x) float32 -50.02 -45.85 ... nan nan
somixhgt (time_counter, y, x) float32 19.84 19.76 ... nan nan
sowindsp (time_counter, y, x) float32 5.591 5.48 ... nan nan
sohefldp (time_counter, y, x) float32 nan nan nan ... nan nan
sowafldp (time_counter, y, x) float32 3.94e-06 ... nan
sobowlin (time_counter, y, x) float32 20.04 20.04 ... nan nan
Attributes:
name: 1_VIKING20X.L46-KKG36107B_1d_19930101_19930704...
description: ocean T grid variables
title: ocean T grid variables
Conventions: CF-1.6
timeStamp: 2019-Sep-11 21:02:06 GMT
uuid: a091f081-2943-4c66-aa77-0917f654b802
history: Thu Sep 12 22:15:40 2019: ncrcat -O -F /gfs1/w...
NCO: 4.4.8
nco_openmp_thread_number: 1
然后我尝试横截面部分（还要注意，二维纬度和经度标记为 'nav_lat' 和 'nav_lon'，而不是纬度和经度）：
start = (40, -40)
end = (50, -30)
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
print(cross)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
165 try:
--> 166 crs_data = data.metpy.pyproj_crs
167 x = data.metpy.x
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/xarray.py in pyproj_crs(self)
252 """Return the coordinate reference system (CRS) as a pyproj object."""
--> 253 return self.crs.to_pyproj()
254
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/xarray.py in crs(self)
232 return self._data_array.coords['metpy_crs'].item()
--> 233 raise AttributeError('crs attribute is not available.')
234
AttributeError: crs attribute is not available.
During handling of the above exception, another exception occurred:
ValueError Traceback (most recent call last)
/tmp/ipykernel_2676704/3597335100.py in <module>
----> 1 cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
2 print(cross)
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
154 if isinstance(data, xr.Dataset):
155 # Recursively apply to dataset
--> 156 return data.map(cross_section, True, (start, end), steps=steps,
157 interp_type=interp_type)
158 elif data.ndim == 0:
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/xarray/core/dataset.py in map(self, func, keep_attrs, args, **kwargs)
5106 if keep_attrs is None:
5107 keep_attrs = _get_keep_attrs(default=False)
-> 5108 variables = {
5109 k: maybe_wrap_array(v, func(v, *args, **kwargs))
5110 for k, v in self.data_vars.items()
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/xarray/core/dataset.py in <dictcomp>(.0)
5107 keep_attrs = _get_keep_attrs(default=False)
5108 variables = {
-> 5109 k: maybe_wrap_array(v, func(v, *args, **kwargs))
5110 for k, v in self.data_vars.items()
5111 }
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
167 x = data.metpy.x
168 except AttributeError:
--> 169 raise ValueError('Data missing required coordinate information. Verify that '
170 'your data have been parsed by MetPy with proper x and y '
171 'dimension coordinates and added crs coordinate of the '
ValueError: Data missing required coordinate information. Verify that your data have been parsed by MetPy with proper x and y dimension coordinates and added crs coordinate of the correct projection for each variable.
我尝试了替代方案，因为可能 parse_cf 导致了问题：
data = ds.metpy.assign_latitude_longitude(force=True).squeeze()
但是我在应用横截面时仍然出现同样的错误。
关于如何解决这个问题有什么想法吗？
再次抱歉缺乏可重复性，但任何想法都会有很大帮助:)
可能跟投影有关？
这是一个示例图像，说明数据如何看待一个时间和深度实例（观察温度）：
metpy.interpolate.cross_section requires that your data include both x and y dimension coordinates and the added metpy_crs coordinate (from either parse_cf or assign_crs). In this situation where these x and y dimension coordinates are missing, but you do have 2D latitude and longitude coordinates, these dimension coordinates can be calculated and added with .metpy.assign_y_x()（而不是你说你试过的 assign_latitude_longitude，它的作用恰恰相反——从 y/x 维度坐标添加 lat/lon 辅助坐标）。
因此，如果您的数据集具有与您的数据投影相对应的有效 CF 网格映射，您将拥有：
data = ds.metpy.parse_cf()
data = data.metpy.assign_y_x()
start = (40, -40)
end = (50, -30)
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
如果没有：
data = ds.metpy.assign_crs({
"grid_mapping_name": "lambert_conformal_conic",
"standard_parallel": 25.0,
"longitude_of_central_meridian": 265.0,
"latitude_of_projection_origin": 25.0,
})
data = data.metpy.assign_y_x()
start = (40, -40)
end = (50, -30)
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
（变化 the projection attributes to those needed for your given projection）
我有 4 维数据（时间、深度、y 和 x），但纬度和经度都是二维数组。 y 和 x 只是索引，所以只是从 0、1...end 等开始的整数。 非常类似于MetPy提供的示例数据集：
https://unidata.github.io/MetPy/latest/examples/cross_section.html
不幸的是，这不是最可重现的，因为它非常特定于数据本身。但是我在横截面部分遇到了麻烦。我可以根据metpy解析数据，但是取横截面的时候报错：
data = ds.metpy.parse_cf().squeeze()
print(data)
<xarray.Dataset>
Dimensions: (y: 1347, x: 1379, deptht: 46, axis_nbounds: 2, time_counter: 12)
Coordinates:
nav_lat (y, x) float32 20.92 20.92 20.92 ... 68.49 68.49 68.48
nav_lon (y, x) float32 -78.95 -78.9 -78.85 ... -3.614 -3.546
* deptht (deptht) float32 3.047 9.454 ... 5.625e+03 5.875e+03
time_centered (time_counter) datetime64[ns] 1993-01-16T12:00:00 ....
* time_counter (time_counter) datetime64[ns] 1993-01-16T12:00:00 ....
Dimensions without coordinates: y, x, axis_nbounds
Data variables: (12/17)
deptht_bounds (deptht, axis_nbounds, y, x) float32 0.0 0.0 ... nan
time_centered_bounds (time_counter, axis_nbounds, y, x) datetime64[ns] 1...
time_counter_bounds (time_counter, axis_nbounds, y, x) datetime64[ns] 1...
votemper (time_counter, deptht, y, x) float32 26.63 ... nan
vosaline (time_counter, deptht, y, x) float32 35.88 ... nan
sosstsst (time_counter, y, x) float32 26.63 26.58 ... nan nan
... ...
sohefldo (time_counter, y, x) float32 -50.02 -45.85 ... nan nan
somixhgt (time_counter, y, x) float32 19.84 19.76 ... nan nan
sowindsp (time_counter, y, x) float32 5.591 5.48 ... nan nan
sohefldp (time_counter, y, x) float32 nan nan nan ... nan nan
sowafldp (time_counter, y, x) float32 3.94e-06 ... nan
sobowlin (time_counter, y, x) float32 20.04 20.04 ... nan nan
Attributes:
name: 1_VIKING20X.L46-KKG36107B_1d_19930101_19930704...
description: ocean T grid variables
title: ocean T grid variables
Conventions: CF-1.6
timeStamp: 2019-Sep-11 21:02:06 GMT
uuid: a091f081-2943-4c66-aa77-0917f654b802
history: Thu Sep 12 22:15:40 2019: ncrcat -O -F /gfs1/w...
NCO: 4.4.8
nco_openmp_thread_number: 1
然后我尝试横截面部分（还要注意，二维纬度和经度标记为 'nav_lat' 和 'nav_lon'，而不是纬度和经度）：
start = (40, -40)
end = (50, -30)
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
print(cross)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
165 try:
--> 166 crs_data = data.metpy.pyproj_crs
167 x = data.metpy.x
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/xarray.py in pyproj_crs(self)
252 """Return the coordinate reference system (CRS) as a pyproj object."""
--> 253 return self.crs.to_pyproj()
254
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/xarray.py in crs(self)
232 return self._data_array.coords['metpy_crs'].item()
--> 233 raise AttributeError('crs attribute is not available.')
234
AttributeError: crs attribute is not available.
During handling of the above exception, another exception occurred:
ValueError Traceback (most recent call last)
/tmp/ipykernel_2676704/3597335100.py in <module>
----> 1 cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
2 print(cross)
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
154 if isinstance(data, xr.Dataset):
155 # Recursively apply to dataset
--> 156 return data.map(cross_section, True, (start, end), steps=steps,
157 interp_type=interp_type)
158 elif data.ndim == 0:
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/xarray/core/dataset.py in map(self, func, keep_attrs, args, **kwargs)
5106 if keep_attrs is None:
5107 keep_attrs = _get_keep_attrs(default=False)
-> 5108 variables = {
5109 k: maybe_wrap_array(v, func(v, *args, **kwargs))
5110 for k, v in self.data_vars.items()
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/xarray/core/dataset.py in <dictcomp>(.0)
5107 keep_attrs = _get_keep_attrs(default=False)
5108 variables = {
-> 5109 k: maybe_wrap_array(v, func(v, *args, **kwargs))
5110 for k, v in self.data_vars.items()
5111 }
~/miniconda3/envs/py3_std_maps/lib/python3.9/site-packages/metpy/interpolate/slices.py in cross_section(data, start, end, steps, interp_type)
167 x = data.metpy.x
168 except AttributeError:
--> 169 raise ValueError('Data missing required coordinate information. Verify that '
170 'your data have been parsed by MetPy with proper x and y '
171 'dimension coordinates and added crs coordinate of the '
ValueError: Data missing required coordinate information. Verify that your data have been parsed by MetPy with proper x and y dimension coordinates and added crs coordinate of the correct projection for each variable.
我尝试了替代方案，因为可能 parse_cf 导致了问题：
data = ds.metpy.assign_latitude_longitude(force=True).squeeze()
但是我在应用横截面时仍然出现同样的错误。
关于如何解决这个问题有什么想法吗？ 再次抱歉缺乏可重复性，但任何想法都会有很大帮助:) 可能跟投影有关？
这是一个示例图像，说明数据如何看待一个时间和深度实例（观察温度）：
metpy.interpolate.cross_section requires that your data include both x and y dimension coordinates and the added metpy_crs coordinate (from either parse_cf or assign_crs). In this situation where these x and y dimension coordinates are missing, but you do have 2D latitude and longitude coordinates, these dimension coordinates can be calculated and added with .metpy.assign_y_x()（而不是你说你试过的 assign_latitude_longitude，它的作用恰恰相反——从 y/x 维度坐标添加 lat/lon 辅助坐标）。
因此，如果您的数据集具有与您的数据投影相对应的有效 CF 网格映射，您将拥有：
data = ds.metpy.parse_cf()
data = data.metpy.assign_y_x()
start = (40, -40)
end = (50, -30)
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
如果没有：
data = ds.metpy.assign_crs({
"grid_mapping_name": "lambert_conformal_conic",
"standard_parallel": 25.0,
"longitude_of_central_meridian": 265.0,
"latitude_of_projection_origin": 25.0,
})
data = data.metpy.assign_y_x()
start = (40, -40)
end = (50, -30)
cross = cross_section(data, start, end).set_coords(('nav_lon', 'nav_lat'))
（变化 the projection attributes to those needed for your given projection）