<jupyter_start><jupyter_text>Plot all regime data<jupyter_code># First we check that they look OK, and have the same orientation.
plt.imshow(np.flip(labels2[0][:][:]))
plt.figure()<jupyter_output><empty_output><jupyter_text>Original pseudo codedef creaCalc(data, areaData):    Get data and areaData into numpy arrays    Calculate area for the data for all the regimes    (This will be an array where only the areas that are covered by the different dynamical regimes are non-nan)    return areaDataRegime0, areaDataRegime1, areaDataRegime2, areaDataRegime3, areaDataRegime4, areaDataRegime5<jupyter_code>def extractionAreaCalc(data, area,chunks=None):
    """Code to extract area of different regimes. In principle this is applicable to any field, not just the area."""
    '''
    input:
    data is the labelled regime data 
    area is the total area
    output:
    mask matrix, indexed by regime labels 0 to 5
    '''
    regimeData=np.array(data)         #May have to add selection of specific data
    areaData=np.array(area)           #May have to add selection of specific data
   # print(areaData.shape)
    #Create dummy variable
    Msk=np.zeros((6,chunks, data.shape[1], data.shape[2]))
    for chunk in np.arange(0,chunks): #this will run for 7 chunks of 20 years
        for regime in np.arange(0,6):     # Loop over regimes
            regimeIndexes=np.where(regimeData[chunk]!=regime)   # Determine where regimes are *NOT* present
            areaSelected=np.copy(areaData[0])
            areaSelected[regimeIndexes]=-9999 # Set areas there the regime is not present to zero. Can be np.nan too, but I'm not sure what will caouse the least trouble!
            Msk[regime][chunk]=areaSelected      # Slot into results matrix
           # plt.figure()                  # Can be commented out
           # plt.title("Chunk-"+(str)(chunk)+"-- Area of regime "+(str)(regime))
            #plt.imshow(np.rot90(Msk[regime][chunk])) #*land,1))       # Can be commented out
            #print(areaSelected[(areaSelected>0)])
            #print("test",Msk[regime][chunk][Msk[regime][chunk] > 0])
    return Msk                        # Decided to return one matric instead of 6 fields

#todo concatenate chunks in regimeData
labels2.shape
areaData=np.array(areaGlobal)
Mask=extractionAreaCalc(labels2, areaGlobal,chunks=7)
Mask.shape
labels2.shape
#new_labels2=np.zeros((1, labels2[0][:][:].shape[0], labels2[0][:][:].shape[1]))

#new_labels2[0] = labels2[0][:]
plt.imshow(labels2[0])
da_regimeData = xr.DataArray(
        data=labels2,
        dims=["chunk","y", "x"],
)

da_regimeData.shape
ds = xr.Dataset(
     data_vars=dict(
        regimeData=(["chunk","y","x"], da_regimeData),
        areaGlobal=(ds2.areacello),
        mask=(["regimeIndex","chunk","y", "x"], Mask),
    ),
     coords=dict(
         x=(ds2.x),
         y=(ds2.y),
         chunk=[0, 1, 2 , 3 , 4, 5, 6]
     )
)
ds
ds.regimeData[0].where(ds.regimeData[0] == 1).plot(cmap=my_cmap)
ds.mask[1][0].where(ds.mask[1][0] != -9999).plot(cmap=my_cmap)#        mask=(["regimeIndex","chunk","y", "x"], Mask),

ds.regimeData[0].where(ds.regimeData[0] == 1).plot(cmap=my_cmap)
regime_nino3 = ds.regimeData.sel(y = slice(-5,5), x = slice(-150, -90))
regime_nino3[0].plot(cmap=my_cmap)
mask_nino3 = ds.mask.sel(y = slice(-5,5), x = slice(-150, -90))

mask_nino3.shape
area_sum = np.zeros((3, 6))
#ds.mask[1][0].where(ds.mask[1][0] != -1).sel(y = slice(-5,5), x = slice(-150, -90)).plot()
c = 0
for chunk in np.arange(4,7):
    for regime in np.arange(0,6): 
        #plt.figure()
        mask_nino3_regime = ds.mask[regime][chunk].where(ds.mask[regime][chunk] != -9999).sel(y = slice(-5,5), x = slice(-150, -90))
        #mask_nino3_regime.plot()
        avgarea = mask_nino3_regime.sum(dim=('x', 'y'))
        area_sum[c][regime] = avgarea
        #plt.title("Regime:"+(str)(regime)+"Average area:"+(str)(avgarea))
    c+=1

print ("Mean and std of 7 chunks of 20 years")
mean_his=np.zeros(6)
std_his=np.zeros(6)
for regime in range(0,6):
    all_chuncks = np.zeros(3)
    for yrs in range(0,3):
        all_chuncks[yrs] = area_sum[yrs][regime]
    mean_his[regime], std_his[regime] = all_chuncks.mean(), all_chuncks.std()
    print("mean= ", "{:e}".format(all_chuncks.mean()), "std= ", "{:e}".format(all_chuncks.std()))<jupyter_output>Mean and std of 7 chunks of 20 years
