<jupyter_start><jupyter_text>CESM subsampling pi control with different 'ensemble' time periods Author: <NAME>, [github](https://github.com/chrisrwp/synthetic-ensemble/SIA/CESM_pi_control.ipynb)**Input:** CESM Large Ensemble preindustrial control, [website](https://www.cesm.ucar.edu/projects/community-projects/LENS/instructions.htmlESG), [CESM LE doi](https://doi.org/10.1175/BAMS-D-13-00255.1)- Coupled (years 0400-2200): `b.e11.B1850C5CN.f09_g16.005.cice.h.aice_nh*.nc` - Slab-ocean (years 0100-1001): `e.e11.E1850C5CN.f09_g16.001.cice.h.aice_nh*.nc` CESM Large Ensemble, historical and RCP8.5- 1920???-1990?????, little trend during that time (20th century near constant - 20C)- 1995???-2045???, approx linear trend for that time (21th century near linear trend - 21C)**Output:** - SIA time series for coupled and slab-ocean preindustrial controls- Detrended SIA time series- Standard deviation of resampled time series for all time period lengths**Method:** - Use the 40 emsemble members from LE or 'simulate' 40 ensemble members in the preindustrial control by using randomly chosen start dates for 40 time series. **Use the same random start dates for all time period lengths**- Allow the time series to be of length between 6 and 298 years with an interval of 4 years- Detrend either by each time series or from the mean trend from the 40 'ensemble' members for each time period- Resample 1000 times with a 2 year block boostrap size - Repeat the analysis for the preindustrial control 10 times to select different randomly chosen start years and average the results**TO DO:** - Decide on start/end dates for 20C and 21C. Make reduced dataset files for those. Should be easy to resample these and detrend based on ensemble or individual trend - Redo the code to select the same random start years for all time periods- Also look at Atmosphere-only control run, f.e11.F1850C5CN.f09_f09.001, years 1-2600 availible?<jupyter_code>%matplotlib inline
# import matplotlib.pyplot as plt
# import matplotlib.path as mpath
# from matplotlib.patches import Patch
import numpy as np
import xarray as xr
import datetime
import dask
# from cftime import num2date, DatetimeNoLeap
print(datetime.datetime.utcnow().strftime("%H:%M UTC %a %Y-%m-%d"))

data_path = '/glade/scratch/cwpowell/Synthetic_ensemble/'
from dask_jobqueue import PBSCluster
from dask.distributed import Client

cluster = PBSCluster(cores    = 8,
                     memory   = '4GB',
                     queue    = 'regular',
                     walltime = '00:44:00')

cluster.scale(8)
client = Client(cluster)
client<jupyter_output><empty_output><jupyter_text>CESM LE 1920-1990 and 1985-2045 Make a reduced dataset for the 20C and 21C time periods of the CESM LE   Resample the two eras with 6-70 or 6-60 year length time series, interval of 2 years<jupyter_code>#compare with existing data from the 1979-2020 resampling<jupyter_output><empty_output><jupyter_text>CESM Preindustrial Control Make reduced datast of SIA for whole of PI control (slab-ocean and coupled)<jupyter_code>#open the pi_control files
pi_control = xr.open_mfdataset('/glade/collections/cdg/data/cesmLE/CESM-CAM5-BGC-LE/ice/proc/tseries/monthly/aice/b.e11.B1850C5CN.f09_g16.005.cice.h.aice_nh*', concat_dim='time', chunks=104)
pi_control_slab = xr.open_mfdataset('/glade/collections/cdg/data/cesmLE/CESM-CAM5-BGC-LE/ice/proc/tseries/monthly/aice/e.e11.E1850C5CN.f09_g16.001.cice.h.aice_nh*', concat_dim='time', chunks=104)
#open the corresponding area file
CESM_area_file = xr.open_dataset('/glade/collections/cdg/data/CLIVAR_LE/cesm_lens/fx/areacello/areacello_fx_CESM1-CAM5_historical_r0i0p0.nc')
area_truc = CESM_area_file['areacello'].where(CESM_area_file['lat']>35,drop=True) #drop the latitudes below 35N
#weight the grid cells
pi_weighted      = pi_control['aice'] * area_truc.values
pi_slab_weighted = pi_control_slab['aice'] * area_truc.values
#sum up the SIA
pi_SIA      = pi_weighted.sum('nj').sum('ni')
pi_slab_SIA = pi_slab_weighted.sum('nj').sum('ni')
save coupled SIA to NetCDF
pi_SIA.load()
pi_SIA = pi_SIA/1e14

pi_SIA.attrs = {'Description': 'SIA for CESM 1 pre-industrial control (1850) coupled runs, northern hemisphere, time 0400-2200.', 
                'Units'      : 'Million square km',
                'Timestamp'  : str(datetime.datetime.utcnow().strftime("%H:%M UTC %a %Y-%m-%d")),
                'Data source': '/glade/collections/cdg/data/cesmLE/CESM-CAM5-BGC-LE/ice/proc/tseries/monthly/aice/b.e11.B1850C5CN.f09_g16.005.cice.h.aice_nh.*.nc (doi:10.1175/BAMS-D-13-00255.1)',
                'Analysis'   : 'Python 3.7.9 - https://github.com/chrisrwp/synthetic-ensemble/SIA/CESM_pi_control.ipynb'}

pi_SIA.to_netcdf(data_path+'Raw_data/CESM_pi_control/CESM_pi_control_1850_coupled_SIA_0400_2200.nc')

#save slab SIA to NetCDF
pi_slab_SIA.load()
pi_slab_SIA = pi_slab_SIA/1e14

pi_slab_SIA.attrs = {'Description': 'SIA for CESM 1 pre-industrial control (1850) slab-ocean runs, northern hemisphere, time 0100-1001.', 
                     'Units'      : 'Million square km',
                     'Timestamp'  : str(datetime.datetime.utcnow().strftime("%H:%M UTC %a %Y-%m-%d")),
                     'Data source': '/glade/collections/cdg/data/cesmLE/CESM-CAM5-BGC-LE/ice/proc/tseries/monthly/aice/e.e11.E1850C5CN.f09_g16.001.cice.h.aice_nh*.nc (doi:10.1175/BAMS-D-13-00255.1)',
                     'Analysis'   : 'Python 3.7.9 - https://github.com/chrisrwp/synthetic-ensemble/SIA/CESM_pi_control.ipynb'}

pi_slab_SIA.to_netcdf(data_path+'Raw_data/CESM_pi_control/CESM_pi_control_1850_slab_SIA_0100_1001.nc')
#CORRECT THE MONTH??????<jupyter_output><empty_output><jupyter_text>Resample the Preindustrial control data We make 40 ensembles like the CESM1 LE The total length of time of the control runs are 900 or 1800 years We want to test resampling of 40 time series of length 6-250 years<jupyter_code>#load the SIA from computations
pi_SIA      = xr.open_dataarray(data_path+'Raw_data/CESM_pi_control/CESM_pi_control_1850_coupled_SIA_0400_2200.nc')
pi_slab_SIA = xr.open_dataarray(data_path+'Raw_data/CESM_pi_control/CESM_pi_control_1850_slab_SIA_0100_1001.nc')

month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 
               'August', 'September', 'October', 'November', 'December']
