<jupyter_start><jupyter_text>February and April 2020 precipitation anomaliesIn this notebook, we will analyze precipitation anomalies of February and April 2020, which seemed to be very contrasting in weather. We use the EOBS dataset.   Import packages<jupyter_code>##This is so variables get printed within jupyter
from IPython.core.interactiveshell import InteractiveShell 
InteractiveShell.ast_node_interactivity = "all"
##import packages
import os
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import matplotlib.ticker as mticker

##We want the working directory to be the UNSEEN-open directory
pwd = os.getcwd() ##current working directory is UNSEEN-open/Notebooks/1.Download
pwd #print the present working directory
os.chdir(pwd+'/../../') # Change the working directory to UNSEEN-open
os.getcwd() #print the working directory<jupyter_output><empty_output><jupyter_text>Load EOBSI downloaded EOBS (from 1950 - 2019) and the most recent EOBS data (2020) [here](https://surfobs.climate.copernicus.eu/dataaccess/access_eobs.php). Note, you have to register as E-OBS user.The data has a daily timestep. I resample the data into monthly average mm/day. I chose not to use the total monthly precipitation because of leap days.<jupyter_code>EOBS = xr.open_dataset('../UK_example/EOBS/rr_ens_mean_0.25deg_reg_v20.0e.nc') ## open the data
EOBS = EOBS.resample(time='1m').mean() ## Monthly averages
# EOBS = EOBS.sel(time=EOBS['time.month'] == 2) ## Select only February
EOBS<jupyter_output>/soge-home/users/cenv0732/.conda/envs/UNSEEN-open/lib/python3.8/site-packages/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Here I define the attributes, that xarray uses when plotting<jupyter_code>EOBS['rr'].attrs = {'long_name': 'rainfall',  ##Define the name
 'units': 'mm/day', ## unit
 'standard_name': 'thickness_of_rainfall_amount'} ## original name, not used
EOBS['rr'].mean('time').plot() ## and show the 1950-2019 average February precipitation<jupyter_output><empty_output><jupyter_text>The 2020 data file is separate and needs the same preprocessing:<jupyter_code>EOBS2020 = xr.open_dataset('../UK_example/EOBS/rr_0.25deg_day_2020_grid_ensmean.nc.1') #open
EOBS2020 = EOBS2020.resample(time='1m').mean() #Monthly mean
EOBS2020['rr'].sel(time='2020-04').plot() #show map
EOBS2020 ## display dataset<jupyter_output>/soge-home/users/cenv0732/.conda/envs/UNSEEN-open/lib/python3.8/site-packages/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Plot the 2020 eventI calculate the anomaly (deviation from the mean in mm/d) and divide this by the standard deviation to obtain the standardized anomalies.<jupyter_code>EOBS2020_anomaly = EOBS2020['rr'].groupby('time.month') - EOBS['rr'].groupby('time.month').mean('time')
EOBS2020_anomaly

EOBS2020_sd_anomaly = EOBS2020_anomaly.groupby('time.month') / EOBS['rr'].groupby('time.month').std('time')

EOBS2020_sd_anomaly.attrs = {
    'long_name': 'Monthly precipitation standardized anomaly',
    'units': '-'
}

EOBS2020_sd_anomaly<jupyter_output>/soge-home/users/cenv0732/.conda/envs/UNSEEN-open/lib/python3.8/site-packages/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>I select February and April (tips on how to select this are appreciated)<jupyter_code>EOBS2020_sd_anomaly
# EOBS2020_sd_anomaly.sel(time = ['2020-02','2020-04']) ## Dont know how to select this by label?
EOBS2020_sd_anomaly[[1,3],:,:] ## Dont know how to select this by label?<jupyter_output><empty_output><jupyter_text>And plot using cartopy!<jupyter_code>EOBS_plots = EOBS2020_sd_anomaly[[1, 3], :, :].plot(
    transform=ccrs.PlateCarree(),
    robust=True,
    col='time',
    cmap=plt.cm.twilight_shifted_r,
    subplot_kws={'projection': ccrs.EuroPP()})

