Python (Cartopy) draw shaded figure inside specific country

Question

I am trying to plot precipitation over Cameroon, shading only the inside of Cameroon's borders.
All other countries will be masked.
Below is the plot and the beginning of the script that I am using:

import numpy as np
from scipy import stats 
from netCDF4 import Dataset
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy
import cartopy.io.shapereader as shapereader
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.font_manager import FontProperties
    

ref_dat  = Dataset("R95ptot_Kmer_chirps-v2.0.1981_2019.days_p05.nc", mode='r')
lat     = ref_dat.variables["latitude"][:]
lon     = ref_dat.variables["longitude"][:]
time     = ref_dat.variables["time"]
pre0     = ref_dat.variables["precip"][:, :, :]    
        
slope      =  np.zeros(pre0[0, :, :].shape)
intercept  =  np.zeros(pre0[0, :, :].shape)
r_value    =  np.zeros(pre0[0, :, :].shape)
p_value    =  np.zeros(pre0[0, :, :].shape)
std_err    =  np.zeros(pre0[0, :, :].shape)

for ilat in range(len(lat)):
    for jlon in range(len(lon)):
         slope[ilat, jlon], intercept[ilat, jlon], r_value[ilat, jlon], p_value[ilat, jlon], std_err[ilat, jlon] = stats.linregress(time, pre0[:, ilat, jlon])
         

[lons2d, lats2d] = np.meshgrid(lon, lat) 

fig, ax = plt.subplots(1, 1, sharex=True, sharey=True, figsize=(8.5, 6.98),
subplot_kw={'projection':ccrs.PlateCarree()})

mymap = ax.contourf(lons2d, lats2d, slope, 
                              transform=ccrs.PlateCarree(), 
                              cmap=plt.cm.Spectral_r, extend='both')

Answer 1

Finally i solved my problem. Here plot I am looking for and the script used.

Here a helpful link that used to solved my problem

Python cartopy map, clip area outside country (polygon)

enter image description here

import numpy as np
from scipy import stats 
from netCDF4 import Dataset
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy
import cartopy.io.shapereader as shapereader
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.font_manager import FontProperties
import geopandas
from matplotlib.ticker import (AutoMinorLocator, MultipleLocator)
import matplotlib.ticker as mticker
from shapely.geometry import Polygon

from cartopy.io import shapereader
import cartopy.io.img_tiles as cimgt
import geopandas


def rect_from_bound(xmin, xmax, ymin, ymax):
    """Returns list of (x,y)'s for a rectangle"""
    xs = [xmax, xmin, xmin, xmax, xmax]
    ys = [ymax, ymax, ymin, ymin, ymax]
    return [(x, y) for x, y in zip(xs, ys)]



ref_dat  = Dataset("R95ptot_Kmer_chirps-v2.0.1981_2019.days_p05.nc", mode='r')
lat     = ref_dat.variables["latitude"][:]
lon     = ref_dat.variables["longitude"][:]
time     = ref_dat.variables["time"]
pre0     = ref_dat.variables["precip"][:, :, :]





# request data for use by geopandas
resolution = '10m'
category = 'cultural'
name = 'admin_0_countries'

shpfilename = shapereader.natural_earth(resolution, category, name)
df = geopandas.read_file(shpfilename)

# get geometry of a country
poly = [df.loc[df['ADMIN'] == 'Cameroon']['geometry'].values[0]]

#stamen_terrain = cimgt.Stamen('terrain-background')

# projections that involved
st_proj = ccrs.PlateCarree()#stamen_terrain.crs  #projection used by Stamen images
ll_proj = ccrs.PlateCarree()  #CRS for raw long/lat

# create fig and axes using intended projection
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(1, 1, 1, projection=st_proj)
ax.add_geometries(poly, crs=ll_proj, facecolor='none', edgecolor='black')

pad1 = .1  #padding, degrees unit

exts = [poly[0].bounds[0] - pad1, poly[0].bounds[2] + pad1, poly[0].bounds[1] - pad1, poly[0].bounds[3] + pad1];

#exts = [min(lon), max(lon), min(lat), max(lat)]

ax.set_extent(exts, crs=ccrs.PlateCarree())

ax.set_extent(exts, crs=ll_proj)

# make a mask polygon by polygon's difference operation
# base polygon is a rectangle, another polygon is simplified switzerland
msk = Polygon(rect_from_bound(*exts)).difference( poly[0].simplify(0.01) )
msk_stm  = st_proj.project_geometry (msk, ll_proj)  # project geometry to the projection used by stamen

# get and plot Stamen images
#ax.add_image(stamen_terrain, 8) # this requests image, and plot




        
slope      =  np.zeros(pre0[0, :, :].shape)
intercept  =  np.zeros(pre0[0, :, :].shape)
r_value    =  np.zeros(pre0[0, :, :].shape)
p_value    =  np.zeros(pre0[0, :, :].shape)
std_err    =  np.zeros(pre0[0, :, :].shape)



for ilat in range(len(lat)):
    for jlon in range(len(lon)):
         slope[ilat, jlon], intercept[ilat, jlon], r_value[ilat, jlon], p_value[ilat, jlon], std_err[ilat, jlon] = stats.linregress(time, pre0[:, ilat, jlon])
         

[lons2d, lats2d] = np.meshgrid(lon, lat) 


mymap = ax.contourf(lons2d, lats2d, slope, 
                              transform=ccrs.PlateCarree(), 
                              cmap=plt.cm.Spectral_r, extend='both')




#ax.set_extent([min(lon), max(lon), min(lat), max(lat)], crs=ccrs.PlateCarree())


#cbarax = fig.add_axes([0.2, 0.95, 0.6, 0.02])
cbar = plt.colorbar(mymap, orientation='vertical')
cbar.set_label('Precipitation [mm]', rotation=90, fontsize=12, labelpad=1)
cbar.ax.tick_params(labelsize=12, length=0)


pfieldm = np.ma.masked_greater(p_value, 0.05)

ax.contourf(lons2d, lats2d, pfieldm, transform = ccrs.PlateCarree(), hatches=["..."], alpha=0.0)

ax.add_geometries( msk_stm, st_proj, zorder=12, facecolor='white', edgecolor='none', alpha=1.)

ax.outline_patch.set_edgecolor('white')

#ax.outline_patch.set_visible(False)


plt.savefig('figure.pdf', format='pdf', dpi=500)
 

plt.show()