python - 尝试将不规则间隔的纬度/经度数据放置在规则间隔的网格上

Question

我有一些卫星数据，正在尝试将其插值到 0.25 度 x 0.25 度网格上。

我正在尝试使用scipy.intepolate.griddata，但我得到了意想不到的结果。

我只需要在卫星的扫描范围内进行插值。我不需要在整个全局范围内进行插值。

这是我的代码:

import numpy as np  
import scipy as sp 
import matplotlib as mpl  
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from scipy.interpolate import griddata
from pyhdf.SD import SD, SDC

hdf = SD(files[0], SDC.READ)
lon = hdf.select('Longitude')[:,:]
lat = hdf.select('Latitude')[:,:]
refl = hdf.select('correctZFactor')[:,:,70]/100

m = Basemap()

lonMin = -180
lonMax = 180
latMin = -40
latMax = 40
res = 1
lonGrid = np.arange(lonMin, lonMax, res)
latGrid = np.arange(latMin, latMax, res)
lonGrid,latGrid = np.meshgrid(lonGrid,latGrid)

reflGrid = griddata((lon.ravel(),lat.ravel()),refl.ravel(),(lonGrid,latGrid), method = 'nearest')

当我在网格化之前绘制数据时，它看起来像这样:

网格化后，如下所示:

这是我正在使用的 HDF 文件: http://www.filedropper.com/2a2520150314987057

显然，生成的图像未正确插值。我可以做什么来解决这个问题？

我的最终目标是获取数千个这样的卫星测绘带，每个卫星测绘带都经过全局不同的路径，并将它们组合成一个数据集。网格化到较粗分辨率的目的是 1. 减少所有数据量，2. 能够导出特定网格点的统计数据。另一件事:理想情况下， strip 之外的点在网格化后将转换为 NaN

Answer 1

正如我在评论中指出的那样，问题在于您的数据沿着非常狭窄的 strip ，使您的数据成为伪一维。如果您尝试从此“插值”到整个地球，那么您实际上是根据几乎不存在的值进行外推，这解释了原始图中的噪声。

由于您在编辑中澄清您只对在数据区域中进行插值感兴趣，所以我看到了不同类型的问题。沿着这条狭窄的经纬度点带的固定规则网格对我来说没有意义。查看原始数据的 pcolormesh 图:

import numpy as np
from pyhdf.SD import SD,SDC
import matplotlib.pyplot as plt
import scipy.interpolate as interp

hdffile = 'your_file_name.hdf'

hdf = SD(hdffile, SDC.READ)
lon = hdf.select('Longitude')[:,:]
lat = hdf.select('Latitude')[:,:]
refl = hdf.select('correctZFactor')[:,:,70]/100

lon[lon<0] += 360 # shift longitude to contiguous block

fig,ax = plt.subplots()
ax.pcolormesh(lon,lat,refl,cmap='viridis')

希望上面的图能够传达我的意思:尝试将此域放在常规网格上对于我能想到的任何合理用途都没有用处。特别是如果您认为给定经度的几度纬度宽度确实接近您预期的 0.25 度分辨率。

所以我的建议是在经度上采用常规网格，并且对于每个经度在域中的纬度上采用常规网格。这意味着您的最终网格不是规则的，但它将是拓扑二维格子的(就像由网格网格生成的一样)，因此它将对于绘图或其他后处理目的很有用。

为了做到这一点，我首先为每个经度的最小和最大纬度值构造两个插值器，然后生成 (lon,lat) 插值网格，然后进行插值:

# these will be overwritten later
lat_from = lat[:,0]
lat_to = lat[:,-1]
lon_from = lon[:,0]
lon_to = lon[:,-1]

# create interpolators for starting and ending latitude vs longitude
# only use a subset of the 9k data points
step = 10
latminfun = interp.interp1d(lon_from[::step],lat_from[::step],fill_value='extrapolate')
latmaxfun = interp.interp1d(lon_to[::step],lat_to[::step],fill_value='extrapolate')

# create interpolating mesh: regular in longitude, locally regular in latitude
nlon = 360 # ~1 degree along longitude
nlat = 10 # 10 points along latitude for each longitude
lon_grid = np.linspace(lon.min(),lon.max(),nlon)[:,None]  # shape (nlon,1)
lat_from = latminfun(lon_grid) # lower side of the latitude grid
lat_to = latmaxfun(lon_grid)   # upper side of the latitude grid
x = np.linspace(0,1,nlat) # to perform linear interpolation in lat with
lat_grid = x*lat_to + (1-x)*lat_from # shape (nlon,nlat)

# now (lon_grid,lat_grid) broadcast together to a grid of shape (nlon,nlat)
refl_grid = interp.griddata((lon.ravel(),lat.ravel()),refl.ravel(),(lon_grid,lat_grid),method='nearest')
fig,ax = plt.subplots()
ax.pcolormesh(np.broadcast_to(lon_grid,lat_grid.shape),lat_grid,refl_grid,cmap='viridis')
# of course we could've overwritten lon_grid with the broadcast version

最终图在视觉上几乎无法与原始数据区分开来:

但它包含此直线经纬度网格上的最近邻插值。我希望这是插入数据的最合理方法，无需了解有关结果的计划的任何详细信息。