mean=  6.916359e+12 std=  1.128402e+10
mean=  1.757886e+11 std=  3.959768e+09
mean=  3.126645e+11 std=  1.472449e+10
mean=  0.000000e+00 std=  0.000000e+00
mean=  0.000000e+00 std=  0.000000e+00
mean=  7.941738e+09 std=  0.000000e+00<jupyter_text>1%CO2<jupyter_code>labels_1pct = np.load("/home/jovyan/DNN4Cli/figures/ESM4/ESM4-1pctCO2-000101-014012-EnsembleMLP.npy")
labels_1pct.shape
plt.imshow(labels_1pct[0])
Mask_1pct=extractionAreaCalc(labels_1pct, areaGlobal,chunks=7)
#new_labels_1pct=np.zeros((1, labels_1pct.shape[1], labels_1pct.shape[2]))

#new_labels_1pct = labels_1pct[:]
plt.imshow(labels_1pct[2])
da_regimeData_1pct = xr.DataArray(
        data=labels_1pct,
        dims=["chunk","y", "x"],
)
ds_1pct = xr.Dataset(
     data_vars=dict(
        regimeData=(["chunk","y","x"], da_regimeData_1pct),
        areaGlobal=(ds2.areacello),
        mask=(["regimeIndex","chunk","y", "x"], Mask_1pct),
    ),
     coords=dict(
         x=(ds2.x),
         y=(ds2.y),
         chunk=range(7)
     )
)
ds_1pct
ds_1pct.regimeData[0].where(ds_1pct.regimeData[0] == 1).plot(cmap = my_cmap)
regime_nino3_1pct = ds_1pct.regimeData.sel(y = slice(-5,5), x = slice(-150, -90))
regime_nino3_1pct[5].plot(cmap=my_cmap)
area_sum_1pct = np.zeros((3,6))
c = 0
for chunk in np.arange(4,7):
    for regime in np.arange(0,6): 
        #plt.figure()
        mask_nino3_regime_1pct = ds_1pct.mask[regime][chunk].where(ds_1pct.mask[regime][chunk] != -9999).sel(y = slice(-5,5), x = slice(-150, -90))
        #mask_nino3_regime.plot()
        avgarea = mask_nino3_regime_1pct.sum(dim=('x', 'y'))
        area_sum_1pct[c][regime] = avgarea
        #plt.title("Regime:"+(str)(regime)+"Average area:"+(str)(avgarea))
    c+=1
print ("Mean and std of 3 chunks of 20 years")
mean_1pct = np.zeros(6)
std_1pct = np.zeros(6)
for regime in range(0,6):
    all_chuncks = np.zeros(7)
    for yrs in range(0,3):
        all_chuncks[yrs] = area_sum_1pct[yrs][regime]
    mean_1pct[regime], std_1pct[regime] = all_chuncks.mean(), all_chuncks.std()
    print("mean= ", "{:e}".format(all_chuncks.mean()), "std= ", "{:e}".format(all_chuncks.std()))<jupyter_output>Mean and std of 3 chunks of 20 years
mean=  2.845594e+12 std=  3.286027e+12
mean=  1.620106e+11 std=  1.931602e+11
mean=  1.658861e+11 std=  1.920101e+11
mean=  0.000000e+00 std=  0.000000e+00
mean=  0.000000e+00 std=  0.000000e+00
mean=  3.403602e+09 std=  3.930141e+09<jupyter_text>Abrupt-4xCO2<jupyter_code>labels_4x = np.load("/home/jovyan/DNN4Cli/figures/ESM4/ESM4-abrupt4xCO2-000101-014012-EnsembleMLP.npy")
labels_4x.shape
plt.imshow(labels_4x[0])
Mask_4x = extractionAreaCalc(labels_4x, areaGlobal,chunks=7)
Mask_4x.shape
#new_labels_4x=np.zeros((1, labels_4x[0][:][:].shape[0], labels_4x[0][:][:].shape[1]))

#new_labels_4x[0] = labels_4x[0][:]
da_regimeData_4x = xr.DataArray(
        data=labels_4x,
        dims=["chunk","y", "x"],
)
ds_4x = xr.Dataset(
     data_vars=dict(
        regimeData=(["chunk","y","x"], da_regimeData_4x),
        areaGlobal=(ds2.areacello),
        mask=(["regimeIndex","chunk","y", "x"], Mask_4x),
    ),
     coords=dict(
         x=(ds2.x),
         y=(ds2.y),
         chunk=np.arange(7)
     )
)
ds_4x.regimeData[0].where(ds_4x.regimeData[0] == 1).plot(cmap=my_cmap)
ds_4x.mask[1][0].where(ds_4x.mask[1][0] != -9999).plot(cmap=my_cmap)#        mask=(["regimeIndex","chunk","y", "x"], Mask),

