1 | #!/usr/bin/env python |
---|
2 | # -*- coding: utf-8 -*- |
---|
3 | import string |
---|
4 | import numpy as np |
---|
5 | import matplotlib.pyplot as plt |
---|
6 | import ffgrid2 |
---|
7 | from pylab import * |
---|
8 | from mpl_toolkits.basemap import Basemap |
---|
9 | from mpl_toolkits.basemap import shiftgrid, cm |
---|
10 | import draw_map |
---|
11 | |
---|
12 | #################################################### |
---|
13 | # from polar (lon, lat) to cartesian coords (x, y) # |
---|
14 | #################################################### |
---|
15 | |
---|
16 | |
---|
17 | |
---|
18 | # associates a (x, y) couple to each point of (longitude, latitude) coordinates |
---|
19 | # origin of the (x, y) grid : (x=0, y=0) <=> (lon=0, lat=0) |
---|
20 | L = len(z) |
---|
21 | #z0 = min(z) |
---|
22 | #z1 = max(z) |
---|
23 | Rt = 6371. |
---|
24 | alpha = (pi*Rt)/180. |
---|
25 | beta = pi/180. |
---|
26 | x = np.zeros([L], float) |
---|
27 | y = np.zeros([L], float) |
---|
28 | for k in range (0, L): |
---|
29 | # je lis longi[k], lati[k] |
---|
30 | if ((longi[k] >= 0.) & (longi[k] <= 90.)): # 4eme quadrant |
---|
31 | theta = (90. - longi[k]) * beta |
---|
32 | x[k] = (90. - lati[k]) * alpha * cos(theta) |
---|
33 | y[k] = (90. - lati[k]) * alpha * sin(-theta) |
---|
34 | else: |
---|
35 | if ((longi[k] > 90.) & (longi[k] <= 180.)): # 1er quadrant |
---|
36 | theta = (longi[k] - 90.) * beta |
---|
37 | x[k] = (90. - lati[k]) * alpha * cos(theta) |
---|
38 | y[k] = (90. - lati[k]) * alpha * sin(theta) |
---|
39 | else: |
---|
40 | if ((longi[k] >= -180.) & (longi[k] < 0.)): # 2eme et 3eme quadrants |
---|
41 | theta = (270. + longi[k]) * beta |
---|
42 | x[k] = (90. - lati[k]) * alpha * cos(theta) |
---|
43 | y[k] = (90. - lati[k]) * alpha * sin(theta) |
---|
44 | |
---|
45 | |
---|
46 | # definition of the new cartesian grid (x, y) # |
---|
47 | Mx = max(x) |
---|
48 | mx = min(x) |
---|
49 | x0 = round(mx) - 50. # xmin |
---|
50 | x1 = round(Mx) + 50. # xmax |
---|
51 | dx = 40. # delta(x) |
---|
52 | xvec = np.arange(x0, x1+dx, dx) |
---|
53 | nx = len(xvec) |
---|
54 | My = max(y) |
---|
55 | my = min(y) |
---|
56 | y0 = round(my) - 50. # ymin |
---|
57 | y1 = round(My) + 50. # ymax |
---|
58 | dy = 40. # delta(y) |
---|
59 | yvec = np.arange(y0, y1+dy, dy) |
---|
60 | ny = len(yvec) |
---|
61 | xgrid_cart, ygrid_cart= np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole) |
---|
62 | |
---|
63 | # counting each grid cell # |
---|
64 | ix = np.zeros([L], int) |
---|
65 | i = 0 |
---|
66 | for i in range (0, L): # boucle sur les points M (abscisses) |
---|
67 | if x[i] == x0: |
---|
68 | ix[i] = 0 |
---|
69 | else: |
---|
70 | ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x vakue to a cell number |
---|
71 | |
---|
72 | i = 0 |
---|
73 | iy = np.zeros([L], int) |
---|
74 | for i in range (0, L):# boucle sur les points M (ordonnees) |
---|
75 | if y[i] == y0: |
---|
76 | iy[i] = 0 |
---|
77 | else: |
---|
78 | iy[i] = math.ceil((y[i] - y0) / dy) - 1 # associates each y vakue to a cell number |
---|
79 | |
---|
80 | # calculation of distances between point M(x,y) and 4 grid points of its belonging cell # |
---|
81 | close_point = np.zeros([L, 2], int) |
---|
82 | k = 0 |
---|
83 | for k in range (0, L): # boucle sur les points M (x et y) |
---|
84 | # print 'x', x[k] |
---|
85 | # print 'y', y[k] |
---|
86 | # print 'xvec', xvec[ix[k]] |
---|
87 | # print 'yvec', yvec[iy[k]] |
---|
88 | d1 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2)) |
---|
89 | # print 'd1', d1 |
---|
90 | d2 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2)) |
---|
91 | # print 'd2', d2 |
---|
92 | d3 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2)) |
---|
93 | # print 'd3', d3 |
---|
94 | d4 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2)) |
---|
95 | # print 'd4', d4 |
---|
96 | d_vec = np.array([d1, d2, d3, d4]) |
---|
97 | ind = np.where(d_vec == min(d_vec)) # finds the smallest distance between the 4 points of the grid |
---|
98 | # print 'grid cell no : ' + str(ind[0][0]+1) |
---|
99 | point_vec = np.array([(ix[k], iy[k]), (ix[k] + 1, iy[k]), (ix[k], iy[k] + 1), (ix[k] + 1, iy[k] + 1)]) |
---|
100 | # print 'chosen point in grid', point_vec[ind[0][0]] # we have chosen which point of the grid is closer to M // point_vec[ind[0][0]] = (cell no of closest x, celle no of closest y) |
---|
101 | close_point[k, :] = point_vec[ind[0][0]]# we have chosen which point of the grid is closer to M // point_vec[ind[0][0]] = (cell no of closest x, celle no of closest y) |
---|
102 | |
---|
103 | # association of z value to closest grid point # |
---|
104 | zgrid = np.zeros([ny, nx], float) # z in new grid |
---|
105 | ngrid = np.zeros([ny, nx], int) # nb of obs per new grid cell |
---|
106 | z2grid = np.zeros([ny, nx], float) # z2 in new grid |
---|
107 | for k in range (0, L): |
---|
108 | zgrid[close_point[k, 1], close_point[k, 0]] = zgrid[close_point[k, 1], close_point[k, 0]] + z[k] |
---|
109 | ngrid[close_point[k, 1], close_point[k, 0]] = ngrid[close_point[k, 1], close_point[k, 0]] + 1 |
---|
110 | z2grid[close_point[k, 1], close_point[k, 0]] = z2grid[close_point[k, 1], close_point[k, 0]] + (z[k] * z[k]) |
---|
111 | |
---|
112 | # z in each point of new grid # |
---|
113 | ZGRID = zgrid / ngrid |
---|
114 | |
---|
115 | # variance of z in each grid cell # |
---|
116 | sigmagrid = np.zeros([ny, nx], float) |
---|
117 | for j in range (0, nx): |
---|
118 | for i in range (0, ny): |
---|
119 | if (ngrid[i, j] > 1): # take away the cells where no obs |
---|
120 | sigmagrid[i, j] = (1 / (ngrid[i, j] - 1) * sqrt(z2grid[i, j] - ngrid[i, j] * ZGRID[i, j] * ZGRID[i, j])) |
---|
121 | else: |
---|
122 | sigmagrid[i, j] = NaN |
---|
123 | |
---|
124 | |
---|
125 | |
---|
126 | |
---|
127 | figure() |
---|
128 | clevs = np.arange(0.6, 1.01, 0.01) |
---|
129 | cs=plt.contourf(xgrid_cart,ygrid_cart, ZGRID, clevs, cmap=cm.s3pcpn_l_r) |
---|
130 | colorbar(cs) |
---|
131 | |
---|
132 | |
---|
133 | |
---|
134 | |
---|
135 | #max((reshape(sigmagrid, size(sigmagrid)))[nonzero(isnan(reshape(sigmagrid, size(sigmagrid)) == False))] |
---|
136 | |
---|
137 | |
---|