def resample_boot2(time_period, data):
    boot_2_first_ind = xr.DataArray(data   = np.random.randint(0,time_period-2, (1000,int(time_period/2))), 
                                    coords = {'resampling':np.arange(1,1001,1), 'year_i':np.arange(1,time_period+1,2)},
                                    dims   = ['resampling', 'year_i'])

    boot_2_second_ind = (boot_2_first_ind+1).copy()
    boot_2_second_ind['year_i'] = np.arange(2,time_period+2,2)

    all_boot_2_ind = xr.concat((boot_2_first_ind, boot_2_second_ind), dim='year_i').sortby('year_i')
    
    ind_base = np.repeat(np.arange(0,time_period*1000,time_period),time_period)
    
    ind_1_d = np.ravel(all_boot_2_ind) + ind_base

    SIA_mem_1000 = np.ravel(np.tile(data,(time_period,1000)))

    resample_boot_2 = xr.DataArray(data = np.reshape(SIA_mem_1000[ind_1_d], (time_period,1000)),
                                   coords = {'year_i':np.arange(1,time_period+1,1), 'resampling':np.arange(1,1001,1)},
                                   dims   = ['year_i', 'resampling'])

    return(resample_boot_2)
def resample_boot1(time_period, data):
    boot_1_ind = xr.DataArray(data   = np.random.randint(0,time_period-1, (1000,int(time_period))), 
                              coords = {'resampling':np.arange(1,1001,1), 'year_i':np.arange(1,time_period+1,1)},
                              dims   = ['resampling', 'year_i'])
    
    ind_base = np.repeat(np.arange(0,time_period*1000,time_period),time_period)
    ind_1_d = np.ravel(boot_1_ind) + ind_base
    
    SIA_mem_1000 = np.ravel(np.tile(data,(time_period,1000)))

    resample_boot_1 = xr.DataArray(data = np.reshape(SIA_mem_1000[ind_1_d], (time_period,1000)),
                                   coords = {'year_i':np.arange(1,time_period+1,1), 'resampling':np.arange(1,1001,1)},
                                   dims   = ['year_i', 'resampling'])

    return(resample_boot_1)
