I have the following issue:
I have cloud cover datasets from different satellites, that I want to regrid onto the grid of a climate model to make comparisons between the model output and the observed satellite data.
For now I was using the interp function from basemap, which works perfectly fine for arrays in the shape of: 1 x longitude x latitude, but it doesn't work for arrays of the shape n x longitude x latitude. What would be the best way to regrid these sort of 3-D arrays?
from mpl_toolkits.basemap import interp
import xarray as xr
def regrid_MAR10km(x_in,y_in,data_in):
mar_10km = xr.open_dataset('/media/..../MAR_10km /MARv3.5.2-10km-ERA-2008.nc')
lat = mar_10km['LAT']
lon = mar_10km['LON']
result = interp(data_in, x_in, y_in,lon,lat)
return result
My problem is, that I get error messages of the following form when I try to use 3-D data, in my case the array is of the form 161 (which is the cloud cover for every month!) x lon x lat
<xarray.DataArray 'Cloud_Fraction_Mask_Total_Mean' (Begin_Date: 161, lat: 28, lon: 110)>
[495880 values with dtype=float64]
Coordinates:
* Begin_Date (Begin_Date) datetime64[ns] 2002-07-01 2002-08-01 2002-09-01 ...
* lat (lat) float32 85.5 84.5 83.5 82.5 81.5 80.5 79.5 78.5 77.5 ...
* lon (lon) float32 -94.5 -93.5 -92.5 -91.5 -90.5 -89.5 -88.5 ...
And that is the error that it produces:
IndexError Traceback (most recent call last)
<ipython-input-19-5117a62e570e> in <module>()
----> 1 cloud_regrid = fc.regrid_MAR10km(longitude,latitude,cloud_data_UD)
/media/sf_Shared/Black_and_bloom/CODE/functions.py in regrid_MAR10km(x_in, y_in, data_in)
20 lat = mar_10km['LAT']
21 lon = mar_10km['LON']
---> 22 result = interp(data_in, x_in, y_in,lon,lat)
23 return result
/home/sh16450/miniconda3/envs/snowflakes/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py in interp(datain, xin, yin, xout, yout, checkbounds, masked, order)
4958 dataout = (1.-delx)*(1.-dely)*datain[yi,xi] + \
4959 delx*dely*datain[yip1,xip1] + \
-> 4960 (1.-delx)*dely*datain[yip1,xi] + \
4961 delx*(1.-dely)*datain[yi,xip1]
4962 elif order == 0:
IndexError: index 41 is out of bounds for axis 1 with size 28
From the Basemap docs the interp method deals exclusively with 2-D arrays where the 1st dimension corresponds to the y dimension (latitudes) and the 2nd to x (longitudes). If you give it anything else I'm afraid you'll end up with poor results.
What's actually happening is that since Basemap assumes your data is 2-D it's trying to perform the interpolation on your time and latitude (161,28) axes as if they were latitude and longitude respectively. Since the longitude array (110,) you passed to the method is now larger than the axis Basemap assumes to be longitudes (28,), you end up with a shape mismatch that leads to this IndexError.
I'm not exactly sure why the (1,28,110) works, but I wager that there's some array flattening or squeezing that NumPy does under the covers which ends up basically ignoring that empty dimension.
I recommend either iterating through each of the time steps and performing the interpolation on that 2-D slice of data or searching through the scipy.interpolate docs to see if they offer a function to do what you're looking for.
Good luck!