使用python制造渗流模型

简介

关于渗流模型的简介我们可以看这个页面：渗流模型

构建最基本的渗流模型

调用包并设定参数

    from pylab import *
    from scipy.ndimage import measurements
    import statsmodels.api as sm

    L = 100
    r = rand(L,L)
    p = 0.4
    z = r < p

绘图

    imshow(z, origin='lower', interpolation='nearest')
    colorbar()
    title("Matrix")

标记不同的cluster

    lw, num = measurements.label(z)
    imshow(lw, origin='lower', interpolation='nearest')
    colorbar()
    title("Labeled clusters")

可以随机洗牌cluster label使得颜色更清楚

    b = arange(lw.max() + 1)
    shuffle(b)
    shuffledLw = b[lw]
    imshow(shuffledLw, origin='lower', interpolation='nearest')
    colorbar()
    title("Labeled clusters")

还可以令cluster的颜色与大小正比

    area = measurements.sum(z, lw, index=arange(lw.max() + 1))
    areaImg = area[lw]
    im3 = imshow(areaImg, origin='lower', interpolation='nearest')
    colorbar()
    title("Clusters by area")

给最大cluster加一个方框

500px

    im3 = imshow(areaImg, origin='lower', interpolation='nearest')
    colorbar()
    title("Clusters by area")
    sliced = measurements.find_objects(areaImg == areaImg.max())
    if(len(sliced) > 0):
        sliceX = sliced[0][1]
        sliceY = sliced[0][0]
        plotxlim=im3.axes.get_xlim()
        plotylim=im3.axes.get_ylim()
        plot([sliceX.start, sliceX.start, sliceX.stop, sliceX.stop, sliceX.start], \
                          [sliceY.start, sliceY.stop, sliceY.stop, sliceY.start, sliceY.start], \
                          color="red")
        xlim(plotxlim)
        ylim(plotylim)

拟合cluster size的分布

    def RankOrderPlot(data):
        d=array(data)
        d = d[d>0]
        t=array(sorted(d,key=lambda x:-x))
        r=array(range(1,len(d)+1))   
        y = log(t)
        x = log(r)
        X = sm.add_constant(x, prepend=True)
        res = sm.OLS(y,X).fit()
        C,beta = res.params
        plt.plot(r,t,"o",color="b")
        plt.plot(r,exp(C)*r**beta,"r-")
        plt.xscale('log')
        plt.yscale('log')
        plt.text(max(r)/2,max(t)/100,"beta = " + str(round(beta,2)))
        print beta
         
    RankOrderPlot(area)