matplotlib / basemap

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

pcolormesh false returns masked array error #553

Closed lexipalmer13 closed 1 year ago

lexipalmer13 commented 1 year ago

Using matplotlib with basemap to plot temperature data and getting an error saying pcolormesh cannot plot masked data. When raised on github (link) response suggested that it's a corner case and arrays get overflowed rather than an issue with the data. Putting here as a reminder as requested - happy to provide more info if helpful!

guziy commented 1 year ago

The problem is that the conic projection is not valid globally, and you are trying to feed it era5 global fields... Probably basemap could be more intelligent about it...

Here is the workaround...

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

data = Dataset('ERA5_1959.nc')

lat = data.variables['latitude'][:]
lon = data.variables['longitude'][:]

lon, lat = np.meshgrid(lon, lat)
t2m = data.variables['t2m'][1,1,:,:] #temp data

#create non-masked arrays
lon2 = lon
lat2 = lat

print(lat2.min(), lat2.max())
print(lon2.min(), lon2.max())

#plot
fig = plt.figure(figsize=(10, 8))
m = Basemap(projection='lcc', resolution='c',
            width=8E6, height=8E6, 
            lat_0=45, lon_0=-100,)
m.shadedrelief(scale=0.5)

i_arr, j_arr = np.where((lat2 >= m.latmin) & (lat2 <= m.latmax))

i_min, i_max = i_arr.min(), i_arr.max()
j_min, j_max = j_arr.min(), j_arr.max()

i_slice = slice(i_min, i_max + 1)
j_slice = slice(j_min, j_max + 1)
t2m = t2m[i_slice, j_slice]
lon2 = lon2[i_slice, j_slice]
lat2 = lat2[i_slice, j_slice]

x, y = m(lon2, lat2)
m.pcolormesh(x, y, t2m, cmap='RdBu_r')
# m.pcolormesh(lon2, lat2, t2m,
#              latlon=True, cmap='RdBu_r')
# plt.clim(-8, 8)
m.drawcoastlines(color='lightgray')

plt.title('January 2014 Temperature Anomaly')
plt.colorbar(label='temperature anomaly (°C)');
plt.show()
guziy commented 1 year ago

I get this warning though, it seems like it does not like something about coordinates...

test-basemap-overflow-coords/test.py:48: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
  m.pcolormesh(x, y, t2m, cmap='RdBu_r')
guziy commented 1 year ago

The warning goes away if I do this:

i_slice = slice(i_min, i_max + 1)
j_slice = slice(j_min, j_max + 1)
t2m = t2m[i_min:i_max, j_min:j_max]
lon2 = lon2[i_slice, j_slice]
lat2 = lat2[i_slice, j_slice]
lexipalmer13 commented 1 year ago

My apologies for not responding to this earlier, but this fix works! Thank you so much for your help - I really appreciate it!