[6] | 1 | #!/usr/bin/env python |
---|
| 2 | # -*- coding: utf-8 -*- |
---|
| 3 | import numpy as np |
---|
| 4 | from numpy import * |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | def newgridns(xx, yy, zz, m, vnx, nscan, pos, z0, z1): |
---|
| 8 | |
---|
| 9 | x=xx |
---|
| 10 | y=yy |
---|
| 11 | z=zz |
---|
| 12 | dx=0.1 |
---|
[7] | 13 | dy=1 |
---|
[6] | 14 | dparx=0.05 |
---|
| 15 | dpary=0.5 |
---|
| 16 | n0=0 |
---|
| 17 | y11= -50 |
---|
| 18 | x0, x1 = -180, 180 |
---|
| 19 | y0, y1 = -90, 90 |
---|
| 20 | dx=0.1 |
---|
| 21 | dy=1.0 |
---|
[7] | 22 | n=size(x) |
---|
[6] | 23 | xvec = arange(x0,x1,dx) |
---|
| 24 | yvec = arange(y0,y11,dy) |
---|
| 25 | nx = size(xvec) |
---|
| 26 | ny = size(yvec) |
---|
[7] | 27 | ix = zeros(n, float) |
---|
| 28 | iy = zeros(n, float) |
---|
[6] | 29 | |
---|
| 30 | for kk in range(0,n): |
---|
| 31 | if x[kk] == -180: |
---|
[9] | 32 | ix[kk] = 3600 |
---|
[6] | 33 | else: |
---|
[9] | 34 | ix[kk] = math.ceil((x[kk] - x0)/dx)-1 |
---|
[6] | 35 | |
---|
| 36 | for kk in range(0,n): |
---|
| 37 | if y[kk] == -90: |
---|
[9] | 38 | iy[kk] = 0 |
---|
[6] | 39 | else: |
---|
[9] | 40 | iy[kk] = math.ceil((y[kk] - y0)/dy)-1 |
---|
[6] | 41 | |
---|
[7] | 42 | |
---|
[6] | 43 | if nscan==1: |
---|
| 44 | ipos = ((pos <=5) | (pos >=26)) |
---|
[7] | 45 | if nscan==2: |
---|
[6] | 46 | ipos = ((pos >=6) & (pos <=10))| ((pos >=21) & (pos<=25)) |
---|
| 47 | if nscan==3: |
---|
| 48 | ipos = ((pos >=11) & (pos<=20)) |
---|
[7] | 49 | |
---|
| 50 | inx = (ix >= 0) & (ix < nx) |
---|
[6] | 51 | iny = (iy >= 0) & (iy < ny) |
---|
[7] | 52 | inz = (z >= z0) & (z <= z1) |
---|
[6] | 53 | inn = inx & iny & inz & ipos |
---|
| 54 | iix = ix[inn] |
---|
| 55 | iiy = iy[inn] |
---|
| 56 | zz = z[inn] |
---|
| 57 | mm = m[inn] |
---|
| 58 | N = size(iix)# how many datapoints are left now? |
---|
| 59 | #calculation of the each arc legth |
---|
| 60 | |
---|
| 61 | dxx=vnx |
---|
| 62 | ngrid = zeros([nx, ny], float)# no of obs per grid cell |
---|
| 63 | nngrid = zeros([nx, ny], float) |
---|
| 64 | zgrid = zeros([nx, ny], float)# z-coordinate |
---|
| 65 | zzgrid = zeros([nx, ny], float)# z-coordinate grouped |
---|
| 66 | zdgrid = zeros([nx, ny], float)# z-coordinate 2 |
---|
| 67 | zz2grid = zeros([nx, ny], float)# z-coordinate 2 grouped |
---|
| 68 | |
---|
| 69 | for i in range(0, N): |
---|
| 70 | zgrid[iix[i], iiy[i]] = zgrid[iix[i],iiy[i]]+zz[i] |
---|
| 71 | zdgrid[iix[i], iiy[i]] = zdgrid[iix[i],iiy[i]]+zz[i]*zz[i] |
---|
[7] | 72 | ngrid[iix[i], iiy[i]] = ngrid[iix[i],iiy[i]]+ 1 |
---|
| 73 | |
---|
[6] | 74 | for i in range(0,ny): |
---|
| 75 | jj=round(3600/dxx[i]) |
---|
| 76 | c=reshape(ngrid[:,i],(jj,dxx[i])) |
---|
| 77 | b=c.sum(axis=1) |
---|
| 78 | d=reshape(zgrid[:,i],(jj,dxx[i])) |
---|
| 79 | e=d.sum(axis=1) |
---|
| 80 | f=reshape(zdgrid[:,i],(jj,dxx[i])) |
---|
[7] | 81 | g=f.sum(axis=1) |
---|
[6] | 82 | for j in range(0, nx): |
---|
| 83 | ii=ceil((j+1)/dxx[i])-1 |
---|
| 84 | nngrid[j,i]=b[ii] |
---|
| 85 | zzgrid[j,i]=e[ii] |
---|
| 86 | zz2grid[j,i]=g[ii] |
---|
| 87 | |
---|
[7] | 88 | zzpgrid = zeros([nx, ny], float) |
---|
[6] | 89 | zzpgrid = zzgrid/nngrid |
---|
| 90 | sigma_grid=sqrt(zz2grid/nngrid-zzpgrid*zzpgrid) |
---|
| 91 | zgrid = np.transpose(zgrid) |
---|
| 92 | sigma_grid=np.transpose(sigma_grid) |
---|
| 93 | zzpgrid = np.transpose(zzpgrid) |
---|
| 94 | ngrid = np.transpose(ngrid) |
---|
| 95 | nngrid = np.transpose(nngrid) |
---|
| 96 | return zgrid, zzpgrid, ngrid, nngrid, sigma_grid, xvec, yvec, zz, mm |
---|