#this function computes resamplings for ONE TIME PERIOD and ONE MONTH
#THIS funtion runs a lot faster as it uses indexing rather than nested for loops
def resample_1period_1month_matrix(SIA, month_, resample_n, time_period, max_t):      
    
    #initialize arrays, shape [time period, member, resampling or largest time series]
    sd_resamp1        = np.empty((40, int(resample_n)))*np.nan
    sd_resamp_no_det1 = np.empty((40, int(resample_n)))*np.nan
    sd_resamp2        = np.empty((40, int(resample_n)))*np.nan
    sd_resamp_no_det2 = np.empty((40, int(resample_n)))*np.nan
    SIA_mem           = np.empty((40, max_t))*np.nan
    mem_anom          = np.empty((40, max_t))*np.nan
    mem_coefs         = np.empty((40, 2))*np.nan
    
    #define the starting years by random
    if str(SIA['time'][0].values)[1] == '4': #then it's the coupled dataset
        start_yrs = np.random.randint(400, 2200-time_period, 40) #generate random start dates
    else:
        start_yrs = np.random.randint(101, 1001-time_period, 40) #generate random start dates

    #################### loop through all the 'members' ###################
    for mem_i in range(40): 
        #define the years to select from the random start date
        yr_list = np.arange(start_yrs[mem_i], start_yrs[mem_i]+time_period, 1)

        #select the years from the SIA time series
        SIA_mem[mem_i][:time_period] = SIA.sel(time=SIA['time.month']==month_).sel(time=slice(str(yr_list[0]).zfill(4),str(yr_list[-1]).zfill(4))).values
        
        #calculate the linear trend coefficients and compute the anomalies from the lienar trend
        mem_coefs[mem_i] = np.polyfit(yr_list, SIA_mem[mem_i][:time_period], 1)
        mem_anom[mem_i][:time_period] = SIA_mem[mem_i][:time_period] - (mem_coefs[mem_i][0] * yr_list + mem_coefs[mem_i][1])

        ###################### compute the resamplings #########################            
        sd_resamp1[mem_i]        = resample_boot1(time_period, mem_anom[mem_i][:time_period]).std('year_i')
        sd_resamp_no_det1[mem_i] = resample_boot1(time_period, SIA_mem[mem_i][:time_period]).std('year_i')
        sd_resamp2[mem_i]        = resample_boot2(time_period, mem_anom[mem_i][:time_period]).std('year_i')
        sd_resamp_no_det2[mem_i] = resample_boot2(time_period, SIA_mem[mem_i][:time_period]).std('year_i')

    output = xr.Dataset(data_vars={'sd_resamp1'        : (('member','resampling'), sd_resamp1),
                                   'sd_resamp_no_det1' : (('member','resampling'), sd_resamp_no_det1),
                                   'sd_resamp2'        : (('member','resampling'), sd_resamp2),
                                   'sd_resamp_no_det2' : (('member','resampling'), sd_resamp_no_det2),
                                   'SIA_mem'           : (('member','year_i'), SIA_mem),
                                   'mem_anom'          : (('member','year_i'), mem_anom),
                                   'mem_coefs'         : (('member', 'coef'), mem_coefs),
                                   'start_yrs'        : (('member'), start_yrs)},
                           
                           coords={'member'     : np.arange(1,41,1), 
                                   'resampling' : np.arange(1,resample_n+1,1), 
                                   'year_i'     : np.arange(1,max_t+1,1),
                                   'coef'       : ['gradient','intercept']})
    
    return(output)<jupyter_output><empty_output><jupyter_text>Use Dask to compute the resampling for all months<jupyter_code>#now try with dask for all months #N.B need to change the function, description, data source, and save path with slab verses coupled
