CODEWITHSUNDEEP
python-2.7matplotlibmatplotlib-basemappython-xarray
jetesdal
jetesdal

Reputation: 118

Interpolation between two xarray datasets with Basemap

I have two different xarray datasets that have different latitude/longitude grid resolutions. I want to regrid the one xarray with lower resolution to the same resolution as the one xarray with higher resolution. I found some examples (e.g., http://earthpy.org/interpolation_between_grids_with_basemap.html), but it does not work for me. Here is one example that I made for testing:

import numpy as np
import xarray as xray
import mpl_toolkits.basemap

var1=xray.DataArray(np.random.randn(len(np.linspace(40.5,49.5,10)),len(np.linspace(-39.5,-20.5,20))),coords=[np.linspace(40.5,49.5,10), np.linspace(-39.5,-20.5,20)],dims=['lat','lon'])

(xlon, xlat)=np.meshgrid(np.linspace(-39.875,-20.125,80),np.linspace(40.125,49.875,40))
var2=xray.DataArray(-xlon**2+xlat**2,coords=[np.linspace(40.125,49.875,40),np.linspace(-39.875,-20.125,80)],dims=['lat','lon'])

mpl_toolkits.basemap.interp(var1,var1.lon,var1.lat,var2.lon,var2.lat,checkbounds=False,masked=False,order=0)

I get following error:

ValueError: xout and yout must have same shape!

Screenshot: see screenshot

Does basemap.interp() require xout and yout to be the same shape? So var2 needs to be a square? This is almost never the case with any of my datasets! How can I regrid var1 to be the same resolution as var2?

Note: After regridding, I want to subsample var1 given some condition related to var2. For example:

var1_subset = var1.where(var2>1000)

So I want to minimize any loss of grid points during the interpolation.

Upvotes: 1

Views: 954

Answers (1)

suvy
suvy

Reputation: 693

basemap.interp will work only when xout and yout are same in number or number of output nlons and nlats are same,

why not generate same length output nlats and nlons and subset it later.

For example:

import numpy as np
import xarray as xray 
import mpl_toolkits.basemap
var1=xray.DataArray(np.random.randn(len(np.linspace(40.5,49.5,10)),len(np.linspace(-39.5,-20.5,20))),coords=[np.linspace(40.5,49.5,10), np.linspace(-39.5,-20.5,20)],dims=['lat','lon'])

(xlon,xlat)=np.meshgrid(np.linspace(-39.875,20.125,80),np.linspace(40.125,49.875,80))
var2=xray.DataArray(-xlon**2+xlat**2,coords[np.linspace(40.125,49.875,80),np.linspace(-39.875,-20.125,80)],dims=['lat','lon'])
mpl_toolkits.basemap.interp(var1,var1.lon,var1.lat,var2.lon,var2.lat,checkbounds=False,masked=False,order=0)

Here is another cool trick with xarray. 

lonreg=var1.groupby_bins('lon',np.linspace(-39.875,20.125,80)).mean(dim='lon')

regridded=lonreg.groupby_bins('lat',np.linspace(-39.5,20.5,20)).mean(dim='lat')

if you want weighted averaged regridding, it is easy to extend this for area averaged regridding by using weights and sum function on groupby object.

Upvotes: 1

Related Questions