In [1]:

    

import sys, os
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import matplotlib.colors as mc
%matplotlib inline


    

In [2]:

    

palette = [
    "#000000", "#0019ff", "#0080ff",
    "#00e5ff", "#00ff4d", "#19fe00",
    "#80ff00", "#e6ff00", "#ffff00",
    "#ffe53c", "#ffdb77", "#ffe0b2",
]
n = len(palette)


    

In [3]:

    

cmap = mc.ListedColormap(palette, 'indexed')
im = plt.imshow([range(n)], interpolation='none', cmap=cmap)
plt.show()


    

In [4]:

    

class Raster(np.ndarray):
    def __new__(cls, input_array, cellsize, xllcorner, yllcorner):
        obj = np.asarray(input_array).view(cls)
        obj.cellsize  = cellsize
        obj.xllcorner = xllcorner
        obj.yllcorner = yllcorner
        
        # 1-D list of the raster's values 
        v = obj.flatten()
        obj.values = v[~np.isnan(v)]
        
        return obj
    
    @classmethod
    def load(cls, filename):
        data  = [ ]
        attrs = { }
        with open(filename) as f:
            for i in range(6):
                key,value = next(f).strip().split()
                attrs[key] = value
            NODATA_value = attrs.get("NODATA_value", "-9999")
        
            for line in f:
                row = list(map(lambda x: np.nan if x == NODATA_value else float(x), line.split()))
                data.append(row)

        return Raster(np.asarray(data), float(attrs["cellsize"])
                                      , float(attrs["xllcorner"])
                                      , float(attrs["yllcorner"]))


    

In [5]:

    

%%time
I = Raster.load("../bb450_final.asc")


    

    




CPU times: user 3.98 s, sys: 293 ms, total: 4.27 s
Wall time: 5.71 s

In [6]:

    

lo    = float(np.min(I.values))
hi    = float(np.max(I.values))
mu    = float(np.mean(I.values))
sigma = float(np.std(I.values))


    

In [7]:

    

print("""
Size:       {} x {}
Cell size:  {}
Min Value:  {}
Max Value:  {}
Mean Value: {}
Std Dev:    {}
Colors:     {}
""".format(I.shape[0], I.shape[1], I.cellsize, lo, hi, mu, sigma, n))


    

    




Size:       2845 x 2483
Cell size:  450.0
Min Value:  0.0
Max Value:  7.128635
Mean Value: 0.13578042048469882
Std Dev:    0.08593180316938359
Colors:     12

In [8]:

    

def plot_histogram(I, divs):
    # print(list(zip(divs, palette)))
    cmap = mc.LinearSegmentedColormap.from_list("", list(zip(divs, palette)))
    nbars  = 300
    width  = (hi - lo) / nbars
    h,x    = np.histogram(I.values - width/2., nbars)  # shift the bars slightly, so they are centered over x coord
    colors = list(map(cmap, x[1:] / hi))  # data values must be [0.,1.]

    plt.figure(figsize=(12, 5))
    ax = plt.bar(x[1:], h, width, color=colors, edgecolor=colors)
    for i in divs:
        plt.axvline(x=i*hi, linestyle="--", color="#cecece")

    plt.yscale("symlog")
    plt.ylabel("Count")
    plt.xlabel("Amps")
    plt.xlim([-0.1,None])


    

In [9]:

    

def plot_map(I, divs):
    cmap = mc.LinearSegmentedColormap.from_list("", list(zip(divs, palette)))
    fig = plt.figure(figsize=(10,10))
    img = plt.imshow(I, origin="upper", cmap=cmap)
    fig.colorbar(img)


    

In [10]:

    

divs = [ i/(n-1.) for i in range(n) ]
divs


    

    
Out[10]:





[0.0,
 0.09090909090909091,
 0.18181818181818182,
 0.2727272727272727,
 0.36363636363636365,
 0.45454545454545453,
 0.5454545454545454,
 0.6363636363636364,
 0.7272727272727273,
 0.8181818181818182,
 0.9090909090909091,
 1.0]

In [11]:

    

plot_histogram(I, divs)


    

In [12]:

    

plot_map(I, divs)


    

In [13]:

    

c = 1.75
divs = [0.] + [ 1/c**(n-i-1) for i in range(1,n) ]
divs


    

    
Out[13]:





[0.0,
 0.003712098683732818,
 0.006496172696532431,
 0.011368302218931755,
 0.01989452888313057,
 0.0348154255454785,
 0.06092699470458737,
 0.10662224073302791,
 0.18658892128279883,
 0.32653061224489793,
 0.5714285714285714,
 1.0]

In [14]:

    

plot_histogram(I, divs)


    

In [15]:

    

plot_map(I, divs)


    

In [16]:

    

divs = [0.] + [ (i*sigma+mu)/hi for i in range(1,n-1) ] + [1.]
divs


    

    
Out[16]:





[0.0,
 0.031101637782560393,
 0.04315609184976731,
 0.055210545916974225,
 0.06726499998418116,
 0.07931945405138806,
 0.09137390811859498,
 0.1034283621858019,
 0.11548281625300882,
 0.12753727032021575,
 0.13959172438742265,
 1.0]

In [17]:

    

plot_histogram(I, divs)


    

In [18]:

    

plot_map(I, divs)


    

In [19]:

    

divs = [0.] + [ a[-1]/hi for a in np.array_split(np.sort(I.values), n-1) ]
divs


    

    
Out[19]:





[0.0,
 0.0068480150828314256,
 0.010583372553090458,
 0.012930946808189786,
 0.014841971850150834,
 0.016628288585402395,
 0.018434244424072772,
 0.020363084938420893,
 0.02266675176944815,
 0.025702957157997287,
 0.03107284914994245,
 1.0]

In [20]:

    

plot_histogram(I, divs)


    

In [21]:

    

plot_map(I, divs)


    

In [ ]: