source: trunk/src/ffgrid2.py @ 56

# griddata.py - 2010-07-11 ccampo
import numpy as np
from numpy import *
import bina

def ffgrid(x, y, z, dx, dy, x0, x1,y0, y1, z0, z1):
#def ffgrid(x, y, z, dx, dy, x0, x1,y0, y1, z0, z1):
    # [zgrid, xvec, yvec] = ffgrid(x, y, z, dx, dy, [x0 x1 y0 y1 z0 z1 n0]);
    #
    #  Fast 'n' Furious data gridding.
    #
    #  Input:  Unevenly spaced data in vectors x, y, and z of equal length.
    #          If dx and dy are omitted, the default is 75 bins in each direction.
    #
    #          If dx or dy is negative, then the variable is taken as the number of
    #          bins rather than a grid resolution. If dx is complex, the imaginary
    #          part is taken as the value with which empty grid points are padded.
    #          The default padding is 10% below minimum value in grid.
    #
    #          The vector containing the limits can be padded with NaNs if only
    #          certain limits are desired, e g if x1 and z0 are wanted:
    #
    #            ffgrid(x, y, z, [.5 nan nan nan 45])
    #
    #          The last parameter, n0, removes outliers from the data set by
    #          ignoring grid points with n0 or less observations. When n0 is
    #          negative it is treated as the percentage of the total number of
    #          data points.
    #
    #  Output: The gridded matrix ZGRID together with vectors XVEC and YVEC.
    #          If no output arguments are given, FFGRID will plot the gridded
    #          function with the prescribed axes using PCOLOR.
    #
    #  Requires bin.m. Tested under MatLab 4.2, 5.0, and 5.1.
    #
    #  See also bin.m for further details about dx, x0, etc, and density.m.
    #

    # 28.7.97, Oyvind.Breivik@gfi.uib.no.
    #
    # Oyvind Breivik
    # Department of Geophysics
    # University of Bergen
    # NORWAY

    dpar=0.5
    n0=0

    xvec = arange(x0,x1+dx*dpar,dx)
    yvec = arange(y0,y1+dy*dpar,dy)

    #dx=(x1-x0)/100
    #dy=(y1-y0)/50

    ix = bina.bina(x, dx, x0, x1+dx*dpar)
    iy = bina.bina(y, dy, y0, y1+dy*dpar)



    print ix[0:5]

   # bin data in (x,y)-space

    #xvec = arange(x0,x1+dx/2,dx);
    #yvec = arange(y0,y1+dy/2,dy);

    #xvec = linspace(x0,x1,100)
    #yvec = linspace(y0,y1,50)
    nx = size(xvec)
    ny = size(yvec)

    inx = (ix >= 1) & (ix <= nx)
    iny = (iy >= 1) & (iy <= ny)
    inz = (z >= z0) & (z <= z1)
    inn = inx & iny & inz
    ix = ix[inn]
    iy = iy[inn]

    print ix[0:5]
    z = z[inn]

    N = size(ix)# how many datapoints are left now?

    ngrid = zeros([nx, ny], float)# no of obs per grid cell
    zgrid = zeros([nx, ny], float)# z-coordinate

    for i in range(0, N):
        zgrid[ix[i]-1, iy[i]-1] = zgrid[ix[i]-1,iy[i]-1]+z[i]
        ngrid[ix[i]-1, iy[i]-1] = ngrid[ix[i]-1,iy[i]-1]+ 1


    Nil = (ngrid == 0)
    ngrid[Nil] = 1# now we don't divide by zero
    zgrid = zgrid /ngrid
    zgrid[Nil] = NaN

    N = sum(sum(ngrid))# how many datapoints are left now?

    #Nil = (ngrid == 0)

        #ngrid[Nil] = 1# now we don't divide by zero

        #for i in range(len(ix)): zgrid[ix[i]-1, iy[i]-1] = zgrid[ix[i]-1, iy[i]-1]/ngrid[ix[i]-1, iy[i]-1]

        #zgrid = linalg.solve(zgrid,ngrid)
        #zgrid = zgrid /ngrid# Make average of values for each point

        #c = np.linalg.lstsq(ngrid.T, zgrid.T)[0].T

        #zgrid[Nil] = NaN# Empty grid points are set to default value
    #rho = nnz(not Nil) / length(Nil(mslice[:]))

    zgrid = np.transpose(zgrid)
    ngrid = np.transpose(ngrid)

    #dum = size(zgrid.cT)

    zzgrid = zgrid
    xxvec = xvec
    yyvec = yvec

    return zgrid, xvec, yvec