尝试将不规则间隔的 lat/lon 数据放置在规则间隔的网格上
Trying to place irregularly spaced lat/lon data on regularly spaced grid
我有一些卫星数据，我正试图将其插入到 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. 能够导出特定网格点的统计数据。另一件事：理想情况下，在网格
之后，条带之外的点将转换为 NaN
正如我在评论中也指出的那样，问题是您的数据沿着一条非常窄的条带分布，使您的数据成为伪 1d。如果您尝试 "interpolate" 从这个到整个地球，您实际上是在根据几乎不存在的值进行外推，这解释了原始图中的噪声。
由于您在编辑中阐明了您只对数据区域中的插值感兴趣，所以我看到了另一种问题。沿着这条狭长的经纬度点的固定规则网格对我来说没有意义。查看原始数据的 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 度分辨率。
因此，我建议取而代之的是采用规则的经度网格，并且对于每个经度，在域中采用规则的纬度网格。这意味着您的最终网格不是规则的，但它在拓扑上是二维格子（就像由 meshgrid 生成的一样），因此它对于绘图或其他后处理目的很有用。
为了做到这一点，我首先为每个经度的最小和最大纬度值构造两个插值器，然后生成 (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
最终图在视觉上与您的原始数据几乎没有区别：
但它包含此直线经纬度网格上的最近邻插值。我希望这是插入数据的最合理方式，无需了解有关您的计划的任何详细信息。
我最终使用 KD-tree/nearest-neighbor 查找找到了解决方案
from scipy import spatial
kdtree = spatial.cKDTree(zip(lon.ravel(),lat.ravel()))
kdtree_gridPts = spatial.cKDTree(zip(lonGrid.ravel(),latGrid.ravel()))
closePts = kdtree_gridPts.query_ball_tree(kdtree, res/2)
reflGrid = ones_like(lonGrid)*nan
for ind,p in enumerate(closePts):
if len(p) > 0:
reflGrid.ravel()[ind] = mean(refl.ravel()[p])
reflGrid = ma.masked_where(isnan(reflGrid), reflGrid)
基本上它是以每个网格框为中心的 0.125 (res/2) 圆内所有点的平均值。至少，这就是我认为它正在做的...
我也将我的域缩小到我感兴趣的领域。它运行得非常快。
向下网格化之前：
网格化后：
预网格化图片顶部的水平线实际上看起来只是 pcolormesh 的人工产物
我有一些卫星数据，我正试图将其插入到 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. 能够导出特定网格点的统计数据。另一件事：理想情况下，在网格
之后，条带之外的点将转换为 NaN正如我在评论中也指出的那样，问题是您的数据沿着一条非常窄的条带分布，使您的数据成为伪 1d。如果您尝试 "interpolate" 从这个到整个地球，您实际上是在根据几乎不存在的值进行外推，这解释了原始图中的噪声。
由于您在编辑中阐明了您只对数据区域中的插值感兴趣，所以我看到了另一种问题。沿着这条狭长的经纬度点的固定规则网格对我来说没有意义。查看原始数据的 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 度分辨率。
因此，我建议取而代之的是采用规则的经度网格，并且对于每个经度，在域中采用规则的纬度网格。这意味着您的最终网格不是规则的，但它在拓扑上是二维格子（就像由 meshgrid 生成的一样），因此它对于绘图或其他后处理目的很有用。
为了做到这一点，我首先为每个经度的最小和最大纬度值构造两个插值器，然后生成 (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
最终图在视觉上与您的原始数据几乎没有区别：
但它包含此直线经纬度网格上的最近邻插值。我希望这是插入数据的最合理方式，无需了解有关您的计划的任何详细信息。
我最终使用 KD-tree/nearest-neighbor 查找找到了解决方案
from scipy import spatial
kdtree = spatial.cKDTree(zip(lon.ravel(),lat.ravel()))
kdtree_gridPts = spatial.cKDTree(zip(lonGrid.ravel(),latGrid.ravel()))
closePts = kdtree_gridPts.query_ball_tree(kdtree, res/2)
reflGrid = ones_like(lonGrid)*nan
for ind,p in enumerate(closePts):
if len(p) > 0:
reflGrid.ravel()[ind] = mean(refl.ravel()[p])
reflGrid = ma.masked_where(isnan(reflGrid), reflGrid)
基本上它是以每个网格框为中心的 0.125 (res/2) 圆内所有点的平均值。至少，这就是我认为它正在做的...
我也将我的域缩小到我感兴趣的领域。它运行得非常快。
向下网格化之前：
网格化后：
预网格化图片顶部的水平线实际上看起来只是 pcolormesh 的人工产物