Interpolate_to_grid 返回所有 nan

Question

使用 MetPy Monday interpolate_to_grid 练习元数据，我成功地让 mslp 网格工作。 转到潜在温度，结果全是 nan。当它“工作”时。当它不起作用时，我会收到一组似乎没有帮助的错误...

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from siphon.catalog import TDSCatalog
from metpy.io import parse_metar_file
from metpy.interpolate import interpolate_to_grid, remove_nan_observations
from metpy.plots import add_metpy_logo, current_weather, sky_cover, StationPlot
from metpy.calc import wind_components, wet_bulb_temperature, altimeter_to_station_pressure,potential_temperature,gradient

from metpy.units import units
from datetime import datetime,timedelta
import pandas as pd

mapcrs = ccrs.LambertConformal(central_longitude=-100.,central_latitude=35.,standard_parallels=(30.,60.))
datacrs = ccrs.PlateCarree()

cat = TDSCatalog('https://thredds-test.unidata.ucar.edu/thredds/catalog/noaaport/text/metar/catalog.xml')
ds = cat.datasets[-4]

dattim = ds.name[6:14]+' '+ds.name[15:19]
ds.download()

df = parse_metar_file(ds.name)
#pandas dataframe
#df.head()
df.columns.values
extent = [-120,-72,24,50]

df = df.dropna(subset=['latitude','longitude','elevation','altimeter','air_temperature','eastward_wind','northward_wind','air_pressure_at_sea_level','dew_point_temperature'])

lon = df['longitude'].values
lat = df['latitude'].values
stn_ids = df['station_id'].values
elev = df['elevation'].values
altimeter = df['altimeter'].values
t2 = df['air_temperature'].values
mslp = df['air_pressure_at_sea_level'].values
#projected coords
xp, yp, _ = mapcrs.transform_points(datacrs,lon,lat).T # x,y returned

#mslp WORKS
x_masked, y_masked, mslp = remove_nan_observations(xp,yp,mslp)

#altgridx,altgridy,alt = interpolate_to_grid(x_masked,y_masked,alt, interp_type='cressman')
altgridx,altgridy,mslp = interpolate_to_grid(x_masked,y_masked,mslp, interp_type='barnes',gamma=.5,kappa_star=10, hres=25000)

#Potential Temperature doesnt work
pres = altimeter_to_station_pressure(altimeter * units('mbar'), elev * units('m'))*33.8639
print(pres)

# theta
x_masked, y_masked, temp = remove_nan_observations(xp,yp,t2*units('degC'))
x_masked, y_masked, pres = remove_nan_observations(xp,yp,pres)
print(np.size(temp))
potemp = potential_temperature(pres, temp)
print(np.size(potemp))
print(np.unique(np.array(potemp)))

grdx = 75000.
thgridx,thgridy,theta = interpolate_to_grid(x_masked,y_masked, potemp, interp_type='barnes',kappa_star=6, gamma=0.5,hres=grdx)

print(np.shape(thgridx))
print(np.unique(theta))

这是从最后一部分返回的内容：

[949.361081708803 993.4468013877739 987.2845093729651 ... 1029.0930108008558 1016.002484792407 930.3708063382303] millibar
5837
5837
[236.32885315 237.21299941 239.04372591 ... 368.37047837 369.20079652
 370.76269267]
---------------------------------------------------------------------------
DimensionalityError                       Traceback (most recent call last)
~/miniconda3/lib/python3.7/site-packages/pint/quantity.py in __float__(self)
    896             return float(self._convert_magnitude_not_inplace(UnitsContainer()))
--> 897         raise DimensionalityError(self._units, "dimensionless")
    898 

DimensionalityError: Cannot convert from 'kelvin' to 'dimensionless'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
/var/folders/5n/sg5k98bx6gg4flb4fskykh4m0000gn/T/ipykernel_41626/379842406.py in <module>
     11 
     12 grdx = 75000.
---> 13 thgridx,thgridy,theta = interpolate_to_grid(x_masked,y_masked, potemp, interp_type='barnes',kappa_star=6, gamma=0.5,hres=grdx)
     14 print(np.shape(thgridx))
     15 print(np.unique(theta))

~/miniconda3/lib/python3.7/site-packages/metpy/pandas.py in wrapper(*args, **kwargs)
     19         kwargs = {name: (v.values if isinstance(v, pd.Series) else v)
     20                   for name, v in kwargs.items()}
---> 21         return func(*args, **kwargs)
     22     return wrapper

~/miniconda3/lib/python3.7/site-packages/metpy/interpolate/grid.py in interpolate_to_grid(x, y, z, interp_type, hres, minimum_neighbors, gamma, kappa_star, search_radius, rbf_func, rbf_smooth, boundary_coords)
    301                                 minimum_neighbors=minimum_neighbors, gamma=gamma,
    302                                 kappa_star=kappa_star, search_radius=search_radius,
--> 303                                 rbf_func=rbf_func, rbf_smooth=rbf_smooth)
    304 
    305     return grid_x, grid_y, img.reshape(grid_x.shape)

~/miniconda3/lib/python3.7/site-packages/metpy/interpolate/points.py in interpolate_to_points(points, values, xi, interp_type, minimum_neighbors, gamma, kappa_star, search_radius, rbf_func, rbf_smooth)
    365             return inverse_distance_to_points(points, values, xi, search_radius, gamma, kappa,
    366                                               min_neighbors=minimum_neighbors,
--> 367                                               kind=interp_type)
    368 
    369     # If this is radial basis function, make the interpolator and apply it

~/miniconda3/lib/python3.7/site-packages/metpy/interpolate/points.py in inverse_distance_to_points(points, values, xi, r, gamma, kappa, min_neighbors, kind)
    268                 img[idx] = cressman_point(dists, values_subset, r)
    269             elif kind == 'barnes':
--> 270                 img[idx] = barnes_point(dists, values_subset, kappa, gamma)
    271 
    272             else:

ValueError: setting an array element with a sequence.

我在单位方面遇到了困难，但我认为现在单位是正确的。这可能是什么原因造成的？ 我尝试了 cressman，尝试了更大的 Barnes 网格，并尝试确保 search_radius 很大。仍然是 nan，当它起作用时。

Answer 1

这个问题是由interpolate_to_grid 在使用 Cressman 或 Barnes 时导致单位窒息引起的——我们definitely need to fix。目前的解决方案是使用不同的插值方法（如interp_type='linear'，默认值），或者在调用之前剥离单位：

thgridx, thgridy, theta = interpolate_to_grid(x_masked, y_masked, potemp.magnitude,
    interp_type='barnes', kappa_star=6, gamma=0.5, hres=grdx)
theta = units.Quantity(theta, 'K')

就NaNs 的问题而言，您可能需要查看search_radius 参数，该参数控制与目标点观测值的最大距离。在某些数据稀疏的地区，这可能会导致您有一些辍学。默认情况下，它使用从一个 ob 点到其最近邻居的平均距离的 5 倍进行猜测。