How could I plot an array with conditions over a contour map?

I have plotted a global map of GPP using the code below:

( 'lon' and 'lat' are both netCDF4 attributes and have a shape of (144, ) and (90, ) respectively, whilst 'gpp_avg' is a numpy array with a shape of (90, 144) )

import numpy as np
import netCDF4 as n4
import matplotlib.pyplot as plt

import cartopy as cart
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from mpl_toolkits.basemap import Basemap

>> gpp_avg = n4.Dataset('decadal_gpp.nc', 'r')
>> lon = gpp_avg.variables['lon'] # 144 grid cells every 2.5 degrees (east-west)
>> lat = gpp_avg.variables['lat'] # 90 grid cells every 2 degrees (north-south)

>> # Plotting data on a map with Cartopy
>> plt.figure()
>> ax = plt.axes(projection=ccrs.PlateCarree())
>> ax.coastlines() # Adding coastlines
>> ax.add_feature(cart.feature.OCEAN, zorder=100, edgecolor='k')
>> cs = ax.contourf(lon[:], lat[:], gpp_avg[:], cmap = 'Spectral')

>> cbar = plt.colorbar(cs, ax=ax) # Additional necessary information
>> cbar.set_label('g[C]/m^2/day')
>> gridl = ax.gridlines(color="black", linestyle="dotted", 
           draw_labels=True) # Adding axis labels - latitude & longitude
>> gridl.xformatter=LONGITUDE_FORMATTER
>> gridl.yformatter=LATITUDE_FORMATTER
>> gridl.xlabels_top = False
>> gridl.ylabels_right = False
>> plt.show()

I have a numpy array 'ci_95_gpp' which has the shape (90, 144) which contains the p-values for each grid cell of the global map. I want to plot points on top of the global contour map where the p-values are greater than 2.

How would I go about doing this? Many thanks.

I generate a set of data for contour plot on a Cartopy map. The data points for contouring are separated into 2 groups, with negative and positive z-values. Numpy maskedarray is used in that operation. I hope that this is useful for the general readers, including the OP.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import pandas as pd

from numpy.random import uniform, seed
from matplotlib.mlab import griddata   
# TODO, use newer scipy.interpolate() instead of `griddata`
import numpy.ma as ma

# make up some data around long,lat: (90, 18)
seed(0)
npts = 200
x0, y0 = 90, 18   # center of map in (long, lat), degrees
x = x0+uniform(-2, 2, npts)
y = y0+uniform(-2, 2, npts)
#z = x*np.exp(-x**2 - y**2)
z = (x-x0)*np.exp(-(x-x0)**2 - (y-y0)**2)  # elevation in meters

# define grid, for points interpolation from the made-up data above
gridx, gridy = 50,50
xi = x0+np.linspace(-2.1, 2.1, gridx)
yi = y0+np.linspace(-2.1, 2.1, gridy)

# interpolate for gridded data of (gridx, gridy) 
zi = griddata(x, y, z, xi, yi, interp='linear')
# xi.shape, yi.shape, zi.shape  => ((50,), (50,), (50, 50))

xig,yig = np.meshgrid(xi, yi)

# projection
useproj = ccrs.PlateCarree()
fig = plt.figure(figsize = (9, 7))
rect = [0.05, 0.05, 0.95, 0.95]  # for map extent
ax = fig.add_axes( rect, projection=useproj )

# contour the gridded data, plotting dots at the nonuniform data points.
CS = ax.contour(xig, yig, zi, 15, linewidths=0.5, colors='k')
CS = ax.contourf(xig, yig, zi, 15,
                  vmax=abs(zi).max(), vmin=-abs(zi).max())
plt.colorbar(CS)  # draw colorbar

# prep points for scatterplot of the gridded points
# make 2 masked-arrays, based on `zi`
mag = ma.masked_greater(zi, 0)  # mask points with +ve zi values
mal = ma.masked_less(zi, 0)     # mask points with -ve zi values

# apply masking to xig,yig; borrowing mask from mag
xig_greater_masked = ma.MaskedArray(xig, mask=mag.mask)  # must have compatible values
yig_greater_masked = ma.MaskedArray(yig, mask=mag.mask)

# apply masking to xig,yig; borrowing mask from mal
xig_less_masked = ma.MaskedArray(xig, mask=mal.mask)
yig_less_masked = ma.MaskedArray(yig, mask=mal.mask)

# for points with -ve z values (result of .masked_greater)
plt.scatter(xig_greater_masked, yig_greater_masked, s=3, color="w", \
            alpha=1, zorder=15, label="masked_greater z")
# for points with +ve z values (result of .masked_less)
ax.scatter(xig_less_masked, yig_less_masked, s=3, color="r", alpha=1, \
           zorder=15, label="masked_less z") 

leg = ax.legend(title='Masked z', framealpha=1.0, facecolor="lightgray")
leg.set_zorder(20)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle='--')

gl.xlabels_top = False
gl.ylabels_right = False

gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'size': 15, 'color': 'gray'}
#gl.xlabel_style = {'color': 'gray', 'weight': 'bold'}

plt.title('Masked data plot on contour')
plt.show()

The resulting plot: