pv_ghi_1.py

import numpy as np

import matplotlib.pyplot as plt

import cartopy.feature as cfeature
import cartopy.crs as ccrs

import rasterio
from rasterio import plot

cmsaf = rasterio.open('IRRData/gh_0_year.asc')
sarah = rasterio.open('IRRData/gh_0_year_sarah.asc')
sarah2 = rasterio.open('IRRData/gh_0_year_sarah2.tif')

kwh_to_w = 1000/8760
vmin, vmax = 800, 2500
cmap ='YlOrRd'
crs = ccrs.PlateCarree()

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

cmsaf_extent = [cmsaf.bounds[0], cmsaf.bounds[2], cmsaf.bounds[1], cmsaf.bounds[3]]
rasterio.plot.show(cmsaf, title='CM SAF Solar Radiation - PVGIS', ax=axes[0], cmap=cmap, vmin=vmin*kwh_to_w, vmax=vmax*kwh_to_w, extent=cmsaf_extent, zorder=2)

sarah_extent = [sarah.bounds[0], sarah.bounds[2], sarah.bounds[1], sarah.bounds[3]]
rasterio.plot.show(sarah, title='SARAH Solar Radiation - PVGIS', ax=axes[1], cmap=cmap, vmin=vmin*kwh_to_w, vmax=vmax*kwh_to_w, extent=sarah_extent, zorder=2)

sarah2_extent = [sarah2.bounds[0], sarah2.bounds[2], sarah2.bounds[1], sarah2.bounds[3]]
rasterio.plot.show(sarah2, title='SARAH-2 Solar Radiation - PVGIS', ax=axes[2], cmap=cmap, vmin=vmin, vmax=vmax, extent=sarah2_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()