source: trunk/src/ffgrid2_gc_nof_pos.py @ 56

#modificado para grilla irregular griddata.py
import numpy as np
from numpy import *

x=xx
y=yy
z=zz
dx=0.1
dy=1
dparx=0.05
dpary=0.5
n0=0
z0=100
z1=300
y11= -50
x0, x1 = -180, 180
y0, y1 = -90, 90
dx=0.1
dy=1.0

n=size(lon)
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

#save txt file
#np.savetxt('ix_py.txt', ix, fmt='%.18e', delimiter=';')


inx = (ix >= 0) & (ix < nx)
iny = (iy >= 0) & (iy < ny)
inz = (z >= z0) & (z <= z1)
ipos = ((pos >=11) & (pos<=20)) #| ((pos >=21) & (pos <=25))
inn = inx & iny & inz & ipos
iix = ix[inn]
iiy = iy[inn]
zz = z[inn]

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

#calculation of the each arc legth
Rt=6371
beta=math.pi/180
R=beta*Rt
L=zeros(ny,float)#arc vector
L[0]=2*math.pi*R

for i in range(1, ny):
    L[i]=2*math.pi*R*(i+1)

nnx=zeros(ny,float)
nnx=L/L[0]
nnx=200/nnx
nnx=nnx.round()

div=(1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,30,36,40,45,48,50,60,
72,75,80,90,100,120,144,150,180,200,225,240,300,360,400,450,
600,720,900,1200,1800,3600)
vnx=zeros(ny,float)

for i in range(0, ny):
    a=nnx[i]
    for j in range(0, 45):
        b=div[j]
        if b<=a:
            vnx[i]=b

vnx[13:40]=12
li=(L/3600)*vnx
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
z2grid = zeros([nx, ny], float)# z-coordinate ²
zz2grid = zeros([nx, ny], float)# z-coordinate ² grouped

for i in range(0, N):
    zgrid[iix[i], iiy[i]] = zgrid[iix[i],iiy[i]]+zz[i]
    z2grid[iix[i], iiy[i]] = z2grid[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(z2grid[:,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)
zzgrid = np.transpose(zzgrid)
zzpgrid = np.transpose(zzpgrid)
zz2grid = np.transpose(zz2grid)
ngrid = np.transpose(ngrid)
nngrid = np.transpose(nngrid)
sigma_grid=np.transpose(sigma_grid)