0

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/file.nc'
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):-

/__init__.py", 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 https://www.dropbox.com/s/ddlpbw5vvj5kz5w/mslp.txt?dl=0

link to lats https://www.dropbox.com/s/lpwjavxtwtt3r13/xlat.txt?dl=0

and the link to lons https://www.dropbox.com/s/1a0q49drfcd2o9h/xlon.txt?dl=0

Community
  • 1
  • 1
sundar_ima
  • 3,604
  • 8
  • 33
  • 52
  • I am not a `netcdf` expert, but I don't see how the code you provided can be run with the plain text file you uploaded. – lanery May 14 '16 at 07:15
  • You can read data as numpy array with this `data = np.loadtxt('mslp.txt',delimiter=',')`. – sundar_ima May 14 '16 at 07:17
  • Sure, but the shape of `data` is (189, 309) and it is not obvious to me what `lons` and `lats` should be. – lanery May 14 '16 at 07:22

1 Answers1

2

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])
m.drawcoastlines()

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

enter image description here

lanery
  • 5,222
  • 3
  • 29
  • 43
  • I have updated the answer with link to data, lats and lons . – sundar_ima May 14 '16 at 08:59
  • Thanks for uploading your data, by using it I think I was able to find the source of your problem. Hopefully my updated answer makes sense and is what you were looking for in a solution. – lanery May 14 '16 at 11:43
  • Thanks for the code. When I try to plot the above code in base map, it runs without an error but no contour lines are visible. Can you plot the same over basemap instead of subplot? – sundar_ima May 14 '16 at 13:55
  • Got it working. But the resolution does not change when I make it coarser or denser. I think it is the problem with converting 2d to 1d array (lats and lons). – sundar_ima May 14 '16 at 15:09
  • Sorry, but I don't have the same issue. To change the resolution you are changing `[::4]` to `[::6]` or `[::10]` for example, yes? – lanery May 14 '16 at 17:37
  • Perhaps when plotting with contour, the change in resolution is just not appearing the way you expect it to. I added plotting with imshow to make it more obvious that the resolution is indeed changing. As further evidence here are plots with [`[::1]`](http://i.stack.imgur.com/LYLOq.png) and [`[::10]`](http://i.stack.imgur.com/KG6Et.png) – lanery May 14 '16 at 18:30
  • @sundar_ima Are you still unsatisfied with the resolution of the contours? Maybe now that I plotted with basemap you can point out what you'd like to be different. – lanery May 17 '16 at 05:34