sundar_ima sundar_ima - 1 year ago 327
Python Question

Regridding high resolution gridded data to low resolution coarser resolution data?

I have a netcdf file with high resolution grided data (.1 x .1 deg) and plotting on the basemap using matplotlib. The plot I am trying to do is contour lines. For a specific reason I would like to plot the data at the interval of 1 x 1 deg resolution. For which I have used the following code sample example from here Regridding regular netcdf data.

See the update 1 for link to actual data.

For the sake of clearity following is the code I am trying to execute for regridding to lower resolution:-

from mpl_toolkits import basemap
from netCDF4 import Dataset

filename = 'path/to/netcdf/'
with Dataset(filename, mode='r') as fh:
lons = fh.variables['LON'][:]
lats = fh.variables['LAT'][:]
data = fh.variables['Data'][:].squeeze()

lons_sub, lats_sub = np.meshgrid(lons[::4], lats[::4], sparse=True)

data_coarse = basemap.interp(data, lons, lats, lons_sub, lats_sub, order=1)

The code seems to be correct. But when I execute the code I get the following error in the line
data_coarse = basemap.interp(data, lons, lats, lons_sub, lats_sub, order=3)

/", line 4897, in interp
if xin[-1]-xin[0] < 0 or yin[-1]-yin[0] < 0:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I have not understood where the problem lies and how to solve this?

Any help is appreciated.

Update 1

link to the actual data

link to lats

and the link to lons

Answer Source

In your code snippet, xin = lons and yin = lats. From the basemap.interp docs,

xin, yin: rank-1 arrays containing x and y of datain grid in increasing order.

meaning xin (lons) and yin (lats) have to be continuously increasing. From your error traceback, it looks like that is not the case since either xin[-1]-xin[0] < 0 and/or yin[-1]-yin[0] < 0.

Without knowing what lons or lats are exactly (since your code does not run with the data you link to), it is tough to elaborate any further.

Edit - plotting with basemap

I think the main problem you were running into was the shape of your latitude and longitude data. Since lons and lats were essentially stacks of the same arrays, I converted them both to 1d arrays as you'll see below.

import numpy as np
from mpl_toolkits import basemap
import matplotlib.pyplot as plt

# load data
data = np.loadtxt('mslp.txt', delimiter=',')
lons = np.loadtxt('xlon.txt', delimiter=',')
lats = np.loadtxt('xlat.txt', delimiter=',')

# take 1 dimensional slice of lons and lats
longitude = lons[0]   # array([31.1, 31.4, 31.7,...,119.1, 119.3, 119.6])
latitude  = lats[:,0] # array([-2.3, -2.1, -1.8,...,39.4, 39.6, 39.8])

# subdivide longitude and latitude into ~1deg grid data
lons_sub, lats_sub = np.meshgrid(longitude[::4], latitude[::4])

# implement basemap.interp to interpolate along ~1deg grid data
data_course = basemap.interp(datain=data, xin=longitude, yin=latitude,
                             xout=lons_sub, yout=lats_sub, order=3)

m = basemap.Basemap(llcrnrlon=longitude[0],  llcrnrlat=latitude[0],
                    urcrnrlon=longitude[-1], urcrnrlat=latitude[-1],
                    projection='merc', resolution='c')

fig = plt.figure()
# ax1 -- contour plot of resampled data
ax1 = fig.add_subplot(211)
m.contour(lons_sub, lats_sub, data_course, linewidths=2.5, latlon=True)
m.fillcontinents(color=[0.7, 0.9, 0.8], lake_color=[0.3, 0.7, 0.9])
m.drawmapboundary(fill_color=[0.3, 0.7, 0.9])

# ax2 -- imshow plot of resampled data
ax2 = fig.add_subplot(212)
m.imshow(data_course, interpolation='nearest',
         extent=[lons_sub[0], lons_sub[-1], lats[0], lats[-1]])

enter image description here