I am trying to map an irregularly gridded dataset (raw satellite data) with associated latitudes and longitudes to a regularly gridded set of latitudes and longitudes given by basemap.makegrid()
. I am using matplotlib.mlab.griddata
with mpl_toolkits.natgrid
installed. Below is a list of the variables being used as output by whos
in ipython and some stats on the variables:
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
and datalon
are the orignal data coordinates
gridlat
and gridlon
are the coordinates to interpolate to
var
contains the actual data
Using these variables, when I call griddata(datalon, datalat, var, gridlon, gridlat)
it has taken as long as 20 minutes to complete and returns an array of nan
. From looking at the data, the latitudes and longitudes appear to be correct with the original coordinates overlapping a portion of the new area and a few data points lying outside of the new area. Does anyone have any suggestions? The nan values suggest that I'm doing something stupid...
It looks like the mlab.griddata
routine may introduce additional constraints on your output data that may not be necessary. While the input locations may be anything, the output locations must be a regular grid - since your example is in lat/lon space, your choice of map projection may violate this (i.e. regular grid in x/y is not a regular grid in lat/lon).
You might try the interpolate.griddata
routine from SciPy as an alternative - you'll need to combine your location variables into a single array, though, since the call signature is different: something like
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')
for nearest-neighbor interpolation. This gets the locations into an array with 2 columns corresponding to your 2 dimensions. You may also want to perform the interpolation in the transformed space of your map projection.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With