<jupyter_start><jupyter_text>Draw regions found by the metrics in the mapThis notebook is used to obtain the map drawings.<jupyter_code>import matplotlib.pyplot as plt
import pandas as pd
import matplotlib.cm as cm
import numpy as np
import geopandas as gpd
from scipy.sparse import coo_matrix,csr_matrix,save_npz,load_npz
%matplotlib inline
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('retina')
import os
import math
import numpy as np
import matplotlib
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['ps.fonttype'] = 42
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.patches import Rectangle
import matplotlib.colors as mpc
import matplotlib.ticker as ticker
from collections import Counter
import matplotlib.patches as mpatches
import matplotlib as mpl
from tqdm.auto import tqdm
os.makedirs("../Figures/", exist_ok=True)
%load_ext Cython
%%cython

cimport numpy as np
import numpy as np
cimport cython
ctypedef np.float64_t DTYPE_t
from math import radians, cos, sin, asin, sqrt

def haversine(lng1, lat1, lng2, lat2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians 
    lng1, lat1, lng2, lat2 = map(radians, [lng1, lat1, lng2, lat2])

    # haversine formula 
    dlng = abs(lng2 - lng1)
    if dlng > np.pi: #warp around longitude
        dlng = 2*np.pi-dlng
    dlat = abs(lat2 - lat1)
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlng/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r

def manhatton(lng1, lat1, lng2, lat2):
    """
    Calculate the Manhattan distance between two points 
    on the cylinder (all gird distance are transformed into equator scale)
    """
    # convert decimal degrees to radians 
    lng1, lat1, lng2, lat2 = map(radians, [lng1, lat1, lng2, lat2])

    # haversine formula 
    dlng = abs(lng2 - lng1)
    if dlng > np.pi: #warp around longitude
        dlng = 2*np.pi-dlng
    dlat = abs(lat2 - lat1)
    c = dlng + dlat
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles
    return c * r

def geoKernel_cython(float[:,::1] data, float threshold):
    
    n = len(data)
    X = np.empty((n,n), dtype=float)
    for i in range(n):
        for j in range(i, n):
            # call the Haversine formula
            dist = manhatton(data[i][0],data[i][1],data[j][0],data[j][1])
            if (dist <= threshold)&(dist>0.1):
                X[i][j] = dist
            else:
                X[i][j] = 0
    return X
# Loads or recalculates and store geokernel
# set reCalculateGeoKernel to True to recalculate GeoKernel
reCalculateGeoKernel = False
def getGeoKernel(dataPath,nx,ny):
    if (reCalculateGeoKernel):
        dfCorrelation = pd.read_csv(dataPath, sep=',')
        pointsCorrelation = gpd.GeoDataFrame(dfCorrelation, geometry=gpd.points_from_xy(dfCorrelation.lng, dfCorrelation.lat))

        data = pointsCorrelation[['lng','lat']].to_numpy()
        data = data.copy(order='C').astype('float32') 
        data.flags
        print(data);

        N, thresh = len(data), 200
        X = geoKernel_cython(data, thresh)

        X[X<0.0000001] = 0
        Xs = X + X.T - np.diag(X.diagonal())
        fullGeo = coo_matrix(Xs)
        save_npz("../Data/geokernel_nx%d_ny%d.npz"%(nx,ny), fullGeo)
    else:
        fullGeo = load_npz("../Data/geokernel_nx%d_ny%d.npz"%(nx,ny))
    return fullGeo;
# efficient sparse geo-weight-dist hadamard-product with input data. Notice this directly cacluates a distance matrix instead of a kernel
def geoWeightedDist(data, geoDist, smoothness):
    # creating masks for input data given the geo-kernel
    value_i = data[geoDist.row]
    value_j = data[geoDist.col]
    # only calculate between points that is non-zero in geo-kernel
    #Xd = np.subtract(value_i,value_j)
    Xshrink = np.multiply(value_i,value_j).reshape(len(value_i),)#Shrinking factor based on value, smaller the value, closer the distance
    #Xd = np.exp(Xd**2/smoothness**2).reshape(len(value_i),)
    #geoSmooth = geoDist.copy()
#     print(Xshrink)
    print(geoDist.data**2)
    print(np.power(Xshrink,smoothness))
    print(geoDist.data.std()**2)
    Xd = np.exp(geoDist.data**2*np.power(Xshrink,smoothness)/geoDist.data.std()**2)
    #Xd = np.multiply(Xd,geoSmooth.data)
#     print(Xd)
    C = coo_matrix((Xd/max(Xd), (geoDist.row, geoDist.col)), shape=(len(data),len(data)))
    return C
def drawNewMap(filename = 'Data/Correlation_d7_l180_v3.csv',nx=280,ny=140,percentile = 0.30,pvalueThreshold = None,smoothness = 0.5, quantile = 0.5, extraEpsilonFactor = 0.5, min_samples =5, min_points = 20,exportName="",useCircular=False,maxPositivePValue = None):
    dfCorrelation = pd.read_csv(filename, sep=',')
  
    gridRegions = [];
    for coords in dfCorrelation[["lat","lng"]].to_numpy():
        region = (coords[0] - 70/ny, coords[0]+70/ny,coords[1]-180/nx,coords[1]+180/nx)
        gridRegions.append(region)
  
  
    dfCorrelation["region"] = gridRegions;
  
  
    # pointsCorrelation = gpd.GeoDataFrame(dfCorrelation, geometry=gpd.points_from_xy(dfCorrelation.lng, dfCorrelation.lat))
    pointsCorrelation = dfCorrelation.reset_index(drop=True)
  
    if(percentile<1.0):
        thresholdCorrelation = np.quantile(pointsCorrelation[['pvalue']],percentile);
    else:
        thresholdCorrelation = 1.0
      
    if(pvalueThreshold is not None):
        thresholdCorrelation = pvalueThreshold;
      
    pointsCorrelation['log-pvalue'] = np.log10(pointsCorrelation['pvalue'])
    pointsCorrelation['thresh-pvalue'] = pointsCorrelation['pvalue']
  
  
    pointsCorrelation.loc[pointsCorrelation['thresh-pvalue']>thresholdCorrelation,'thresh-pvalue']=1
  
    points = pointsCorrelation;
  
    print("Chosen p-value: %.2g"%thresholdCorrelation);
    print("Chosen correlation: %.2g"%np.max(pointsCorrelation.loc[pointsCorrelation['thresh-pvalue']>thresholdCorrelation,'value']));
  
    tab20b = cm.get_cmap('tab20b', 20)
    tab20c = cm.get_cmap('tab20c', 20)
    newcmp = mpc.ListedColormap(tab20b(np.linspace(0.25, 0.75, 256)))
    a = np.concatenate([tab20b(range(20)),tab20c(range(20))])
    np.random.shuffle(a)
    newcmp = mpc.ListedColormap(a)
    from sklearn.cluster import DBSCAN
    def identity_scale(minval, maxval):
        def scalar(val):
            return 2
        return scalar
  
    # pthreshold=0.05
    fullGeo = getGeoKernel(filename,nx,ny)
    if(quantile<1.0):
        plotGraph = geoWeightedDist(points[['thresh-pvalue']].to_numpy(),fullGeo,smoothness)
        
        clustering = DBSCAN(eps=np.quantile(plotGraph.data, quantile)*extraEpsilonFactor, min_samples=min_samples, metric='precomputed').fit(plotGraph)
        labels = clustering.labels_
  
        smallerGroups = set();
        allLabels = list(map(str,(labels)));
        for label,count in Counter(allLabels).items():
            if(count<min_points):
                smallerGroups.add(label)
  
        #labels2 = spectral_clustering(graph2, n_clusters=10, eigen_solver='arpack')
        # plt.hist(labels, normed=True, bins=30)
        toPlotPoints = points.copy()
        toPlotPoints['labels'] = [entry if entry not in smallerGroups else "-1" for entry in allLabels]
        toPlotData = toPlotPoints[(toPlotPoints["labels"]!="-1") * (toPlotPoints["thresh-pvalue"]<thresholdCorrelation)]
    else:
        toPlotData = points.copy();
  
  
    ENSORegion = [[-5 , -170] ,[5, -120]];
  
    cdict1 = {'red':   ((0.0, 0.0, 0.0),
                       (0.5, 0.5, 1.0),
                       (1.0, 1.0, 0.0)),
  
             'green': ((0.0, 0.0, 0.0),
                       (0.5, 1.0, 1.0),
                       (1.0, 0.0, 0.0)),
  
             'blue':  ((0.0, 0.0, 1.0),
                       (0.5, 1.0, 0.0),
                       (1.0, 0.0, 0.0))
            }
  
    ccmap = matplotlib.cm.get_cmap("autumn");#LinearSegmentedColormap("bluered", cdict1)
  
    if(maxPositivePValue is None):
        maxPositivePValue = np.min(toPlotData["log-pvalue"][np.isfinite(toPlotData["log-pvalue"])])
    print(maxPositivePValue)
  #   maxPositivePValue = -40
  
    if(True):
        if(useCircular):
            projection = ccrs.Mollweide(central_longitude=150)
        else:
            projection = ccrs.PlateCarree(central_longitude=150);
  
        fig = plt.figure(figsize=(10,8),tight_layout=True)
        ax = plt.axes(projection=projection)
        # ax.set_extent([-30, 150, -70, 70])
        # ax.stock_img()
        ax.set_global()
#         if(not useCircular):
#             ax.set_extent([-30, 150, -70, 70])
        ax.add_feature(cfeature.LAND,facecolor='gray',alpha=1.0)
        ax.add_feature(cfeature.OCEAN,facecolor=(0.7,0.7,0.7),alpha=1.0)
        # ax.add_feature(cfeature.OCEAN)
  
      #Arrows lines 
  
        transform = projection._as_mpl_transform(ax)
      
        def plotRectangle(ax,region,name,color,transform = ccrs.Geodetic(),lw=0,ec=None):
            if(transform):
                ax.add_patch(mpatches.Rectangle(xy=[region[2], region[0]], width=(region[3]-region[2]), height=(region[1]-region[0]),
                                              fc=color,
                                              ec=None,
                                              lw=0,
                                              transform=transform))
            else:
                ax.add_patch(mpatches.Rectangle(xy=[region[2], region[0]], width=(region[3]-region[2]), height=(region[1]-region[0]),
                                              fc=color,
                                              lw=0,
                                              ec=None))
            if(name):
                ax.text((region[3]+region[2])*0.5, (region[0]+region[1])*0.5, name,horizontalalignment='center',verticalalignment='center', transform=ccrs.Geodetic(),zorder=6)
  
      
        norm = mpl.colors.Normalize(vmin=maxPositivePValue, vmax=np.log10(thresholdCorrelation))
      
        for region,pvalue in  toPlotData[["region","thresh-pvalue"]].to_numpy():
            regionPvalue = pvalue;
            regionPvalueLog = np.log10(regionPvalue);
            lat,lng = np.mean(region[0:1]),np.mean(region[2:3]);
            if(np.isfinite(regionPvalueLog) and regionPvalueLog>=maxPositivePValue):
                color=bluered(norm(regionPvalueLog));
            else:
                color=(0.5,0.0,0.0,1.0);
  
            plotRectangle(ax,region,"",color,transform=ccrs.Geodetic());
  
        plotRectangle(ax,[-5 , 5 ,-170, -120],"ENSO","#CCCCCC",transform=ccrs.Geodetic(),lw=0.2,ec=(0.1,0.1,0.1));
  
  
        ax.add_feature(cfeature.COASTLINE,edgecolor='black',alpha=0.5)
        ax.add_feature(cfeature.BORDERS,edgecolor='black',alpha=0.5)
  
  
        import json
  
        with open("../Data/ElNidoShapes.json","r") as fd:
            shapes = [shape[0] for shape in json.load(fd)];
  
        with open("../Data/LaNinaShapes.json","r") as fd:
            shapes += [shape[0] for shape in json.load(fd)];
  
  
  
        shapesData = [np.array(shape) for shape in shapes]  
  
        for i in range(len(shapesData)):
            polygon = mpatches.Polygon(shapesData[i],
              closed=False,
              facecolor='white',
              alpha=0.2,
              transform=ccrs.Geodetic())
            ax.add_patch(polygon)
  
  
  
        print(maxPositivePValue);
        if(exportName):
            fig.savefig("../Figures/"+exportName+".pdf");
        plt.show()

    fig, ax = plt.subplots(figsize=(6, 1))
    fig.subplots_adjust(bottom=0.5)
    
    def fmt(x, pos):
        return '$10^{%d}$'%(x);
  
    if(-maxPositivePValue<10):
        bounds=np.arange(np.ceil(maxPositivePValue),round(np.log10(thresholdCorrelation))+1,1)
    elif(-maxPositivePValue<50):
        bounds=np.arange(np.ceil(maxPositivePValue),round(np.log10(thresholdCorrelation))+1,5)
    elif(-maxPositivePValue<100):
        bounds=np.arange(np.ceil(maxPositivePValue),round(np.log10(thresholdCorrelation))+1,10)
    elif(-maxPositivePValue<200):
        bounds=np.arange(np.ceil(maxPositivePValue),round(np.log10(thresholdCorrelation))+1,20)
    else:
        bounds=np.arange(np.ceil(maxPositivePValue),round(np.log10(thresholdCorrelation))+1,50)
      
    print(bounds);
    ccmap = matplotlib.cm.get_cmap("autumn")
    ccmap.set_under(color=(0.5,0.0,0.0,1.0))
    norm = mpl.colors.Normalize(vmin=maxPositivePValue, vmax=np.log10(thresholdCorrelation))
    cb1 = mpl.colorbar.ColorbarBase(ax, cmap=ccmap,
                                    norm=norm,
  #                                   boundaries=[-maxPositivePValue] + bounds + [maxPositivePValue],
                                    extend='min',
                                    extendfrac=0.05,
                                    spacing='uniform',
                                    ticks=bounds,
                                    orientation='horizontal',
                                   format=ticker.FuncFormatter(fmt))
  
    if(exportName):
        plt.savefig("../Figures/cm_%s.pdf"%exportName)
    fig.show()

#This creates maps comparing G Causality and Correlation
nx = 288
ny = 140
# nx = 36
# ny = 14
#Use to select the decade or the whole data (0,100)
dataRange = (0,100) #percentage
# 0,26 1980s
# 26,52 1990s
# 52,78 2000s
# 78,100 2010s
windowDays = 7
maxLag = 90
highResolution = (288,140)
lowResolution = (36,14)

dataRangeLags = [
#     ((0,26),90,highResolution),  # 1980s
#     ((26,52),90,highResolution), # 1990s
#     ((52,78),90,highResolution), # 2000s
#     ((78,100),90,highResolution),# 2010s
#     ((0,100),90,highResolution), # All data lag=90 days
    ((0,100),180,highResolution), # All data lag=180 days
#     ((0,100),90,lowResolution),
]

np.seterr(divide='ignore')
for dataRange,maxLag,resolution in tqdm(dataRangeLags):
    nx=resolution[0]
    ny=resolution[1]
    percentile = 0.30
    pvalueThreshold = None
    smoothness = 0.5
    quantile = 0.5
    extraEpsilonFactor = 0.5
    min_samples =5
    min_points = 20
    exportName=""
    useCircular=False
    maxPositivePValue = None
    smoothness = 0.5

    extraEpsilonFactor = 0.5
    min_samples = 5
    min_points = 30

    for pvalueThreshold in [0.001,0.005,0.01,0.05]:
        exportName="region_differences_7-30_l%d_nx%d_ny%d_R%d-%d_p%.5f_v8"%(maxLag,nx,ny,dataRange[0],dataRange[1],pvalueThreshold)
        maxPositivePValue=maxPositivePValue

        dfsAllDataCausality = []
        dfsAllDataCorrelation = []
        for index,windowDays in enumerate(tqdm(range(7,31))):
            filenameCausality = "../Data/Causality_d%d_l%d_nx%d_ny%d_R%d-%d_v8.csv"%(windowDays,maxLag,nx,ny,dataRange[0],dataRange[1])
            filenameCorrelation = "../Data/LaggedCorrelation_d%d_l%d_nx%d_ny%d_R%d-%d_v8.csv"%(windowDays,maxLag,nx,ny,dataRange[0],dataRange[1])
            dfCausality = pd.read_csv(filenameCausality, sep=',')
            dfCorrelation = pd.read_csv(filenameCorrelation, sep=',')
            if(index==0):
                gridRegions = [];
                for coords in dfCausality[["lat","lng"]].to_numpy():
                    region = (coords[0] - 70/ny, coords[0]+70/ny,coords[1]-180/nx,coords[1]+180/nx)
                    gridRegions.append(region)
            dfCausality["region"] = gridRegions;

            dfCorrelation["region"] = gridRegions;

            pointsCausality = dfCausality.reset_index(drop=True)
            pointsCorrelation = dfCorrelation.reset_index(drop=True)

            pointsCausality['log-pvalue'] = np.log10(pointsCausality['pvalue'])
            pointsCausality['thresh-pvalue'] = pointsCausality['pvalue']

            pointsCorrelation['log-pvalue'] = np.log10(pointsCorrelation['pvalue'])
            pointsCorrelation['thresh-pvalue'] = pointsCorrelation['pvalue']


            if(percentile<1.0):
                thresholdCausality = np.quantile(pointsCausality[['pvalue']],percentile);
            else:
                thresholdCausality = 1.0
            if(pvalueThreshold is not None):
                thresholdCausality = pvalueThreshold;


            if(percentile<1.0):
                thresholdCorrelation = np.quantile(pointsCorrelation[['pvalue']],percentile);
            else:
                thresholdCorrelation = 1.0
            if(pvalueThreshold is not None):
                thresholdCorrelation = pvalueThreshold;


            pointsCausality.loc[pointsCausality['thresh-pvalue']>thresholdCausality,'thresh-pvalue']=1
            pointsCorrelation.loc[pointsCorrelation['thresh-pvalue']>thresholdCorrelation,'thresh-pvalue']=1

            dfsAllDataCausality.append(pointsCausality)
            dfsAllDataCorrelation.append(pointsCorrelation)

        countsCausality = np.zeros(len(dfsAllDataCausality[0]))
        numberOfDays = len(dfsAllDataCausality)
        for df in dfsAllDataCausality:
            countsCausality[df["thresh-pvalue"]<=thresholdCausality] +=1
        #     counts += df["log-pvalue"]/numberOfDays

        countsCorrelation = np.zeros(len(dfsAllDataCorrelation[0]))
        numberOfDays = len(dfsAllDataCorrelation)
        for df in dfsAllDataCorrelation:
            countsCorrelation[df["thresh-pvalue"]<=thresholdCorrelation] +=1
        #     counts += df["log-pvalue"]/numberOfDays


        dfsAllDataCausality[0]["counts"] = countsCausality
        dfsAllDataCorrelation[0]["counts"] = countsCorrelation
        countThreshold=5

        saveFilenameCausality = "../Data/%s_dAll_l%d_nx%d_ny%d_R%d-%d_p%.5f_v8.csv"%("Causality",maxLag,nx,ny,dataRange[0],dataRange[1],pvalueThreshold)
        dfsAllDataCorrelation[0].to_csv(saveFilenameCausality,sep=",")

        saveFilenameCorrelation = "../Data/%s_dAll_l%d_nx%d_ny%d_R%d-%d_p%.5f_v8.csv"%("LaggedCorrelation",maxLag,nx,ny,dataRange[0],dataRange[1],pvalueThreshold)
        dfsAllDataCorrelation[0].to_csv(saveFilenameCorrelation,sep=",")

        dfsAllDataCausality[0]["diff"] = countsCausality-countsCorrelation
        dfsAllDataCausality[0]["absolute"] = np.max([countsCausality,countsCorrelation],axis=0)

        import copy
        toPlotData = dfsAllDataCausality[0]
        tab20b = cm.get_cmap('tab20b', 20)
        tab20c = cm.get_cmap('tab20c', 20)
        newcmp = mpc.ListedColormap(tab20b(np.linspace(0.25, 0.75, 256)))
        a = np.concatenate([tab20b(range(20)),tab20c(range(20))])
        np.random.shuffle(a)
        newcmp = mpc.ListedColormap(a)
        from sklearn.cluster import DBSCAN
        def identity_scale(minval, maxval):
            def scalar(val):
                return 2
            return scalar

        useCircular= True
        ENSORegion = [[-5 , -170] ,[5, -120]];

        cdict1 = {'red':   ((0.0, 0.0, 0.0),
                           (0.5, 0.5, 1.0),
                           (1.0, 1.0, 0.0)),

                 'green': ((0.0, 0.0, 0.0),
                           (0.5, 1.0, 1.0),
                           (1.0, 0.0, 0.0)),

                 'blue':  ((0.0, 0.0, 1.0),
                           (0.5, 1.0, 0.0),
                           (1.0, 0.0, 0.0))
                }

        vmin = -len(dfsAllDataCausality)
        vmax = len(dfsAllDataCausality)
        norm = mpl.colors.Normalize(vmin=vmin-0.5, vmax=vmax+0.5)
        normOpacity = mpl.colors.Normalize(vmin=-.05, vmax=vmax+0.5)

        ccmap = copy.copy(matplotlib.cm.get_cmap("Spectral_r",vmax-vmin+1))

        if(True):
            if(useCircular):
                projection = ccrs.Robinson(central_longitude=150)
            else:
                projection = ccrs.PlateCarree(central_longitude=150);

            fig = plt.figure(figsize=(10,8),tight_layout=True)
            ax = plt.axes(projection=projection)
            # ax.set_extent([-30, 150, -70, 70])
            # ax.stock_img()
            ax.set_global()
        #         if(not useCircular):
        #             ax.set_extent([-30, 150, -70, 70])
            ax.add_feature(cfeature.LAND,facecolor='gray',alpha=1.0)
            ax.add_feature(cfeature.OCEAN,facecolor=(0.7,0.7,0.7),alpha=1.0)
            # ax.add_feature(cfeature.OCEAN)

          #Arrows lines 

            transform = projection._as_mpl_transform(ax)

            def plotRectangle(ax,region,name,color,transform = ccrs.Geodetic(),lw=0,ec=None,alpha=1.0):
                if(transform):
                    ax.add_patch(mpatches.Rectangle(xy=[region[2], region[0]], width=(region[3]-region[2]), height=(region[1]-region[0]),
                                                  fc=color,
                                                  ec=ec,
                                                  lw=lw,
                                                  alpha=alpha,
                                                  transform=transform))
                else:
                    ax.add_patch(mpatches.Rectangle(xy=[region[2], region[0]], width=(region[3]-region[2]), height=(region[1]-region[0]),
                                                  fc=color,
                                                  lw=lw,
                                                  alpha=alpha,
                                                  ec=ec))
                if(name):
                    ax.text((region[3]+region[2])*0.5, (region[0]+region[1])*0.5, name,horizontalalignment='center',verticalalignment='center', transform=ccrs.Geodetic(),zorder=6)



            for region,values,difference in  toPlotData[["region","absolute","diff"]].to_numpy():

                if(values>countThreshold):
                    lat,lng = np.mean(region[0:1]),np.mean(region[2:3]);
                    colornorm = mpl.colors.Normalize(vmin=-values-0.5, vmax=values+0.5)
                    color=ccmap(colornorm(difference));
                    a = normOpacity(values)
    #                 color = np.array((0.43,0.43,0.43,1.0))*(1.0-a)+np.array(color)*a
                    color = np.power(np.power((0.43,0.43,0.43,1.0),0.45)*(1.0-a)+np.power(color,0.45)*a,1/0.45)

                    plotRectangle(ax,region,"",color,transform=ccrs.Geodetic(),alpha=1.0, ec=None,lw=0);

            plotRectangle(ax,[-5 , 5 ,-170, -120],"ENSO","#CCCCCC",transform=ccrs.Geodetic(),lw=0.2,ec=(0.1,0.1,0.1));


            ax.add_feature(cfeature.COASTLINE,edgecolor='black',alpha=0.5)
            ax.add_feature(cfeature.BORDERS,edgecolor='black',alpha=0.5)


            import json

            with open("../Data/ElNidoShapes.json","r") as fd:
                shapes = [shape[0] for shape in json.load(fd)];

            with open("../Data/LaNinaShapes.json","r") as fd:
                shapes += [shape[0] for shape in json.load(fd)];



            shapesData = [np.array(shape) for shape in shapes]  

            newPolygons = []
            for d in shapesData:
                newPolygons += split_polygon(geojson= {
                    "type": "Polygon",
                    "coordinates": [d]

                }, output_format= "polygons")
            newShape = unary_union(newPolygons)
            ax.add_patch(PolygonPatch(newShape, fc="white", ec=None,linewidth=0, alpha=0.2,transform=ccrs.Geodetic()))




            if(exportName):
                fig.savefig("../Figures/"+exportName+".png",dpi=300);
                fig.savefig("../Figures/"+exportName+".pdf");
#             plt.show()
            plt.close()



#plots the colorbar

# colornorm = mpl.colors.Normalize(vmin=-values-0.5, vmax=values+0.5)
# color=ccmap(colornorm(difference));

grid = np.random.rand(4, 4)

fig = plt.figure(figsize=(6,3))
ax = fig.gca()

maxValue = len(dfsAllData)
normOpacity = mpl.colors.Normalize(vmin=0-0.5, vmax=maxValue+0.5)
m = np.zeros((maxValue+1,2*maxValue+1,4))
for i in range(m.shape[0]):
    colornorm = mpl.colors.Normalize(vmin=-i-0.5, vmax=i+0.05)
    for j in range(m.shape[1]):
        
        color=ccmap(colornorm(j-maxValue));
        a = normOpacity(i)
        color = np.power(np.power((0.43,0.43,0.43,1.0),0.45)*(1.0-a)+np.power(color,0.45)*a,1/0.45)
        m[i,j] = color
        if(abs(j-maxValue)>i or i<5):
            m[i,j] = "nan"
ax.imshow(m,extent=(-maxValue-0.5, maxValue+0.5, maxValue+0.5, -0.5))
plt.savefig("../Figures/triangle_colormap.pdf")
plt.show()
# Plots the comparison among 90 days and 180 days maximum lag
nx = 288
ny = 140
# nx = 36
# ny = 14
dataRange = (0,100) #percentage
# 0
# 26
# 52
# 78
# 100
windowDays = 7
maxLag = 90
highResolution = (288,140)
lowResolution = (36,14)
countThreshold=0
dataRangeLags = [
#     ((0,26),90,highResolution),
#     ((26,52),90,highResolution),
#     ((52,78),90,highResolution),
#     ((78,100),90,highResolution),
#     ((0,100),90,highResolution),
#     ((0,100),"LaggedCorrelation",highResolution),
    ((0,100),"Causality",highResolution),
]


for dataRange,dataName,resolution in tqdm(dataRangeLags):
    nx=resolution[0]
    ny=resolution[1]
    percentile = 0.30
    pvalueThreshold = None
    smoothness = 0.5
    quantile = 0.5
    extraEpsilonFactor = 0.5
    min_samples =5
    min_points = 20
    exportName=""
    useCircular=False
    maxPositivePValue = None
    smoothness = 0.5

    extraEpsilonFactor = 0.5
    min_samples = 5
    min_points = 30

    for pvalueThreshold in [0.05,0.01,0.005,0.001]:
        exportName="%s_region_differences_7-30_l90x180_nx%d_ny%d_R%d-%d_p%.5f_v8"%(dataName,nx,ny,dataRange[0],dataRange[1],pvalueThreshold)
        maxPositivePValue=maxPositivePValue

        dfsAllDataCausality = []
        dfsAllDataCorrelation = []
        for index,windowDays in enumerate(tqdm(range(7,31))):
            filenameCausality = "../Data/%s_d%d_l%d_nx%d_ny%d_R%d-%d_v8.csv"%(dataName,windowDays,90,nx,ny,dataRange[0],dataRange[1])
            filenameCorrelation = "../Data/%s_d%d_l%d_nx%d_ny%d_R%d-%d_v8.csv"%(dataName,windowDays,180,nx,ny,dataRange[0],dataRange[1])
            dfCausality = pd.read_csv(filenameCausality, sep=',')
            dfCorrelation = pd.read_csv(filenameCorrelation, sep=',')
            if(index==0):
                gridRegions = [];
                for coords in dfCausality[["lat","lng"]].to_numpy():
                    region = (coords[0] - 70/ny, coords[0]+70/ny,coords[1]-180/nx,coords[1]+180/nx)
                    gridRegions.append(region)
            dfCausality["region"] = gridRegions;

            dfCorrelation["region"] = gridRegions;

            pointsCausality = dfCausality.reset_index(drop=True)
            pointsCorrelation = dfCorrelation.reset_index(drop=True)

            pointsCausality['log-pvalue'] = np.log10(pointsCausality['pvalue'])
            pointsCausality['thresh-pvalue'] = pointsCausality['pvalue']

            pointsCorrelation['log-pvalue'] = np.log10(pointsCorrelation['pvalue'])
            pointsCorrelation['thresh-pvalue'] = pointsCorrelation['pvalue']


            if(percentile<1.0):
                thresholdCausality = np.quantile(pointsCausality[['pvalue']],percentile);
            else:
                thresholdCausality = 1.0
            if(pvalueThreshold is not None):
                thresholdCausality = pvalueThreshold;


            if(percentile<1.0):
                thresholdCorrelation = np.quantile(pointsCorrelation[['pvalue']],percentile);
            else:
                thresholdCorrelation = 1.0
            if(pvalueThreshold is not None):
                thresholdCorrelation = pvalueThreshold;


            pointsCausality.loc[pointsCausality['thresh-pvalue']>thresholdCausality,'thresh-pvalue']=1
            pointsCorrelation.loc[pointsCorrelation['thresh-pvalue']>thresholdCorrelation,'thresh-pvalue']=1

            dfsAllDataCausality.append(pointsCausality)
            dfsAllDataCorrelation.append(pointsCorrelation)

        countsCausality = np.zeros(len(dfsAllDataCausality[0]))
        numberOfDays = len(dfsAllDataCausality)
        for df in dfsAllDataCausality:
            countsCausality[df["thresh-pvalue"]<=thresholdCausality] +=1
        #     counts += df["log-pvalue"]/numberOfDays

        countsCorrelation = np.zeros(len(dfsAllDataCorrelation[0]))
        numberOfDays = len(dfsAllDataCorrelation)
        for df in dfsAllDataCorrelation:
            countsCorrelation[df["thresh-pvalue"]<=thresholdCorrelation] +=1
        #     counts += df["log-pvalue"]/numberOfDays


        dfsAllDataCausality[0]["counts"] = countsCausality
        dfsAllDataCorrelation[0]["counts"] = countsCorrelation
        

        saveFilenameCausality = "../Data/%s_dAll_l%d_nx%d_ny%d_R%d-%d_p%.5f_v8.csv"%("Causality",maxLag,nx,ny,dataRange[0],dataRange[1],pvalueThreshold)
        dfsAllDataCorrelation[0].to_csv(saveFilenameCausality,sep=",")

        saveFilenameCorrelation = "../Data/%s_dAll_l%d_nx%d_ny%d_R%d-%d_p%.5f_v8.csv"%("LaggedCorrelation",maxLag,nx,ny,dataRange[0],dataRange[1],pvalueThreshold)
        dfsAllDataCorrelation[0].to_csv(saveFilenameCorrelation,sep=",")

        dfsAllDataCausality[0]["diff"] = countsCausality-countsCorrelation
        dfsAllDataCausality[0]["absolute"] = np.max([countsCausality,countsCorrelation],axis=0)

        import copy
        toPlotData = dfsAllDataCausality[0]
        tab20b = cm.get_cmap('tab20b', 20)
        tab20c = cm.get_cmap('tab20c', 20)
        newcmp = mpc.ListedColormap(tab20b(np.linspace(0.25, 0.75, 256)))
        a = np.concatenate([tab20b(range(20)),tab20c(range(20))])
        np.random.shuffle(a)
        newcmp = mpc.ListedColormap(a)
        from sklearn.cluster import DBSCAN
        def identity_scale(minval, maxval):
            def scalar(val):
                return 2
            return scalar

        useCircular= True
        ENSORegion = [[-5 , -170] ,[5, -120]];

        cdict1 = {'red':   ((0.0, 0.0, 0.0),
                           (0.5, 0.5, 1.0),
                           (1.0, 1.0, 0.0)),

                 'green': ((0.0, 0.0, 0.0),
                           (0.5, 1.0, 1.0),
                           (1.0, 0.0, 0.0)),

                 'blue':  ((0.0, 0.0, 1.0),
                           (0.5, 1.0, 0.0),
                           (1.0, 0.0, 0.0))
                }

        vmin = -len(dfsAllDataCausality)
        vmax = len(dfsAllDataCausality)
        norm = mpl.colors.Normalize(vmin=vmin-0.5, vmax=vmax+0.5)
        normOpacity = mpl.colors.Normalize(vmin=-.05, vmax=vmax+0.5)

        ccmap = copy.copy(matplotlib.cm.get_cmap("PiYG_r",vmax-vmin+1))

        if(True):
            if(useCircular):
                projection = ccrs.Robinson(central_longitude=150)
            else:
                projection = ccrs.PlateCarree(central_longitude=150);

            fig = plt.figure(figsize=(10,8),tight_layout=True)
            ax = plt.axes(projection=projection)
            # ax.set_extent([-30, 150, -70, 70])
            # ax.stock_img()
            ax.set_global()
        #         if(not useCircular):
        #             ax.set_extent([-30, 150, -70, 70])
            ax.add_feature(cfeature.LAND,facecolor='gray',alpha=1.0)
            ax.add_feature(cfeature.OCEAN,facecolor=(0.7,0.7,0.7),alpha=1.0)
            # ax.add_feature(cfeature.OCEAN)

          #Arrows lines 

            transform = projection._as_mpl_transform(ax)

            def plotRectangle(ax,region,name,color,transform = ccrs.Geodetic(),lw=0,ec=None,alpha=1.0):
                if(transform):
                    ax.add_patch(mpatches.Rectangle(xy=[region[2], region[0]], width=(region[3]-region[2]), height=(region[1]-region[0]),
                                                  fc=color,
                                                  ec=ec,
                                                  lw=lw,
                                                  alpha=alpha,
                                                  transform=transform))
                else:
                    ax.add_patch(mpatches.Rectangle(xy=[region[2], region[0]], width=(region[3]-region[2]), height=(region[1]-region[0]),
                                                  fc=color,
                                                  lw=lw,
                                                  alpha=alpha,
                                                  ec=ec))
                if(name):
                    ax.text((region[3]+region[2])*0.5, (region[0]+region[1])*0.5, name,horizontalalignment='center',verticalalignment='center', transform=ccrs.Geodetic(),zorder=6)



            for region,values,difference in  toPlotData[["region","absolute","diff"]].to_numpy():

                if(values>countThreshold):
                    lat,lng = np.mean(region[0:1]),np.mean(region[2:3]);
                    colornorm = mpl.colors.Normalize(vmin=-values-0.5, vmax=values+0.5)
                    color=ccmap(colornorm(difference));
                    a = normOpacity(values)
    #                 color = np.array((0.43,0.43,0.43,1.0))*(1.0-a)+np.array(color)*a
                    color = np.power(np.power((0.43,0.43,0.43,1.0),0.45)*(1.0-a)+np.power(color,0.45)*a,1/0.45)

                    plotRectangle(ax,region,"",color,transform=ccrs.Geodetic(),alpha=1.0, ec=None,lw=0);

            plotRectangle(ax,[-5 , 5 ,-170, -120],"ENSO","#CCCCCC",transform=ccrs.Geodetic(),lw=0.2,ec=(0.1,0.1,0.1));


            ax.add_feature(cfeature.COASTLINE,edgecolor='black',alpha=0.5)
            ax.add_feature(cfeature.BORDERS,edgecolor='black',alpha=0.5)


            import json

            with open("../Data/ElNidoShapes.json","r") as fd:
                shapes = [shape[0] for shape in json.load(fd)];

            with open("../Data/LaNinaShapes.json","r") as fd:
                shapes += [shape[0] for shape in json.load(fd)];



            shapesData = [np.array(shape) for shape in shapes]  

            newPolygons = []
            for d in shapesData:
                newPolygons += split_polygon(geojson= {
                    "type": "Polygon",
                    "coordinates": [d]

                }, output_format= "polygons")
            newShape = unary_union(newPolygons)
            ax.add_patch(PolygonPatch(newShape, fc="white", ec=None,linewidth=0, alpha=0.2,transform=ccrs.Geodetic()))




            if(exportName):
                fig.savefig("../Figures/"+exportName+".png",dpi=300);
#                 fig.savefig("../Figures/"+exportName+".pdf"); #enable to plot pdfs
#             plt.show()
            plt.close()
#plots the colormap

grid = np.random.rand(4, 4)

fig = plt.figure(figsize=(6,3))
ax = fig.gca()

maxValue = len(dfsAllDataCausality)
normOpacity = mpl.colors.Normalize(vmin=0-0.5, vmax=maxValue+0.5)
m = np.zeros((maxValue+1,2*maxValue+1,4))
for i in range(m.shape[0]):
    colornorm = mpl.colors.Normalize(vmin=-i-0.5, vmax=i+0.05)
    for j in range(m.shape[1]):
        
        color=ccmap(colornorm(j-maxValue));
        a = normOpacity(i)
        color = np.power(np.power((0.43,0.43,0.43,1.0),0.45)*(1.0-a)+np.power(color,0.45)*a,1/0.45)
        m[i,j] = color
        if(abs(j-maxValue)>i or i<5):
            m[i,j] = "nan"
ax.imshow(m,extent=(-maxValue-0.5, maxValue+0.5, maxValue+0.5, -0.5))
plt.savefig("../Figures/triangle_colormap_90x180.pdf")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14253606076986544, 'fraction_numerical': 0.040228740504708524, 'mean_word_length': 2.321050642479214, 'pattern_counts': {'":': 4, '<': 25, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 58, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '335774', 'license': 'MIT', 'n_tokens_mistral': 12782, 'n_tokens_neox': 11847, 'n_words': 1726, 'path': 'Notebooks/DrawMaps.ipynb', 'repo_name': 'filipinascimento/teleconectionsgranger', 'size': '35149'}