source: trunk/src/scripts_Laura/new_grid_bad.py @ 45

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import string
import numpy as np
import matplotlib.pyplot as plt
import ffgrid2
from pylab import *
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import shiftgrid, cm
import draw_map

####################################################
# from polar (lon, lat) to cartesian coords (x, y) #
####################################################



# associates a (x, y) couple to each point of (longitude, latitude) coordinates 
# origin of the (x, y) grid : (x=0, y=0) <=> (lon=0, lat=0) 
L = len(z)
#z0 = min(z)
#z1 = max(z)
Rt = 6371.
alpha = (pi*Rt)/180.
beta = pi/180.
x = np.zeros([L], float)
y = np.zeros([L], float)
for k in range (0, L):
   # je lis longi[k], lati[k]
    if ((longi[k] >= 0.) & (longi[k] <= 90.)): # 4eme quadrant
        theta = (90. - longi[k]) * beta
        x[k] = (90. - lati[k]) * alpha * cos(theta)
        y[k] = (90. - lati[k]) * alpha * sin(-theta)
    else:
        if ((longi[k] > 90.) & (longi[k] <= 180.)): # 1er quadrant
            theta = (longi[k] - 90.) * beta
            x[k] = (90. - lati[k]) * alpha * cos(theta)
            y[k] = (90. - lati[k]) * alpha * sin(theta)
        else: 
            if ((longi[k] >= -180.) & (longi[k] < 0.)): # 2eme et 3eme quadrants
                theta = (270. + longi[k]) * beta
                x[k] = (90. - lati[k]) * alpha * cos(theta)
                y[k] = (90. - lati[k]) * alpha * sin(theta)


# definition of the new cartesian grid (x, y) #
Mx = max(x)
mx = min(x)
x0 = round(mx) - 50. # xmin
x1 = round(Mx) + 50. # xmax
dx = 40. # delta(x)
xvec = np.arange(x0, x1+dx, dx)
nx = len(xvec)
My = max(y)
my = min(y)
y0 = round(my) - 50. # ymin
y1 = round(My) + 50. # ymax
dy = 40. # delta(y)
yvec = np.arange(y0, y1+dy, dy)
ny = len(yvec)
xgrid_cart, ygrid_cart= np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole)  

# counting each grid cell #
ix = np.zeros([L], int)
i = 0
for i in range (0, L): # boucle sur les points M (abscisses) 
    if x[i] == x0:
        ix[i] = 0
    else:
        ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x vakue to a cell number

i = 0
iy = np.zeros([L], int) 
for i in range (0, L):# boucle sur les points M (ordonnees) 
    if y[i] == y0:
        iy[i] = 0
    else:
        iy[i] = math.ceil((y[i] - y0) / dy) - 1 # associates each y vakue to a cell number

# calculation of distances between point M(x,y) and 4 grid points of its belonging cell #
close_point = np.zeros([L, 2], int)
k = 0
for k in range (0, L): # boucle sur les points M (x et y)
    #    print 'x', x[k]
    #    print 'y', y[k]
    #    print 'xvec', xvec[ix[k]]
    #    print 'yvec', yvec[iy[k]]
    d1 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2))
     #   print 'd1', d1
    d2 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2))
    #    print 'd2', d2
    d3 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2))
     #   print 'd3', d3
    d4 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2))
    #    print 'd4', d4
    d_vec = np.array([d1, d2, d3, d4])
    ind = np.where(d_vec == min(d_vec)) # finds the smallest distance between the 4 points of the grid
    #    print 'grid cell no : ' + str(ind[0][0]+1)
    point_vec = np.array([(ix[k], iy[k]), (ix[k] + 1, iy[k]), (ix[k], iy[k] + 1), (ix[k] + 1, iy[k] + 1)])
    #    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)
    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)

# association of z value to closest grid point #
zgrid = np.zeros([ny, nx], float) # z in new grid
ngrid = np.zeros([ny, nx], int) # nb of obs per new grid cell
z2grid = np.zeros([ny, nx], float) # z2 in new grid
for k in range (0, L):
    zgrid[close_point[k, 1], close_point[k, 0]] = zgrid[close_point[k, 1], close_point[k, 0]] + z[k]
    ngrid[close_point[k, 1], close_point[k, 0]] = ngrid[close_point[k, 1], close_point[k, 0]] + 1
    z2grid[close_point[k, 1], close_point[k, 0]] = z2grid[close_point[k, 1], close_point[k, 0]] + (z[k] * z[k])
    
# z in each point of new grid #
ZGRID = zgrid / ngrid

# variance of z in each grid cell #
sigmagrid = np.zeros([ny, nx], float)
for j in range (0, nx):
    for i in range (0, ny):
        if (ngrid[i, j] > 1): # take away the cells where no obs 
            sigmagrid[i, j] = (1 / (ngrid[i, j] - 1) * sqrt(z2grid[i, j] - ngrid[i, j] * ZGRID[i, j] * ZGRID[i, j]))
        else:
            sigmagrid[i, j] = NaN
        



figure()
clevs = np.arange(0.6, 1.01, 0.01)
cs=plt.contourf(xgrid_cart,ygrid_cart, ZGRID, clevs, cmap=cm.s3pcpn_l_r)
colorbar(cs)




#max((reshape(sigmagrid, size(sigmagrid)))[nonzero(isnan(reshape(sigmagrid, size(sigmagrid)) == False))]