关于metpy横截面坐标的问题

如何解决关于metpy横截面坐标的问题

我按照代码在这个网站上绘制了横截面图。 (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 不明白为什么会出现这样的结果。

我的代码

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],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,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 保持不变，并且纬度和经度像示例一样位于开始和结束范围内。我真的不明白为什么会这样。

感谢您的关注。

作为参考，我附上数据来源和初始数据。数据来源：来自 KMA 的 LDAPS 数据

完整的NC数据（初始数据）

<xarray.Dataset>
Dimensions:            (time: 1,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,x) float32 ...
DZDT_975mb         (time,x) float32 ...
...
DZDT_70mb          (time,x) float32 ...
DZDT_50mb          (time,x) float32 ...
UGrd_1000mb        (time,x) float32 ...
VGrd_1000mb        (time,x) float32 ...
UGrd_975mb         (time,x) float32 ...
VGrd_975mb         (time,x) float32 ...
...
UGrd_70mb          (time,x) float32 ...
VGrd_70mb          (time,x) float32 ...
UGrd_50mb          (time,x) float32 ...
VGrd_50mb          (time,x) float32 ...
HGT_1000mb         (time,x) float32 ...
HGT_975mb          (time,x) float32 ...
...
HGT_70mb           (time,x) float32 ...
HGT_50mb           (time,x) float32 ...
TMP_1000mb         (time,x) float32 ...
TMP_975mb          (time,x) float32 ...
...
TMP_70mb           (time,x) float32 ...
TMP_50mb           (time,x) float32 ...
var0_1_194_1000mb  (time,x) float32 ...
var0_1_194_975mb   (time,x) float32 ...
...
var0_1_194_70mb    (time,x) float32 ...
var0_1_194_50mb    (time,x) float32 ...
RH_1000mb          (time,x) float32 ...
RH_975mb           (time,x) float32 ...
...
RH_70mb            (time,x) float32 ...
RH_50mb            (time,x) float32 ...
Attributes:
Conventions:          CF-1.0
History:              created by wgrib2
GRIB2_grid_template:  30

解决方法

请注意，metpy_crs 的数据集中标识的 .metpy.parse_cf() 被指定为 Projection: latitude_longitude，这是不正确的，因为您的数据集具有 2D 纬度/经度坐标和 1D y/x 坐标在投影网格空间中。这可能是由于您的数据集的元数据不符合 CF 标准而触发的。现在，因此，MetPy 将您的 x 坐标视为经度，将 y 视为纬度（注意这些坐标中的范围），即使它们不是。这里的结果绝对不是理想的（没有错误消息或任何东西），所以我建议在 MetPy 的问题跟踪器上提交问题。

要直接解决您的问题，您需要手动指定数据的 CRS/投影。为此（在 MetPy v1.0 及更高版本中），请使用 assign_crs method on MetPy's Dataset accessor 而不是使用 parse_cf（仅适用于符合 CF 的数据集）。