Plotting oceans in maps using basemap and python

Question

I am plotting the netCDF file available here: https://goo.gl/QyUI4J

Using the code below, the map looks like this:

However, I want the oceans to be in white color. Better still, I want to be able to specify what color the oceans show up in. How do I change the code below to do that? Right now, the issue is that the oceans are getting plotted on the data scale. (please note that the netCDF file is huge ~3.5 GB).

import pdb, os, glob, netCDF4, numpy
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap

def plot_map(path_nc, var_name):
    """
    Plot var_name variable from netCDF file
    :param path_nc: Name of netCDF file
    :param var_name: Name of variable in netCDF file to plot on map
    :return: Nothing, side-effect: plot an image
    """

    nc = netCDF4.Dataset(path_nc, 'r', format='NETCDF4')
    tmax  = nc.variables['time'][:]

    m = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0)

    m.drawcoastlines()
    m.drawcountries()

    # find x,y of map projection grid.
    lons, lats = get_latlon_data(path_nc)
    lons, lats = numpy.meshgrid(lons, lats)
    x, y = m(lons, lats)
    nc_vars = numpy.array(nc.variables[var_name])

    # Plot!
    m.drawlsmask(land_color='white',ocean_color='white')
    cs = m.contourf(x,y,nc_vars[len(tmax)-1,:,:],numpy.arange(0.0,1.0,0.1),cmap=plt.cm.RdBu)

    # add colorbar
    cb = m.colorbar(cs,"bottom", size="5%", pad='2%')
    cb.set_label('Land cover percentage '+var_name+' in '+os.path.basename(path_nc))

    plt.show()


plot_map('perc_crops.nc','LU_Corn.nc')

Sematic

what is get_latlon_data ? – Gandy

Answer 1

You need to use maskoceans on your nc_vars dataset

Before contourf, insert this

nc_new = maskoceans(lons,lats,nc_vars[len(tmax)-1,:,:])

and then call contourf with the newly masked dataset i.e.

cs = m.contourf(x,y,nc_new,numpy.arange(0.0,1.0,0.1),cmap=plt.cm.RdBu)

To dictate the ocean colour, you can either drop the call to drawslmask if you want white oceans or specify an ocean colour in that call - e.g. insert m.drawlsmask(land_color='white',ocean_color='cyan').

I've given the working code with as few alterations to yours as possible below. Uncomment the call to drawslmask to see cyan coloured oceans.

Output

Full working version of code

import pdb, os, glob, netCDF4, numpy
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap, maskoceans

def plot_map(path_nc, var_name):
    """
    Plot var_name variable from netCDF file
    :param path_nc: Name of netCDF file
    :param var_name: Name of variable in netCDF file to plot on map
    :return: Nothing, side-effect: plot an image
    """

    nc = netCDF4.Dataset(path_nc, 'r', format='NETCDF4')
    tmax  = nc.variables['time'][:]

    m = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0)

    m.drawcoastlines()
    m.drawcountries()    

    # find x,y of map projection grid.
    lons, lats = nc.variables['lon'][:],nc.variables['lat'][:]
    # N.B. I had to substitute the above for unknown function get_latlon_data(path_nc)
    # I guess it does the same job

    lons, lats = numpy.meshgrid(lons, lats)
    x, y = m(lons, lats)
    nc_vars = numpy.array(nc.variables[var_name])
    
    #mask the oceans in your dataset
    nc_new = maskoceans(lons,lats,nc_vars[len(tmax)-1,:,:])
    
    #plot!
    #optionally give the oceans a colour with the line below
    #Note - if land_color is omitted it will default to grey
    #m.drawlsmask(land_color='white',ocean_color='cyan')
    cs = m.contourf(x,y,nc_new,numpy.arange(0.0,1.0,0.1),cmap=plt.cm.RdBu)

    # add colorbar
    cb = m.colorbar(cs,"bottom", size="5%", pad='2%')
    cb.set_label('Land cover percentage '+var_name+' in '+os.path.basename(path_nc))

    plt.show()


plot_map('perc_crops.nc','LU_Corn.nc')

---

P.S. That's a big file to test!!

Minivet

Answer 2

The lawful good solution is to use the utility function maskoceans, which takes in a data array and masks all the points in oceans and lakes.

Instead, you could take the easy way out. Draw your contour plot first, then use drawlsmask, which allows transparent colors:

# Colors can be RGBA tuples
m.drawlsmask(land_color=(0, 0, 0, 0), ocean_color='deeppink', lakes=True)

Land is transparent, which lets the contour plot shows through.

Shanel

Answer 3

The colors that youbsee in the map are related to the colormap cm.plt.RdBu which is passed to the contourcf function. You need to change this color map to achieve the results that you want. Here you can find a tutorial for the basemap colormap.

Impolitic

thanks @rflkortekaas, this is not related to the colormap. The problem is that oceans have value 0.0 or less in this netCDF. When I plot the data, they get the same color as the land pixels with 0.0 values. I want oceans to be drawn separately (as white or some other color of my choice) – Sematic

It's related to the colormap but your information above was not clear. There is an option on the basemap to mask oceans in the array, maybe you can achieve your requirement with that. – Impolitic

I am setting ocean_color to white in drawlsmask, but it does not seem to work – Sematic

That doesn't work due to the fact that contourf overlays the drawlsmask. So your contourf function draws the oceans with a red color from the colour map and draws that over the white ocean from drawlsmask. – Impolitic

true, but moving drawslmask below contourf doesn't help either. – Sematic

I understand that. I doesn't have time to test it but you definitely need to look into the maskoceans function of the basemap. That way you can define the colot with drawlsmask and the oceans in the data array are masked and not drawn by the contourf function. – Impolitic