关于metpy横截面坐标的问题
Question about the coordinates of metpy cross section
我按照代码在这个网站上画了横截面图。
(https://unidata.github.io/MetPy/latest/examples/cross_section.html#sphx-glr-examples-cross-section-py)
在示例中，运行使用以下代码将产生此结果（交叉）。
cross = cross_section(data, start, end).set_coords(('lat', 'lon'))
enter image description here
但是，与示例不同的是，运行宁横截面代码后结果值的x,y坐标发生了变化。并且经纬度是固定的。
enter image description here
一直不明白为什么会出现这样的结果
my code
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import os
import xarray as xr
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.interpolate import cross_section
os.chdir ("D:/PPL/CrossSection/20190404/")
ncfile = r'./2019040400pres.nc'
NC = xr.open_dataset(ncfile)
pres_list = ['1000','975','950','925','900','875','850','800','750','700','650','600','550','500','450','400','350','300','250','200','150','100','70','50'];
p_list = np.array(pres_list, dtype = np.float64)
NC = NC.metpy.parse_cf().squeeze()
NC = NC.assign_coords(isobaric = p_list).set_coords('isobaric')
我的代码必须根据 nc 文件除以等压层进行分析，因此我合并了为每一层划分的值。
像这样，
tmp1 = NC['TMP_1000mb']
tmp2 = NC['TMP_975mb']
...
tmp23 = NC['TMP_70mb']
tmp24 = NC['TMP_50mb']
TMP = xr.concat([tmp1,tmp2],'isobaric')
TMP = xr.concat([TMP,tmp3],'isobaric')
...
TMP = xr.concat([TMP,tmp23],'isobaric')
TMP = xr.concat([TMP,tmp24],'isobaric')
TMP['isobaric']=p_list
NC['Temperature'] = TMP
所以，我打印 NC。结果是
<xarray.Dataset>
Dimensions: (isobaric: 24, x: 602, y: 781)
Coordinates:
* y (y) float64 0.0 1.5e+03 3e+03 ... 1.168e+06 1.17e+06
* x (x) float64 0.0 1.5e+03 3e+03 ... 9e+05 9.015e+05
latitude (y, x) float64 32.26 32.26 32.26 ... 42.94 42.94 42.93
longitude (y, x) float64 121.8 121.9 121.9 ... 132.5 132.5 132.5
time datetime64[ns] 2019-04-04
metpy_crs object Projection: latitude_longitude
* isobaric (isobaric) float64 1e+03 975.0 950.0 ... 100.0 70.0 50.0
Data variables:
...
...
Temperature (isobaric, y, x) float32 282.87515 282.87515 ... nan nan
Attributes:
Conventions: CF-1.0
History: created by wgrib2
GRIB2_grid_template: 30
但是我运行横截面的时候，x,y改变了起点和终点，原来的经纬度是固定的
start = (37.5, 126.63)
end = (37.8 ,128.86)
cross = cross_section(NC, start, end)
#I didn't write .set_codes here because unlike the example, coordinates are given
#latitude and longitude.
print(cross)
<xarray.Dataset>
Dimensions: (index: 100, isobaric: 24)
Coordinates:
latitude (index) float64 32.26 32.26 32.26 ... 32.26 32.26 32.26
longitude (index) float64 121.8 121.8 121.8 ... 121.8 121.8 121.8
time datetime64[ns] 2019-04-04
metpy_crs object Projection: latitude_longitude
x (index) float64 126.6 126.7 126.7 ... 128.8 128.8 128.9
y (index) float64 37.5 37.5 37.51 37.51 ... 37.79 37.8 37.8
* index (index) int32 0 1 2 3 4 5 6 7 ... 92 93 94 95 96 97 98 99
* isobaric (isobaric) float64 1e+03 975.0 950.0 ... 100.0 70.0 50.0
Data variables:
...
...
Temperature (isobaric, index) float64 282.9 282.9 ... 209.6 209.6
Attributes:
Conventions: CF-1.0
History: created by wgrib2
GRIB2_grid_template: 30
我希望 x,y 保持不变，latitude 和 longitude 保持不变示例中的开始和结束范围。
我真的不明白为什么会这样。
感谢您的关注。
供参考，我附上数据源和初始数据。
数据来源：来自 KMA
的 LDAPS 数据
Full NC data(initial data)
<xarray.Dataset>
Dimensions: (time: 1, x: 602, y: 781)
Coordinates:
* y (y) float64 0.0 1.5e+03 3e+03 ... 1.168e+06 1.17e+06
* x (x) float64 0.0 1.5e+03 3e+03 ... 9e+05 9.015e+05
latitude (y, x) float64 ...
longitude (y, x) float64 ...
* time (time) datetime64[ns] 2019-04-04
Data variables:
DZDT_1000mb (time, y, x) float32 ...
DZDT_975mb (time, y, x) float32 ...
...
DZDT_70mb (time, y, x) float32 ...
DZDT_50mb (time, y, x) float32 ...
UGRD_1000mb (time, y, x) float32 ...
VGRD_1000mb (time, y, x) float32 ...
UGRD_975mb (time, y, x) float32 ...
VGRD_975mb (time, y, x) float32 ...
...
UGRD_70mb (time, y, x) float32 ...
VGRD_70mb (time, y, x) float32 ...
UGRD_50mb (time, y, x) float32 ...
VGRD_50mb (time, y, x) float32 ...
HGT_1000mb (time, y, x) float32 ...
HGT_975mb (time, y, x) float32 ...
...
HGT_70mb (time, y, x) float32 ...
HGT_50mb (time, y, x) float32 ...
TMP_1000mb (time, y, x) float32 ...
TMP_975mb (time, y, x) float32 ...
...
TMP_70mb (time, y, x) float32 ...
TMP_50mb (time, y, x) float32 ...
var0_1_194_1000mb (time, y, x) float32 ...
var0_1_194_975mb (time, y, x) float32 ...
...
var0_1_194_70mb (time, y, x) float32 ...
var0_1_194_50mb (time, y, x) float32 ...
RH_1000mb (time, y, x) float32 ...
RH_975mb (time, y, x) float32 ...
...
RH_70mb (time, y, x) float32 ...
RH_50mb (time, y, x) float32 ...
Attributes:
Conventions: CF-1.0
History: created by wgrib2
GRIB2_grid_template: 30
请注意，在 .metpy.parse_cf() 的数据集中标识的 metpy_crs 给出为 Projection: latitude_longitude，这是不正确的，因为您的数据集具有 2D latitude/longitude 坐标和 1D y/x 投影网格中的坐标 space。这可能是由于您的数据集的元数据不符合 CF 标准而触发的。现在，正因为如此，MetPy 将您的 x 坐标视为经度，将 y 视为纬度（注意这些坐标中的范围），即使它们不是。这里的结果肯定远非理想（没有错误消息或任何东西），所以我建议在 MetPy 的问题跟踪器上提交问题。
要直接解决您的问题，您需要手动指定数据CRS/projection。为此（在 MetPy v1.0 及更高版本中），请使用 assign_crs method on MetPy's Dataset accessor 而不是使用 parse_cf（仅适用于符合 CF 的数据集）。
我按照代码在这个网站上画了横截面图。 (https://unidata.github.io/MetPy/latest/examples/cross_section.html#sphx-glr-examples-cross-section-py)
在示例中，运行使用以下代码将产生此结果（交叉）。
cross = cross_section(data, start, end).set_coords(('lat', 'lon'))
enter image description here
但是，与示例不同的是，运行宁横截面代码后结果值的x,y坐标发生了变化。并且经纬度是固定的。 enter image description here 一直不明白为什么会出现这样的结果
my code
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import os
import xarray as xr
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.interpolate import cross_section
os.chdir ("D:/PPL/CrossSection/20190404/")
ncfile = r'./2019040400pres.nc'
NC = xr.open_dataset(ncfile)
pres_list = ['1000','975','950','925','900','875','850','800','750','700','650','600','550','500','450','400','350','300','250','200','150','100','70','50'];
p_list = np.array(pres_list, dtype = np.float64)
NC = NC.metpy.parse_cf().squeeze()
NC = NC.assign_coords(isobaric = p_list).set_coords('isobaric')
我的代码必须根据 nc 文件除以等压层进行分析，因此我合并了为每一层划分的值。 像这样，
tmp1 = NC['TMP_1000mb']
tmp2 = NC['TMP_975mb']
...
tmp23 = NC['TMP_70mb']
tmp24 = NC['TMP_50mb']
TMP = xr.concat([tmp1,tmp2],'isobaric')
TMP = xr.concat([TMP,tmp3],'isobaric')
...
TMP = xr.concat([TMP,tmp23],'isobaric')
TMP = xr.concat([TMP,tmp24],'isobaric')
TMP['isobaric']=p_list
NC['Temperature'] = TMP
所以，我打印 NC。结果是
<xarray.Dataset>
Dimensions: (isobaric: 24, x: 602, y: 781)
Coordinates:
* y (y) float64 0.0 1.5e+03 3e+03 ... 1.168e+06 1.17e+06
* x (x) float64 0.0 1.5e+03 3e+03 ... 9e+05 9.015e+05
latitude (y, x) float64 32.26 32.26 32.26 ... 42.94 42.94 42.93
longitude (y, x) float64 121.8 121.9 121.9 ... 132.5 132.5 132.5
time datetime64[ns] 2019-04-04
metpy_crs object Projection: latitude_longitude
* isobaric (isobaric) float64 1e+03 975.0 950.0 ... 100.0 70.0 50.0
Data variables:
...
...
Temperature (isobaric, y, x) float32 282.87515 282.87515 ... nan nan
Attributes:
Conventions: CF-1.0
History: created by wgrib2
GRIB2_grid_template: 30
但是我运行横截面的时候，x,y改变了起点和终点，原来的经纬度是固定的
start = (37.5, 126.63)
end = (37.8 ,128.86)
cross = cross_section(NC, start, end)
#I didn't write .set_codes here because unlike the example, coordinates are given
#latitude and longitude.
print(cross)
<xarray.Dataset>
Dimensions: (index: 100, isobaric: 24)
Coordinates:
latitude (index) float64 32.26 32.26 32.26 ... 32.26 32.26 32.26
longitude (index) float64 121.8 121.8 121.8 ... 121.8 121.8 121.8
time datetime64[ns] 2019-04-04
metpy_crs object Projection: latitude_longitude
x (index) float64 126.6 126.7 126.7 ... 128.8 128.8 128.9
y (index) float64 37.5 37.5 37.51 37.51 ... 37.79 37.8 37.8
* index (index) int32 0 1 2 3 4 5 6 7 ... 92 93 94 95 96 97 98 99
* isobaric (isobaric) float64 1e+03 975.0 950.0 ... 100.0 70.0 50.0
Data variables:
...
...
Temperature (isobaric, index) float64 282.9 282.9 ... 209.6 209.6
Attributes:
Conventions: CF-1.0
History: created by wgrib2
GRIB2_grid_template: 30
我希望 x,y 保持不变，latitude 和 longitude 保持不变示例中的开始和结束范围。 我真的不明白为什么会这样。
感谢您的关注。
供参考，我附上数据源和初始数据。 数据来源：来自 KMA
的 LDAPS 数据Full NC data(initial data)
<xarray.Dataset>
Dimensions: (time: 1, x: 602, y: 781)
Coordinates:
* y (y) float64 0.0 1.5e+03 3e+03 ... 1.168e+06 1.17e+06
* x (x) float64 0.0 1.5e+03 3e+03 ... 9e+05 9.015e+05
latitude (y, x) float64 ...
longitude (y, x) float64 ...
* time (time) datetime64[ns] 2019-04-04
Data variables:
DZDT_1000mb (time, y, x) float32 ...
DZDT_975mb (time, y, x) float32 ...
...
DZDT_70mb (time, y, x) float32 ...
DZDT_50mb (time, y, x) float32 ...
UGRD_1000mb (time, y, x) float32 ...
VGRD_1000mb (time, y, x) float32 ...
UGRD_975mb (time, y, x) float32 ...
VGRD_975mb (time, y, x) float32 ...
...
UGRD_70mb (time, y, x) float32 ...
VGRD_70mb (time, y, x) float32 ...
UGRD_50mb (time, y, x) float32 ...
VGRD_50mb (time, y, x) float32 ...
HGT_1000mb (time, y, x) float32 ...
HGT_975mb (time, y, x) float32 ...
...
HGT_70mb (time, y, x) float32 ...
HGT_50mb (time, y, x) float32 ...
TMP_1000mb (time, y, x) float32 ...
TMP_975mb (time, y, x) float32 ...
...
TMP_70mb (time, y, x) float32 ...
TMP_50mb (time, y, x) float32 ...
var0_1_194_1000mb (time, y, x) float32 ...
var0_1_194_975mb (time, y, x) float32 ...
...
var0_1_194_70mb (time, y, x) float32 ...
var0_1_194_50mb (time, y, x) float32 ...
RH_1000mb (time, y, x) float32 ...
RH_975mb (time, y, x) float32 ...
...
RH_70mb (time, y, x) float32 ...
RH_50mb (time, y, x) float32 ...
Attributes:
Conventions: CF-1.0
History: created by wgrib2
GRIB2_grid_template: 30
请注意，在 .metpy.parse_cf() 的数据集中标识的 metpy_crs 给出为 Projection: latitude_longitude，这是不正确的，因为您的数据集具有 2D latitude/longitude 坐标和 1D y/x 投影网格中的坐标 space。这可能是由于您的数据集的元数据不符合 CF 标准而触发的。现在，正因为如此，MetPy 将您的 x 坐标视为经度，将 y 视为纬度（注意这些坐标中的范围），即使它们不是。这里的结果肯定远非理想（没有错误消息或任何东西），所以我建议在 MetPy 的问题跟踪器上提交问题。
要直接解决您的问题，您需要手动指定数据CRS/projection。为此（在 MetPy v1.0 及更高版本中），请使用 assign_crs method on MetPy's Dataset accessor 而不是使用 parse_cf（仅适用于符合 CF 的数据集）。