for ax in EOBS_plots.axes.flat:
    ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
    ax.coastlines(resolution='50m')
    gl = ax.gridlines(crs=ccrs.PlateCarree(),
                      draw_labels=False,
                      linewidth=1,
                      color='gray',
                      alpha=0.5,
                      linestyle='--')

# plt.savefig('graphs/February_April_2020_precipAnomaly.png', dpi=300)
# EOBS = EOBS.sel(time=EOBS['time.month'] == 2) ## Select only February

EOBS2020_anomaly_f = EOBS2020['rr'].sel(time='2020-02') - EOBS['rr'].sel(time=EOBS['time.month'] == 2).mean('time')
EOBS2020_sd_anomaly_f = EOBS2020_anomaly / EOBS['rr'].sel(time=EOBS['time.month'] == 2).std('time')

EOBS2020_sd_anomaly_f.attrs = {
    'long_name': 'Monthly precipitation standardized anomaly',
    'units': '-'
}


EOBS2020_anomaly_a = EOBS2020['rr'].sel(time='2020-04') - EOBS['rr'].sel(time=EOBS['time.month'] == 4).mean('time')
EOBS2020_sd_anomaly_a = EOBS2020_anomaly_a / EOBS['rr'].sel(time=EOBS['time.month'] == 4).std('time')

EOBS2020_sd_anomaly_a.attrs = {
    'long_name': 'Monthly precipitation standardized anomaly',
    'units': '-'
}

EOBS2020_sd_anomaly_f
EOBS2020_sd_anomaly_a
fig, axs = plt.subplots(1, 2, subplot_kw={'projection': ccrs.EuroPP()})

EOBS2020_sd_anomaly_f.plot(transform=ccrs.PlateCarree(),
                           robust=True,
                           cmap=plt.cm.twilight_shifted_r,
                           ax=axs[0])  #,cmap=plt.cm.Blues,

EOBS2020_sd_anomaly_a.plot(transform=ccrs.PlateCarree(),
                           robust=True,
                           vmin=-3,
                           vmax=3,
                           cmap=plt.cm.twilight_shifted_r,
                           ax=axs[1])  #,cmap=plt.cm.Blues,

for ax in axs.flat:
    ax.coastlines(resolution='50m')

