0

尝试从 THREDDS 数据服务器绘制表面观测结果时，我的脚本挂起或返回KeyError. 下面是代码：

metarfile = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/'
                       'nws/metar/ncdecoded/files/latest.xml')
latestmetar = metarfile.datasets[0]

file = latestmetar.remote_access(use_xarray=True)

parsed_temp = file.metpy.parse_cf('air_temperature')
parsed_td = file.metpy.parse_cf('dew_point_temperature')
parsed_mslp = file.metpy.parse_cf('air_pressure_at_sea_level')

parsed_temp = parsed_temp * 9/5 + 32
parsed_td = parsed_td * 9/5 + 32

prj = ccrs.LambertConformal(central_latitude=35,
                            central_longitude=-98,
                            standard_parallels=(30, 60))

points = prj.transform_points(ccrs.PlateCarree(), file['longitude'].values,
                                                  file['latitude'].values)

file[mpcalc.reduce_point_density(points, 100000)]

fig = plt.figure(1, figsize=(10,10))
ax = fig.add_subplot(1, 1, 1, projection=prj)

ax.set_extent((-104.1, -95.5, 32.1, 39.1))

ax.add_feature(cfeature.STATES.with_scale('50m'),linewidth=.5,
                                                 edgecolor='black',
                                                 zorder=5)

metarplots = StationPlot(ax, file['longitude'].values, file['latitude'].values,
                             clip_on=True, transform=ccrs.PlateCarree(),
                             fontsize=12)

metarplots.plot_parameter('NW', file['air_temperature'], color='red')
metarplots.plot_parameter('SW', file['dew_point_temperature'], color='green')
metarplots.plot_parameter('NE', file['air_pressure_at_sea_level'],
                                color='yellow')

当散列 MetPyreduce_point_density以解决数据点过多的问题时，代码挂起且没有错误消息。此外，当从中打印数据时，parsed_td会返回我想要的值。

<xarray.DataArray 'dew_point_temperature' (recNum: 15331)>
array([69.8 , 41.  , 32.  , ..., 78.98, 60.8 , 62.96], dtype=float32)
Dimensions without coordinates: recNum

我真的不确定从哪里开始解决这个问题，因为我之前没有绘制点数据。任何帮助都会有所帮助！

编辑：包括错误信息。

Traceback (most recent call last):
  File "/anaconda3/lib/python3.7/site-packages/xarray/core/dataset.py", line 945, in _copy_listed
    variables[name] = self._variables[name]
KeyError: True

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "mslp_vis.py", line 63, in <module>
    file[mpcalc.reduce_point_density(points, 100000)]
  File "/anaconda3/lib/python3.7/site-packages/xarray/core/dataset.py", line 1060, in __getitem__
    return self._copy_listed(np.asarray(key))
  File "/anaconda3/lib/python3.7/site-packages/xarray/core/dataset.py", line 948, in _copy_listed
    self._variables, name, self._level_coords, self.dims)
  File "/anaconda3/lib/python3.7/site-packages/xarray/core/dataset.py", line 68, in _get_virtual_variable
    raise KeyError(key)
KeyError: True
4

1 回答 1

0

这里的问题是如何编写 netCDF。它采用离散采样几何 (DSG) 格式，不同于典型的网格格式。因此，此处计算的掩码仅对 Lat 和 Long 变量有效。对于数据变量，必须根据计算初始掩码后剩余的有效台站计算新掩码。下面的代码改编了原始代码来演示这个完整的过程。请注意，循环很慢，也许可以重构。

metarfile = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/'
                       'nws/metar/ncdecoded/files/latest.xml')
latestmetar = metarfile.datasets[0]

file = latestmetar.remote_access(use_xarray=True)

parsed_temp = file.metpy.parse_cf('air_temperature')
parsed_td = file.metpy.parse_cf('dew_point_temperature')
parsed_mslp = file.metpy.parse_cf('air_pressure_at_sea_level')

parsed_temp = parsed_temp * 9/5 + 32
parsed_td = parsed_td * 9/5 + 32

prj = ccrs.LambertConformal(central_latitude=35,
                            central_longitude=-98,
                            standard_parallels=(30, 60))

points = prj.transform_points(ccrs.PlateCarree(), file['longitude'].values,
                                                  file['latitude'].values)

points[2477, 0:2] = [0,0]
mask = mpcalc.reduce_point_density(points, 100000)

mask_stations = []
for x in file['parent_index'].values:
    if x in file['station'][mask]:
        mask_stations.append(True)
    else:
        mask_stations.append(False)

fig = plt.figure(1, figsize=(10,10))
ax = fig.add_subplot(1, 1, 1, projection=prj)

ax.set_extent((-104.1, -95.5, 32.1, 39.1))

ax.add_feature(cfeature.STATES.with_scale('50m'),linewidth=.5,
                                                 edgecolor='black',
                                                 zorder=5)

metarplots = StationPlot(ax, file['longitude'].values[mask], file['latitude'].values[mask],
                             clip_on=True, transform=ccrs.PlateCarree(),
                             fontsize=12)

metarplots.plot_parameter('NW', file['air_temperature'][mask_stations].metpy.convert_units('degF'), color='red')
metarplots.plot_parameter('SW', file['dew_point_temperature'][mask_stations], color='green')
metarplots.plot_parameter('NE', file['air_pressure_at_sea_level'][mask_stations],
                                color='blue')
于 2019-10-07T19:21:33.677 回答