[4] | 1 | # griddata.py - 2010-07-11 ccampo |
---|
| 2 | import numpy as np |
---|
| 3 | from numpy import * |
---|
| 4 | import bina |
---|
| 5 | |
---|
| 6 | def ffgrid(x, y, z, dx, dy, x0, x1,y0, y1, z0, z1): |
---|
| 7 | #def ffgrid(x, y, z, dx, dy, x0, x1,y0, y1, z0, z1): |
---|
| 8 | # [zgrid, xvec, yvec] = ffgrid(x, y, z, dx, dy, [x0 x1 y0 y1 z0 z1 n0]); |
---|
| 9 | # |
---|
| 10 | # Fast 'n' Furious data gridding. |
---|
| 11 | # |
---|
| 12 | # Input: Unevenly spaced data in vectors x, y, and z of equal length. |
---|
| 13 | # If dx and dy are omitted, the default is 75 bins in each direction. |
---|
| 14 | # |
---|
| 15 | # If dx or dy is negative, then the variable is taken as the number of |
---|
| 16 | # bins rather than a grid resolution. If dx is complex, the imaginary |
---|
| 17 | # part is taken as the value with which empty grid points are padded. |
---|
| 18 | # The default padding is 10% below minimum value in grid. |
---|
| 19 | # |
---|
| 20 | # The vector containing the limits can be padded with NaNs if only |
---|
| 21 | # certain limits are desired, e g if x1 and z0 are wanted: |
---|
| 22 | # |
---|
| 23 | # ffgrid(x, y, z, [.5 nan nan nan 45]) |
---|
| 24 | # |
---|
[9] | 25 | # The last parameter, n0, removes outliers from the data set by |
---|
| 26 | # ignoring grid points with n0 or less observations. When n0 is |
---|
| 27 | # negative it is treated as the percentage of the total number of |
---|
[4] | 28 | # data points. |
---|
| 29 | # |
---|
[9] | 30 | # Output: The gridded matrix ZGRID together with vectors XVEC and YVEC. |
---|
[4] | 31 | # If no output arguments are given, FFGRID will plot the gridded |
---|
| 32 | # function with the prescribed axes using PCOLOR. |
---|
| 33 | # |
---|
| 34 | # Requires bin.m. Tested under MatLab 4.2, 5.0, and 5.1. |
---|
| 35 | # |
---|
| 36 | # See also bin.m for further details about dx, x0, etc, and density.m. |
---|
| 37 | # |
---|
| 38 | |
---|
| 39 | # 28.7.97, Oyvind.Breivik@gfi.uib.no. |
---|
| 40 | # |
---|
| 41 | # Oyvind Breivik |
---|
| 42 | # Department of Geophysics |
---|
| 43 | # University of Bergen |
---|
| 44 | # NORWAY |
---|
[9] | 45 | |
---|
[4] | 46 | dpar=0.5 |
---|
| 47 | n0=0 |
---|
[9] | 48 | |
---|
[4] | 49 | xvec = arange(x0,x1+dx*dpar,dx) |
---|
| 50 | yvec = arange(y0,y1+dy*dpar,dy) |
---|
[9] | 51 | |
---|
[4] | 52 | #dx=(x1-x0)/100 |
---|
| 53 | #dy=(y1-y0)/50 |
---|
[9] | 54 | |
---|
[4] | 55 | ix = bina.bina(x, dx, x0, x1+dx*dpar) |
---|
| 56 | iy = bina.bina(y, dy, y0, y1+dy*dpar) |
---|
[9] | 57 | |
---|
| 58 | |
---|
| 59 | |
---|
[4] | 60 | print ix[0:5] |
---|
[9] | 61 | |
---|
[4] | 62 | # bin data in (x,y)-space |
---|
| 63 | |
---|
| 64 | #xvec = arange(x0,x1+dx/2,dx); |
---|
| 65 | #yvec = arange(y0,y1+dy/2,dy); |
---|
[9] | 66 | |
---|
[4] | 67 | #xvec = linspace(x0,x1,100) |
---|
| 68 | #yvec = linspace(y0,y1,50) |
---|
| 69 | nx = size(xvec) |
---|
| 70 | ny = size(yvec) |
---|
| 71 | |
---|
[9] | 72 | inx = (ix >= 1) & (ix <= nx) |
---|
[4] | 73 | iny = (iy >= 1) & (iy <= ny) |
---|
[9] | 74 | inz = (z >= z0) & (z <= z1) |
---|
[4] | 75 | inn = inx & iny & inz |
---|
| 76 | ix = ix[inn] |
---|
| 77 | iy = iy[inn] |
---|
[9] | 78 | |
---|
[4] | 79 | print ix[0:5] |
---|
| 80 | z = z[inn] |
---|
| 81 | |
---|
| 82 | N = size(ix)# how many datapoints are left now? |
---|
| 83 | |
---|
| 84 | ngrid = zeros([nx, ny], float)# no of obs per grid cell |
---|
| 85 | zgrid = zeros([nx, ny], float)# z-coordinate |
---|
[9] | 86 | |
---|
[4] | 87 | for i in range(0, N): |
---|
| 88 | zgrid[ix[i]-1, iy[i]-1] = zgrid[ix[i]-1,iy[i]-1]+z[i] |
---|
| 89 | ngrid[ix[i]-1, iy[i]-1] = ngrid[ix[i]-1,iy[i]-1]+ 1 |
---|
| 90 | |
---|
[9] | 91 | |
---|
[4] | 92 | Nil = (ngrid == 0) |
---|
| 93 | ngrid[Nil] = 1# now we don't divide by zero |
---|
| 94 | zgrid = zgrid /ngrid |
---|
| 95 | zgrid[Nil] = NaN |
---|
[9] | 96 | |
---|
[4] | 97 | N = sum(sum(ngrid))# how many datapoints are left now? |
---|
| 98 | |
---|
| 99 | #Nil = (ngrid == 0) |
---|
| 100 | |
---|
[9] | 101 | #ngrid[Nil] = 1# now we don't divide by zero |
---|
[4] | 102 | |
---|
[9] | 103 | #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] |
---|
[4] | 104 | |
---|
[9] | 105 | #zgrid = linalg.solve(zgrid,ngrid) |
---|
| 106 | #zgrid = zgrid /ngrid# Make average of values for each point |
---|
| 107 | |
---|
| 108 | #c = np.linalg.lstsq(ngrid.T, zgrid.T)[0].T |
---|
| 109 | |
---|
| 110 | #zgrid[Nil] = NaN# Empty grid points are set to default value |
---|
[4] | 111 | #rho = nnz(not Nil) / length(Nil(mslice[:])) |
---|
| 112 | |
---|
| 113 | zgrid = np.transpose(zgrid) |
---|
| 114 | ngrid = np.transpose(ngrid) |
---|
| 115 | |
---|
| 116 | #dum = size(zgrid.cT) |
---|
| 117 | |
---|
| 118 | zzgrid = zgrid |
---|
| 119 | xxvec = xvec |
---|
| 120 | yyvec = yvec |
---|
[9] | 121 | |
---|
[4] | 122 | return zgrid, xvec, yvec |
---|