python - Half of the world masked when using maskoceans in Basemap

Question

I would like to mask oceans when plotting the data from a netCDF dataset. I followed the great instructions given in the answer to this question. It works great for half of the world, but somehow, everything west of Greenwich is masked as well, both ocean and land. Here is my code:

import netCDF4
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import mpl_toolkits
from mpl_toolkits import basemap
from mpl_toolkits.basemap import Basemap, maskoceans

filename = 'myfile.nc'
vmin = 0.
vmax = 1

nc = netCDF4.Dataset(filename, 'r')
data = nc.variables['sum'][:]  
lats_1d = nc.variables['lat'][:]
lons_1d = nc.variables['lon'][:]
lons, lats = np.meshgrid(lons_1d, lats_1d)

labels = ['DJF', 'MAM', 'JJA', 'SON']
cmap = cm.RdYlBu
cmap.set_over('#00FF00')

my_dpi = 96

fig = plt.figure(figsize=(1200/my_dpi, 800./my_dpi))
for season in range(4):
    ax = fig.add_subplot(2, 2, season+1)
    map1 = basemap.Basemap(resolution='c', projection='kav7', lon_0=0)
    map1.drawcoastlines()
    map1.drawcountries()

    nc_new = maskoceans(lons,lats,data[season,:,:],resolution='c', grid = 1.25)
    datapc = map1.pcolormesh(lons, lats, nc_new,  vmin=vmin, vmax=vmax, cmap=cmap, latlon=True)

    plt.title(labels[season])

fig.tight_layout(pad=1, w_pad=1, h_pad=4)
ax = fig.add_axes([0.05, 0.52, 0.9, 0.025])
cb = plt.colorbar(cax=ax, orientation='horizontal', cmap=cmap,
                  extend='max', format="%.2f",
                  ticks=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1])

plt.show()

I know that a somewhat similar issue was raised here but never got answered, and it appears that in the end, the problem was mixing up lat-long coordinates with x-y ones. I tried switching to x-y coordinates but got the same half-map. Any idea of what can be happening here?

N.B. when plotting the unmasked data using datapc = map1.pcolormesh(lons, lats, data[season,:,:], vmin=vmin, vmax=vmax, cmap=cmap, latlon=True) the whole world is plotted (land + oceans).

Answer 1

As you've identified, the points with longitudes -180 to 0 are not being plotted. Assuming they're in your data, they must be being masked or discarded for some reason.

My intuition was that the dataset longitudes ran 0-360 instead of -180 to 180, which was confirmed in the comments.

The quick fix for this is to add

lons_1d[lons_1d>180]-=360

just after you pull out lons_1d from nc. This works because lons_1d is a numpy array and it uses numpy boolean array indexing (often called "fancy" indexing) to conditionally select the longitude values greater than 180 and subtract 360 from them.

As you note that the pcolormesh plot works if you omit the mask, this looks like a bug with wrapping in the maskoceans function, or at least unexpected behaviour.

For reference - I do not think you are the first to experience similar "wrapping" type issues with masks, I think this issue on the matplotlib github looks rather similar.