bootsize = 2
resamp_n = 1000
start_t = 6
end_t = 300
inc_t = 4

month_ = 10 

for run_ in np.arange(1,11,1):
    print(datetime.datetime.now(), run_)

    result_list = []

    for t_period in np.arange(start_t, end_t, inc_t):                  #pi_SIA pi_slab_SIA
        result_list.append(dask.delayed(resample_1period_1month_matrix)(pi_slab_SIA, month_, resamp_n, t_period, end_t))
        
    results_computed = dask.compute(*result_list)
    
    results_computed = xr.concat((results_computed), dim='time_period')
    results_computed['time_period'] = np.arange(start_t, end_t, inc_t) 
    
    results_computed['sd_resamp1'].attrs        = {'Description': 'Standard deviation with respect to time for each detrended resampling, bootstrap size of 1 year'}
    results_computed['sd_resamp_no_det1'].attrs = {'Description': 'Standard deviation with respect to time for each resampling without detrending, bootstrap size of 1 year'}
    results_computed['sd_resamp2'].attrs        = {'Description': 'Standard deviation with respect to time for each detrended resampling, bootstrap size of 2 years'}
    results_computed['sd_resamp_no_det2'].attrs = {'Description': 'Standard deviation with respect to time for each resampling without detrending, bootstrap size of 2 years'}
    results_computed['SIA_mem'].attrs           = {'Description': 'The original ramdonly chosen SIA time series for each member without detrending or resampling'}
    results_computed['mem_anom'].attrs          = {'Description': 'Anomalies of the SIA time series relative to a linear trend'}
    results_computed['mem_coefs'].attrs         = {'Description': 'Gradient and y-intercept of the linear trend for each SIA time series'}
    results_computed['start_yrs'].attrs         = {'Description': 'Randomly chosen start years for the SIA time series'}

                                                                                   #coupled slab-ocean      
    results_computed.attrs = {'Description': 'CESM 1 pre-industrial control (1850) slab-ocean run, northern hemisphere sea ice area, month of {}. 40 - so called - members are chosen from random start dates. The length of the times series used varies 6,10,14..298 years. 1000 resamplings are taken for each member'.format(month_names[month_-2]), 
                              'Units'      : 'Million square km',
                              'Timestamp'  : str(datetime.datetime.utcnow().strftime("%H:%M UTC %a %Y-%m-%d")),                #b.e11.B1850C5CN.f09_g16.005 e.e11.E1850C5CN.f09_g16.001
                              'Data source': '/glade/collections/cdg/data/cesmLE/CESM-CAM5-BGC-LE/ice/proc/tseries/monthly/aice/e.e11.E1850C5CN.f09_g16.001.cice.h.aice_nh.*.nc (doi:10.1175/BAMS-D-13-00255.1)',
                              'Analysis'   : 'Python 3.7.9 - https://github.com/chrisrwp/obs-ensemble/CESM_pi_control.ipynb'}
    
    if month_ == 1:
        true_month = '12'
    else:
        true_month = str(month_-1)
                                                                                         #coupled slab                  
    results_computed.to_netcdf('/glade/work/cwpowell/CLIVAR_area_extent_files/pi_control_slab_CESM_SD_time_month_{}_resampling_1000_run_{}.nc'.format(true_month.zfill(2), run_))<jupyter_output>2021-06-10 17:47:12.226544 1
2021-06-10 17:51:38.702983 2
2021-06-10 17:56:06.151409 3
2021-06-10 18:00:23.517518 4
2021-06-10 18:05:04.253485 5
2021-06-10 18:09:15.912360 6
2021-06-10 18:13:57.209075 7
2021-06-10 18:17:59.071629 8
2021-06-10 18:22:28.026551 9
2021-06-10 18:26:29.918819 10
{'fraction_non_alphanumeric': 0.12558229813664595, 'fraction_numerical': 0.057971014492753624, 'mean_word_length': 3.211716621253406, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 22, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '305425', 'license': 'MIT', 'n_tokens_mistral': 6027, 'n_tokens_neox': 5326, 'n_words': 1114, 'path': 'SIA/CESM_pi_control.ipynb', 'repo_name': 'chrisrwp/synthetic-ensemble', 'size': '15469'}