source: trunk/src/python_script/newgrid1.py @ 6

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
from numpy import *
import bina

def newgrid(x, y, z, vnx):


    dx=0.1
    dy=1        
    dparx=0.05
    dpary=0.5
    n0=0
    z0=0
    z1=1
    y11= -50
    x0, x1 = -180, 180
    y0, y1 = -90, 90
    dx=0.1
    dy=1.0
    n=size(x)    
    xvec = arange(x0,x1,dx)
    yvec = arange(y0,y11,dy)
    nx = size(xvec)
    ny = size(yvec)
    ix = zeros(n, float)    
    iy = zeros(n, float)    

    for kk in range(0,n):
        if x[kk] == -180:
           ix[kk] = 3600
        else:
           ix[kk] = math.ceil((x[kk] - x0)/dx)-1

    for kk in range(0,n):
        if y[kk] == -90:
           iy[kk] = 0
        else:
           iy[kk] = math.ceil((y[kk] - y0)/dy)-1

    inx = (ix >= 0) & (ix < nx) 
    iny = (iy >= 0) & (iy < ny)
    inz = (z >= z0) & (z <= z1) 
    inn = inx & iny & inz
    iix = ix[inn]
    iiy = iy[inn]
    zz = z[inn]
    N = size(iix)# how many datapoints are left now?
    #calculation of the each arc legth

    dxx=vnx
    ngrid = zeros([nx, ny], float)# no of obs per grid cell
    nngrid = zeros([nx, ny], float)
    zgrid = zeros([nx, ny], float)# z-coordinate
    zzgrid = zeros([nx, ny], float)# z-coordinate grouped
    zdgrid = zeros([nx, ny], float)# z-coordinate 2
    zz2grid = zeros([nx, ny], float)# z-coordinate 2 grouped

    for i in range(0, N):
        zgrid[iix[i], iiy[i]] = zgrid[iix[i],iiy[i]]+zz[i]
        zdgrid[iix[i], iiy[i]] = zdgrid[iix[i],iiy[i]]+zz[i]*zz[i]
        ngrid[iix[i], iiy[i]] = ngrid[iix[i],iiy[i]]+ 1 
        
    for i in range(0,ny):
        jj=round(3600/dxx[i])
        c=reshape(ngrid[:,i],(jj,dxx[i]))
        b=c.sum(axis=1)
        d=reshape(zgrid[:,i],(jj,dxx[i]))
        e=d.sum(axis=1)
        f=reshape(zdgrid[:,i],(jj,dxx[i]))
        g=f.sum(axis=1)  
        for j in range(0, nx):
            ii=ceil((j+1)/dxx[i])-1
            nngrid[j,i]=b[ii]
            zzgrid[j,i]=e[ii]
            zz2grid[j,i]=g[ii]

    zzpgrid = zeros([nx, ny], float)    
    zzpgrid = zzgrid/nngrid
    sigma_grid=sqrt(zz2grid/nngrid-zzpgrid*zzpgrid)
    zgrid = np.transpose(zgrid)
    sigma_grid=np.transpose(sigma_grid)
    zzpgrid = np.transpose(zzpgrid)
    ngrid = np.transpose(ngrid)
    nngrid = np.transpose(nngrid)
    return zgrid, zzpgrid, ngrid, nngrid, sigma_grid, xvec, yvec, zz