regime_nino3_4x = ds_4x.regimeData.sel(y = slice(-5,5), x = slice(-150, -90))
regime_nino3_4x[6].plot(cmap=my_cmap)
area_sum_4x =np.zeros((3,6))
c=0
for chunk in np.arange(4,7):
    for regime in np.arange(0,6): 
        #plt.figure()
        mask_nino3_regime_4x = ds_4x.mask[regime][chunk].where(ds_4x.mask[regime][chunk] != -9999).sel(y = slice(-5,5), x = slice(-150, -90))
        #mask_nino3_regime.plot()
        avgarea = mask_nino3_regime_4x.sum(dim=('x', 'y'))
        area_sum_4x[c][regime] = avgarea
        #plt.title("Regime:"+(str)(regime)+"Average area:"+(str)(avgarea))
    c+=1

area_sum_4x
print ("Mean and std of 7 chunks of 20 years")
mean_4x = np.zeros(6)
std_4x = np.zeros(6)
for regime in range(0,6):
    all_chuncks = np.zeros(7)
    for yrs in range(0,3):
        all_chuncks[yrs] = area_sum_4x[yrs][regime]
    mean_4x[regime], std_4x[regime] = all_chuncks.mean(), all_chuncks.std()
    print("mean= ", "{:e}".format(all_chuncks.mean()), "std= ", "{:e}".format(all_chuncks.std()))
import pandas as pd
df = pd.DataFrame({"Mean_his": mean_his,
                  "Std_his": std_his,
                  "Mean_1%": mean_1pct,
                  "Std_1%": std_1pct,
                  "Mean_4x": mean_4x,
                  "Std_4x": std_4x
                 })
df
area_sum_1pct
area_sum_4x
MD = area_sum_1pct[:,0]

S_NV =area_sum_1pct[:,1]

N_NV = area_sum_1pct[:,2]

TR = area_sum_1pct[:,5]
areass = [MD, S_NV, N_NV, TR]
MD4x = area_sum_4x[:,0]
S_NV4x =area_sum_4x[:,1]
N_NV4x = area_sum_4x[:,2]
TR4x = area_sum_4x[:,5]
areass4x = [MD4x, S_NV4x, N_NV4x, TR4x]
areass
areass4x
labels = ['Area MD', 'Area SSV', 'Area NSV', 'Area NL'
         ]
area_sum = area_sum[2].tolist()

area_sum
#area_sum.pop(3)
area_sum
labels
plt.plot(1, area_sum[0],'*', markersize=15)
plt.plot(2, area_sum[1],'*', markersize=15)
plt.plot(3, area_sum[2],'*', markersize=15)
plt.plot(4, area_sum[3],'*', markersize=15)
plt.figure(figsize=(8,6))
bplot1 = plt.boxplot(areass,
                      vert=True,  # vertical box alignment
                  patch_artist=True,# fill with color
                     widths = 0.5,
                labels=labels) 
bplot2 = plt.boxplot(areass4x,
                      vert =True,  # vertical box alignment
                  patch_artist=True,# fill with color
                     widths = 0.5,
                labels=labels)
#plothis = plt.plot(labels, area_sum, '*', "green")
plt.plot(1, area_sum[0],'*', markersize=15,color = "green")
plt.plot(2, area_sum[1],'*', markersize=15 ,color = "green")
plt.plot(3, area_sum[2],'*', markersize=15,color = "green")
hist = plt.plot(4, area_sum[3],'*', markersize=15,color = "green")

plt.yscale('log')
#plt.ylim (9.0e+11, 6.0e+12)
#(['A simple line'])
colors = ['darkblue','green', 'pink', 'lightblue']
colors = ["blue", "blue", "blue", "blue"]
#for bplot in (bplot1, bplot2):
    #for patch, color in zip(bplot['boxes'], colors):
        #patch.set_facecolor(color)
for patch in bplot1["boxes"]:
    patch.set(facecolor = 'blue')
for patch in bplot2["boxes"]:
    patch.set(facecolor = 'red')
plt.title("Different Regimes Area StD ESM4")
plt.legend([hist, bplot1["boxes"][0], bplot2["boxes"][0]], ["Historical",'1%CO'+'\N{SUBSCRIPT TWO}', '4xCO'+'\N{SUBSCRIPT TWO}'], loc='upper right')
plt.show()

#plt.ylim (10.0e+11,)<jupyter_output><ipython-input-90-0054a318c980>:31: UserWarning: Legend does not support [<matplotlib.lines.Line2D object at 0x7ff26c9dd3d0>] instances.
A proxy artist may be used instead.
See: https://matplotlib.org/users/legend_guide.html#creating-artists-specifically-for-adding-to-the-legend-aka-proxy-artists
  plt.legend([hist, bplot1["boxes"][0], bplot2["boxes"][0]], ["Historical",'1%CO'+'\N{SUBSCRIPT TWO}', '4xCO'+'\N{SUBSCRIPT TWO}'], loc='upper right')
{'fraction_non_alphanumeric': 0.17637314734088927, 'fraction_numerical': 0.06242371403661726, 'mean_word_length': 3.815701091519731, 'pattern_counts': {'":': 6, '<': 16, '<?xml version=': 0, '>': 18, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '204301', 'license': 'MIT', 'n_tokens_mistral': 5014, 'n_tokens_neox': 4688, 'n_words': 838, 'path': 'NINOArea_Calc/esm4_areaForRegimes.ipynb', 'repo_name': 'aradhakrishnanGFDL/DNN4Cli', 'size': '11470'}