matplotlib / basemap

Plot on map projections (with coastlines and political boundaries) using matplotlib
MIT License
780 stars 392 forks source link

contourf on cylindrical projection with irregular latlon grid (tri=True) is calling shiftdata when it shouldn't #492

Open Kuhron opened 4 years ago

Kuhron commented 4 years ago

I ran into this problem while trying to plot interpolated data NOT on a regular latlon grid. Orthographic projection works but cylindrical does not.

Relevant code:

fig = plt.figure()
ax = fig.add_subplot(111)
levels = 20
m = Basemap(projection="cyl")  # causes ValueError
# m = Basemap(projection="ortho", lat_0=0., lon_0=0., resolution='l')  # works, does not cause ValueError
MC = m.contourf(lons_deg, lats_deg, vals, levels, ax=ax, tri=True, latlon=True)  # latlon=True interprets first two args as lon and lat respectively

Note that lats_deg, lons_deg, and vals are 1d np arrays with irregular point spacing (no meshgrid). I used tri=True to allow interpolation. This works fine on orthographic projection. But with cylindrical I get the following error:

Traceback (most recent call last):
  File "IcosahedralGeodesicLattice.py", line 146, in <module>
    test_lattice.plot_data(data)
  File "/home/wesley/programming/Mapping/Lattice.py", line 86, in plot_data
    MC = m.contourf(lons_deg, lats_deg, vals, levels, ax=ax, tri=True, latlon=True)  # latlon=True interprets first two args as LON and LAT RESPECTIVELY
  File "/home/wesley/.local/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py", line 544, in with_transform
    fix_wrap_around=plotfunc.__name__ not in ["scatter"])
  File "/home/wesley/.local/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py", line 4797, in shiftdata
    if fix_wrap_around and itemindex:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Looking at the relevant file, the docstring for shiftdata says "Only valid for cylindrical/pseudo-cylindrical global projections and data on regular lat/lon grids. longitudes and data can be 1-d or 2-d, if 2-d it is assumed longitudes are 2nd (rightmost) dimension."

So should it not be being called with tri=True? Sorry if I have misunderstood something, let me know. Thanks!

chengdang commented 3 years ago

@Kuhron Hey, have you found any solution to this issue? I'm having the same problem.

Kuhron commented 3 years ago

@chengdang Nope, I just stopped using Basemap altogether. Instead I wrote my own code to do the triangulation for equirectangular projection. Here is the relevant code, if you want to use parts of it:

import matplotlib.pyplot as plt
import matplotlib.tri as tri  # interpolation of irregularly spaced data

x, y = lons_deg, lats_deg  # actual data
z = vals
# how many x/y to put on rectangle projection
n_grid_x = 1000
n_grid_y = 500
xi = np.linspace(-180, 180, n_grid_x)
yi = np.linspace(-90, 90, n_grid_y)
triang = tri.Triangulation(x, y)
interpolator = tri.LinearTriInterpolator(triang, z)
Xi, Yi = np.meshgrid(xi, yi)
zi = interpolator(Xi, Yi)
MC = plt.contourf(xi, yi, zi, levels=contourf_levels, cmap=cmap)
clb = plt.colorbar(MC, ax=plt.gca())  # without these args, it will say it can't find a mappable object for colorbar
if contour_lines:
    plt.contour(xi, yi, zi, levels=contour_line_levels, linewidths=0.5, colors='k')
clb.ax.set_title(key_str)