pv_ghi_2.py

cams_eu = rasterio.open('IRRData/Europe_Average_GHI_yearly_gf.tif')
cams_af = rasterio.open('IRRData/Africa_Average_GHI_yearly_gf.tif')

# Necessary to mask zeros, as 'nodata' attribute is set incorrectly in the tif file
cams_eu_masked = np.ma.masked_where(cams_eu.read(1)==0, cams_eu.read(1))
cams_af_masked = np.ma.masked_where(cams_af.read(1)==0, cams_af.read(1))

vmin, vmax = 800, 2500
cmap ='YlOrRd'
crs = ccrs.PlateCarree()

fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(15,10), sharex=True, subplot_kw={'projection': crs})

cams_ue_extent = [cams_eu.bounds[0], cams_eu.bounds[2], cams_eu.bounds[1], cams_eu.bounds[3]]
rasterio.plot.show(cams_eu_masked, title='CAMS Radiation Europe', ax=axes[0], cmap=cmap, vmin=vmin, vmax=vmax, extent=cams_ue_extent, zorder=2)

cams_af_extent = [cams_af.bounds[0], cams_af.bounds[2], cams_af.bounds[1], cams_af.bounds[3]]
rasterio.plot.show(cams_af_masked, title='CAMS Radiation Africa', ax=axes[1], cmap=cmap, vmin=vmin, vmax=vmax, extent=cams_af_extent, zorder=2)

for i, ax in enumerate(axes.flatten()):
    ax.add_feature(cfeature.LAND, facecolor='lightgrey', zorder=0)
    ax.add_feature(cfeature.OCEAN, linewidth=0, facecolor='white', zorder=2)
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5, zorder=3)
    ax.set_extent([-180, 180, -90, 90])

plt.show()