# axs[0].set_title('SEAS5')
# axs[1].set_title('EOBS')
# axs[2].set_title('EOBS 2020')
EOBS['rr'].groupby('time.month').mean('time')
EOBS['rr'].groupby('time.month').mean('time').plot(
    transform=ccrs.PlateCarree(),
    robust=True,
    col='month',
    col_wrap=3,
    subplot_kws={'projection': ccrs.EuroPP()})<jupyter_output>/soge-home/users/cenv0732/.conda/envs/UNSEEN-open/lib/python3.8/site-packages/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)<jupyter_text>Concatenate SEAS5The SEAS5 data consists of 175 netcdf files that we downloaded from CDS ([1.2](reference). We concatenate those into one netcdf file here. We create a netcdf file containing the dimensions lat, lon, number (25 ensembles), time (35 years) and leadtime (5 initialization months). To generate this, we loop over lead times, and open all 35 years of the lead time and then concatenate those leadtimes.<jupyter_code>init_months = np.append(np.arange(9,13),1) ## Initialization months 9-12,1 (Sep-Jan)

SEAS5_ld1= xr.open_mfdataset('../UK_example/SEAS5/' + '*01.nc', combine='by_coords') # Load the first lead time 
SEAS5= SEAS5_ld1 # Create the xarray dataset to concatenate over
for init_month in init_months[0:-1][::-1]: ## Remove the first that we already have and invert so it is lead time 1-5
    print(init_month)
    SEAS5_ld = xr.open_mfdataset('../UK_example/SEAS5/' + '*' +  "%.2i" % init_month + '.nc', combine='by_coords') # Load the first lead time to concatenate over
    SEAS5=xr.concat([SEAS5,SEAS5_ld],dim='leadtime') 

SEAS5 = SEAS5.assign_coords(leadtime = np.arange(5) + 2) # assign leadtime coordinates
SEAS5['tprate'] = SEAS5['tprate'] * 1000 * 3600 * 24 ## From m/s to mm/d
SEAS5['tprate'].attrs = {'long_name': 'rainfall',
 'units': 'mm/day',
 'standard_name': 'thickness_of_rainfall_amount'}
SEAS5.load()
# SEAS5.to_netcdf('../UK_example/SEAS5/SEAS5.nc') # for future import<jupyter_output><empty_output><jupyter_text>Compare SEAS5 to EOBS<jupyter_code>fig, axs = plt.subplots(1, 2,subplot_kw={'projection': ccrs.OSGB()})

SEAS5['tprate'].mean(dim = ['time','leadtime','number']).plot(transform=ccrs.PlateCarree(), vmin = 0, vmax = 8, cmap=plt.cm.Blues, ax=axs[0])  #,cmap=plt.cm.Blues,
EOBS['rr'].mean('time').plot(transform=ccrs.PlateCarree(), vmin = 0, vmax = 8, cmap=plt.cm.Blues, ax=axs[1])
# EOBS2020['rr'].plot(transform=ccrs.PlateCarree(), vmin = 0, vmax = 8, cmap=plt.cm.Blues, ax=axs[2])

for ax in axs.flat:
    ax.coastlines(resolution='110m')
    
axs[0].set_title('SEAS5')
axs[1].set_title('EOBS')
# axs[2].set_title('EOBS 2020')
fig, axs = plt.subplots(1, 2,subplot_kw={'projection': ccrs.OSGB()})

(SEAS5['tprate'].mean(dim = ['time','leadtime','number'])
#  .sel(latitude=slice(),longitude=slice())
 .plot(transform=ccrs.PlateCarree(), vmin = 0, vmax = 8, cmap=plt.cm.Blues, ax=axs[0]))  #,cmap=plt.cm.Blues,
EOBS['rr'].mean('time').plot(transform=ccrs.PlateCarree(), vmin = 0, vmax = 8, cmap=plt.cm.Blues, ax=axs[1])

for ax in axs.flat:
    ax.coastlines(resolution='110m')
    ax.gridlines()
#     ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
#                       linewidth=2, color='gray', alpha=0.5, linestyle='--')
axs[0].set_title('SEAS5')
axs[1].set_title('EOBS')
# np.arange('1950', '2017', dtype='datetime64[Y]')
EOBS['rr'].mean(dim = ['latitude','longitude']).plot(x='time')

Quantiles = SEAS5['tprate'].mean(dim = ['latitude','longitude']).quantile([2.5/100, 0.5, 97.5/100], dim=['number','leadtime'])
Quantiles.plot.line(x='time')
# SEAS5['tprate'].mean(dim = ['latitude','longitude','number','leadtime']).plot.line(x = 'time')
Quantiles = SEAS5['tprate'].mean(dim = ['latitude','longitude']).quantile([2.5/100, 0.5, 97.5/100], dim=['number','leadtime'])
Quantiles.plot.line(x='time')
# plt.plot(Quantiles.sel(quantile=0.5),label = 'SEAS5 Median')
# plt.fill_between(Quantiles.sel(quantile=0.025), Quantiles.sel(quantile=0.975), color='green', alpha=0.2,label = '97.5% percentile')
SEAS5['tprate'].mean(dim = ['latitude','longitude','number']).plot.line(x = 'time')

Quantiles = Tapajos.assign_coords(time=np.tile(time_year,400)).groupby("time.month").quantile([2.5/100, 0.5, 97.5/100, 99/100, 0.999, 0.9999], dim="time")
plt.plot(np.arange(1,13), Quantiles.sel(quantile=0.5),label = 'Median of 2000 years')
plt.fill_between(np.arange(1,13), Quantiles.sel(quantile=0.025), Quantiles.sel(quantile=0.975), color='green', alpha=0.2,label = '97.5% percentile')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16241636769126344, 'fraction_numerical': 0.04819160283137787, 'mean_word_length': 4.930994824611846, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '284174', 'license': 'MIT', 'n_tokens_mistral': 4226, 'n_tokens_neox': 3860, 'n_words': 711, 'path': 'doc/_build/.doctrees/nbsphinx/Notebooks/dev/2020_contrasting_weather-Copy1.ipynb', 'repo_name': 'esowc/UNSEEN-open', 'size': '10313'}