source: trunk/src/scripts_Laura/cartesian_grid_test.py @ 55

#!/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) #
####################################################

def new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy):

    # 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):
        if ((longitude[k] >= 0.) & (longitude[k] <= 90.)): # 4eme quadrant
            theta = (90. - longitude[k]) * beta
            x[k] = (90. - latitude[k]) * alpha * cos(theta)
            y[k] = (90. - latitude[k]) * alpha * sin(-theta)
        else:
            if ((longitude[k] > 90.) & (longitude[k] <= 180.)): # 1er quadrant
                theta = (longitude[k] - 90.) * beta
                x[k] = (90. - latitude[k]) * alpha * cos(theta)
                y[k] = (90. - latitude[k]) * alpha * sin(theta)
            else: 
                if ((longitude[k] >= -180.) & (longitude[k] < 0.)): # 2eme et 3eme quadrants
                    theta = (270. + longitude[k]) * beta
                    x[k] = (90. - latitude[k]) * alpha * cos(theta)
                    y[k] = (90. - latitude[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 = 50. # 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 = 50. # 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
        

    return ZGRID, ngrid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart 







#emisGRID, nemis, sigmaemis, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(lon_zon, lat_zon, emis_zon, emis_zon.min(), emis_zon.max(), 40, 40)