Rainfall

import matplotlib
matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
import matplotlib.pyplot as plt
import matplotlib.pyplot as P
from mpl_toolkits.basemap import Basemap,shiftgrid,maskoceans
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import shapefile
import numpy as N
from netCDF4 import Dataset, num2date,date2num
import datetime
from scipy import signal,stats
import matplotlib.colors as colors
from matplotlib.offsetbox import AnchoredText
import matplotlib as mpl


plt.style.use('seaborn-talk')

params = {      'legend.fontsize': 8,\
                'font.family': 'serif',\
                }
plt.rcParams.update(params)

city = {

   'Amman':[35.930359,31.963158],
   'Aqaba':[35.00778,29.52667],
   'WadiMusa':[35.480125,30.321635],
}


#file = '../../data/ClimateExplorer/v2p0chirps_25_34-40E_29-34N.cdf'
#file = '../../data/ClimateExplorer/tp_daily_uerra_19810101-20181231_rm.cdf'
file = '../../data/ClimateExplorer/era5_tp_daily_af_34-40E_29-34N_su.cdf'
#file = '../../data/ClimateExplorer/pr_W5E5v2.0_19810101-20191231.cdf'
#file = '../../data/ClimateExplorer/chirps-v2.0.1981-2022.days_p25.cdf'


nc = Dataset(file,'r')
lon = N.array(nc.variables['lon'][:])
lat = N.array(nc.variables['lat'][:])
dat = N.array(nc.variables['tp'][:])
tim = nc.variables['time']

dat[dat>1000] = N.nan

tim = num2date(tim[:],units=tim.units,calendar=tim.calendar) 

jj = []
mm = []
dd = []

for it in tim:

    jj.append(it.year)
    mm.append(it.month)
    dd.append(it.day)

jj = N.array(jj)
mm = N.array(mm)
dd = N.array(dd)

nc.close()

nd = dat.shape[0]
ny = dat.shape[1]
nx = dat.shape[2]

for d in range(nd):

    date = '%04i-%02i-%02i'%(jj[d],mm[d],dd[d])

    tmp = N.reshape(dat[d,:,:],(ny*nx))
    rx = N.nanmax(tmp)

    print (date,rx)

    fig = P.figure(figsize=(8,4))
    ax = P.subplot(111)

    m = Basemap(projection='cyl',llcrnrlat=29,urcrnrlat=33.5,llcrnrlon=34.5,urcrnrlon=39.5,resolution='h')

    xx,yy = m(*N.meshgrid(lon,lat))

    lev = [5,10,20,30,40,50,75,100]

    cl = m.contourf(xx,yy,dat[d,:,:],levels=lev,cmap=P.get_cmap('cool'),extend='max',zorder=4)
    cs = m.contour(xx,yy,dat[d,:,:],levels=lev,colors='k',linestyles='solid',zorder=5,linewidths=0.5)
    cb = fig.colorbar(cl,location='right',shrink=0.5,pad=-0.15,anchor=(0.1, 0.1))

    cb.set_label('mm/d',labelpad=-20,y=1.2,rotation=0,fontsize=10,weight='bold')

    at = AnchoredText('ERA5: %s'%date,prop=dict(size=8,weight='bold'),frameon=True,loc='upper left')
    at.patch.set_boxstyle("round,pad=0.,rounding_size=0.1")
    ax.add_artist(at)

    #at = AnchoredText('%.1f mm/d'%(rx),prop=dict(size=8,weight='bold'),frameon=True,loc='lower right')
    #at.patch.set_boxstyle("round,pad=0.,rounding_size=0.1")
    #ax.add_artist(at)

    for c in city:

        x,y = m(city[c][0],city[c][1])

        m.scatter(x,y,s=30,c='r',ec='k',lw=0.5,zorder=9)
        plt.text(x,y,c,fontsize=8,weight='bold',zorder=10)    

    sf = shapefile.Reader("../../data/shp/JOR_adm0")

    for shape_rec in sf.shapeRecords():
            if shape_rec.record[1] == 'JOR':
                vertices = []
                codes = []
                pts = shape_rec.shape.points
                prt = list(shape_rec.shape.parts) + [len(pts)]
                for i in range(len(prt) - 1):
                    for j in range(prt[i], prt[i+1]):
                        vertices.append(m(pts[j][0], pts[j][1]))
                    codes += [Path.MOVETO]
                    codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
                    codes += [Path.CLOSEPOLY]
                clip = Path(vertices, codes)
                clip = PathPatch(clip, transform=ax.transData,fill=False,lw=0.5)

    for contour in cl.collections:
           contour.set_clip_path(clip)
           ax.add_patch(clip)

    for contour in cs.collections:
           contour.set_clip_path(clip)
           ax.add_patch(clip)

    ax.axis('off')

    plt.savefig('../../data/png/era5/%s.png'%(date),dpi=120,transparent=True,bbox_inches='tight',pad_inches=0.0)
    plt.close()