1 | #modificado para grilla irregular griddata.py |
---|
2 | import numpy as np |
---|
3 | from numpy import * |
---|
4 | |
---|
5 | |
---|
6 | |
---|
7 | x=xx |
---|
8 | y=yy |
---|
9 | z=zz |
---|
10 | dx=0.1 |
---|
11 | dy=1 |
---|
12 | dparx=0.05 |
---|
13 | dpary=0.5 |
---|
14 | n0=0 |
---|
15 | z0=100 |
---|
16 | z1=300 |
---|
17 | y11= -50 |
---|
18 | x0, x1 = -180, 180 |
---|
19 | y0, y1 = -90, 90 |
---|
20 | dx=0.1 |
---|
21 | dy=1.0 |
---|
22 | |
---|
23 | n=size(lon) |
---|
24 | xvec = arange(x0,x1,dx) |
---|
25 | yvec = arange(y0,y11,dy) |
---|
26 | nx = size(xvec) |
---|
27 | ny = size(yvec) |
---|
28 | |
---|
29 | ix = zeros(n, float) |
---|
30 | iy = zeros(n, float) |
---|
31 | |
---|
32 | for kk in range(0,n): |
---|
33 | if x[kk] == -180: |
---|
34 | ix[kk] = 3600 |
---|
35 | else: |
---|
36 | ix[kk] = math.ceil((x[kk] - x0)/dx)-1 |
---|
37 | |
---|
38 | for kk in range(0,n): |
---|
39 | if y[kk] == -90: |
---|
40 | iy[kk] = 0 |
---|
41 | else: |
---|
42 | iy[kk] = math.ceil((y[kk] - y0)/dy)-1 |
---|
43 | |
---|
44 | #save txt file |
---|
45 | #np.savetxt('ix_py.txt', ix, fmt='%.18e', delimiter=';') |
---|
46 | |
---|
47 | |
---|
48 | inx = (ix >= 0) & (ix < nx) |
---|
49 | iny = (iy >= 0) & (iy < ny) |
---|
50 | inz = (z >= z0) & (z <= z1) |
---|
51 | ipos = ((pos >=11) & (pos<=20)) #| ((pos >=21) & (pos <=25)) |
---|
52 | inn = inx & iny & inz & ipos |
---|
53 | iix = ix[inn] |
---|
54 | iiy = iy[inn] |
---|
55 | zz = z[inn] |
---|
56 | |
---|
57 | N = size(iix)# how many datapoints are left now? |
---|
58 | |
---|
59 | #calculation of the each arc legth |
---|
60 | Rt=6371 |
---|
61 | beta=math.pi/180 |
---|
62 | R=beta*Rt |
---|
63 | L=zeros(ny,float)#arc vector |
---|
64 | L[0]=2*math.pi*R |
---|
65 | |
---|
66 | for i in range(1, ny): |
---|
67 | L[i]=2*math.pi*R*(i+1) |
---|
68 | |
---|
69 | nnx=zeros(ny,float) |
---|
70 | nnx=L/L[0] |
---|
71 | nnx=200/nnx |
---|
72 | nnx=nnx.round() |
---|
73 | |
---|
74 | div=(1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,30,36,40,45,48,50,60, |
---|
75 | 72,75,80,90,100,120,144,150,180,200,225,240,300,360,400,450, |
---|
76 | 600,720,900,1200,1800,3600) |
---|
77 | vnx=zeros(ny,float) |
---|
78 | |
---|
79 | for i in range(0, ny): |
---|
80 | a=nnx[i] |
---|
81 | for j in range(0, 45): |
---|
82 | b=div[j] |
---|
83 | if b<=a: |
---|
84 | vnx[i]=b |
---|
85 | |
---|
86 | vnx[13:40]=12 |
---|
87 | li=(L/3600)*vnx |
---|
88 | dxx=vnx |
---|
89 | |
---|
90 | ngrid = zeros([nx, ny], float)# no of obs per grid cell |
---|
91 | nngrid = zeros([nx, ny], float) |
---|
92 | zgrid = zeros([nx, ny], float)# z-coordinate |
---|
93 | zzgrid = zeros([nx, ny], float)# z-coordinate grouped |
---|
94 | z2grid = zeros([nx, ny], float)# z-coordinate ² |
---|
95 | zz2grid = zeros([nx, ny], float)# z-coordinate ² grouped |
---|
96 | |
---|
97 | for i in range(0, N): |
---|
98 | zgrid[iix[i], iiy[i]] = zgrid[iix[i],iiy[i]]+zz[i] |
---|
99 | z2grid[iix[i], iiy[i]] = z2grid[iix[i],iiy[i]]+zz[i]*zz[i] |
---|
100 | ngrid[iix[i], iiy[i]] = ngrid[iix[i],iiy[i]]+ 1 |
---|
101 | |
---|
102 | for i in range(0,ny): |
---|
103 | jj=round(3600/dxx[i]) |
---|
104 | c=reshape(ngrid[:,i],(jj,dxx[i])) |
---|
105 | b=c.sum(axis=1) |
---|
106 | d=reshape(zgrid[:,i],(jj,dxx[i])) |
---|
107 | e=d.sum(axis=1) |
---|
108 | f=reshape(z2grid[:,i],(jj,dxx[i])) |
---|
109 | g=f.sum(axis=1) |
---|
110 | for j in range(0, nx): |
---|
111 | ii=ceil((j+1)/dxx[i])-1 |
---|
112 | nngrid[j,i]=b[ii] |
---|
113 | zzgrid[j,i]=e[ii] |
---|
114 | zz2grid[j,i]=g[ii] |
---|
115 | |
---|
116 | zzpgrid = zeros([nx, ny], float) |
---|
117 | zzpgrid = zzgrid/nngrid |
---|
118 | sigma_grid=sqrt(zz2grid/nngrid-zzpgrid*zzpgrid) |
---|
119 | |
---|
120 | |
---|
121 | zgrid = np.transpose(zgrid) |
---|
122 | zzgrid = np.transpose(zzgrid) |
---|
123 | zzpgrid = np.transpose(zzpgrid) |
---|
124 | zz2grid = np.transpose(zz2grid) |
---|
125 | ngrid = np.transpose(ngrid) |
---|
126 | nngrid = np.transpose(nngrid) |
---|
127 | sigma_grid=np.transpose(sigma_grid) |
---|
128 | |
---|