5

我正在尝试将具有相关纬度和经度的不规则网格数据集(原始卫星数据)映射到由basemap.makegrid(). 我正在使用matplotlib.mlab.griddata安装mpl_toolkits.natgridwhos以下是ipython中用作输出的变量列表以及变量的一些统计信息:

Variable   Type       Data/Info
-------------------------------
datalat    ndarray    666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)
datalon    ndarray    666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)
gridlat    ndarray    1200x1000: 1200000 elems, type `float64`, 9600000 bytes (9 Mb)
gridlon    ndarray    1200x1000: 1200000 elems, type `float64`, 9600000 bytes (9 Mb)
var        ndarray    666x1081: 719946 elems, type `float32`, 2879784 bytes (2 Mb)

In [11]: var.min()
Out[11]: -30.0

In [12]: var.max()
Out[12]: 30.0

In [13]: datalat.min()
Out[13]: 27.339874

In [14]: datalat.max()
Out[14]: 47.05302

In [15]: datalon.min()
Out[15]: -137.55658

In [16]: datalon.max()
Out[16]: -108.41629

In [17]: gridlat.min()
Out[17]: 30.394031556984299

In [18]: gridlat.max()
Out[18]: 44.237140350357713

In [19]: gridlon.min()
Out[19]: -136.17646180595321

In [20]: gridlon.max()
Out[20]: -113.82353819404671

datalat是原始datalon数据坐标

gridlat并且gridlon是要插值的坐标

var包含实际数据

使用这些变量,当我调用griddata(datalon, datalat, var, gridlon, gridlat)它时,它需要长达 20 分钟才能完成并返回一个nan. 通过查看数据,纬度和经度似乎是正确的,原始坐标与新区域的一部分和位于新区域之外的一些数据点重叠。有没有人有什么建议?nan 值表明我在做一些愚蠢的事情......

4

4 回答 4

2

很可能,griddata 太难了。它旨在处理随机采样的数据。您的数据几乎肯定会定期采样——只是不在与目标输出网格相同的网格上。

如果地球的拓扑或曲率影响 yoru 结果,请查看更简单的方法,例如仿射变换或小芯片上的一系列仿射变换。

有一些开箱即用的解决方案可能会有所帮助。 GDAL就是一个很好的例子。

此外,此类问题经常在 GIS 中讨论。看:

https://gis.stackexchange.com/questions/10430/changeing-image-projection-using-python

于 2011-09-15T01:13:26.160 回答
2

看起来该mlab.griddata例程可能会对您的输出数据引入额外的约束,这些约束可能不是必需的。虽然输入位置可以是任何东西,但输出位置必须是常规网格- 因为您的示例位于 lat/lon 空间,您选择的地图投影可能会违反这一点(即 x/y 中的常规网格不是 lat 中的常规网格/隆)。

您可以尝试使用SciPyinterpolate.griddata中的例程作为替代方案 - 但是,您需要将位置变量组合到一个数组中,因为调用签名不同:类似于

import scipy.interpolate
data_locations = np.vstack(datalon.ravel(), datalat.ravel()).T
grid_locations = np.vstack(gridlon.ravel(), gridlat.ravel()).T
grid_data      = scipy.interpolate.griddata(data_locations, val.ravel(),
                                            grid_locations, method='nearest')

用于最近邻插值。这会将位置放入一个数组中,其中 2 列对应于您的 2 个维度。您可能还想在地图投影的变换空间中执行插值。

于 2011-12-15T16:40:15.503 回答
1

如果您的数据位于网格上,因此点处的数据点位于(datalon[i], datalat[j])data[i,j],那么您可以使用scipy.interpolate.RectBivariateSpline代替griddata。不过,一些特定于地理的库可能会提供更多功能。

于 2011-09-20T21:16:26.220 回答
1

如果您使用 pclomesh,则不必进行任何类型的插值。pcolormesh 很乐意按照您在此处给出的方式接受数据结构:

from mpl_toolkits.basemap import Basemap
m = Basemap(-----)
x,y = m(datalon, datalat)
m.pcolormesh(x,y,var)
plt.show()

请使用它并告诉我这是否有效。

但是,当轨道数据重叠时,pcolormesh 会出现一些问题。请参考我的这个问题,你可能会发现一些有用的东西。

使用 pcolormesh 绘制轨道数据

于 2012-09-07T11:26:32.220 回答