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

Last change on this file since 55 was 46, checked in by lahlod, 10 years ago

new scripts
File size: 7.2 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3import string
4import numpy as np
5import matplotlib.pyplot as plt
6import ffgrid2
7from pylab import *
8from mpl_toolkits.basemap import Basemap
9from mpl_toolkits.basemap import shiftgrid, cm
10
11
12####################################################
13# from polar (lon, lat) to cartesian coords (x, y) #
14####################################################
15
16def new_cartesian_grid(nb_days, jour, month, longi, lati, z_ini, z0, z1, dx, dy):
17
18
19    ############################
20    # definition of input data #
21    ############################
22    # jour = vector of days (1D-array)
23    # month = number or days in month (integer)
24    # longi = longitude of data points (1D-array)
25    # lati = latitude of data points (1D-array)
26    # z_ini = data to re-grid (1D-array)
27    # z0 = min value of data (float)
28    # z1 = max value of data (float)
29
30
31    ###############################################
32    # definition of the new cartesian grid (x, y) #
33    ###############################################
34    x0 = -3000. # min limit of grid
35    x1 = 2500. # max limit of grid
36    xvec = np.arange(x0, x1+dx, dx)
37    nx = len(xvec) 
38    y0 = -3000. # min limit of grid
39    y1 = 3000. # max limit of grid
40    yvec = np.arange(y0, y1+dy, dy)
41    ny = len(yvec)
42    xgrid_cart, ygrid_cart= np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole)
43
44
45    ###################################
46    # calculation of the new gridding #
47    ###################################
48    zgrid_output = np.zeros([ny, nx, nb_days], float)
49    ngrid_output = np.zeros([ny, nx, nb_days], float)
50    z2grid_output = np.zeros([ny, nx, nb_days], float)
51    sigmagrid_output = np.zeros([ny, nx, nb_days], float)
52    bblat = nonzero(lati >= 65.) # only select high latitude values of z
53    for ijr in range (0, nb_days): # loop on time - here time = days
54        bbjour = nonzero(jour[bblat] == ijr + 1.)[0]
55        longitude = longi[bblat][bbjour]
56        latitude = lati[bblat][bbjour]
57        z = z_ini[bblat][bbjour]   
58        #################################################################################
59        # associates a (x, y) couple to each point of (longitude, latitude) coordinates #
60        #################################################################################
61        # origin of the (x, y) grid : (x=0, y=0) <=> (lon=0, lat=0)
62        L = len(z)
63        Rt = 6371.
64        alpha = (pi*Rt)/180.
65        beta = pi/180.
66        x = np.zeros([L], float)
67        y = np.zeros([L], float)
68        for k in range (0, L):
69                if ((longitude[k] >= 0.) & (longitude[k] <= 90.)): # 4eme quadrant
70                        theta = (90. - longitude[k]) * beta
71                        x[k] = (90. - latitude[k]) * alpha * cos(theta)
72                        y[k] = (90. - latitude[k]) * alpha * sin(-theta)
73                else:
74                        if ((longitude[k] > 90.) & (longitude[k] <= 180.)): # 1er quadrant
75                                theta = (longitude[k] - 90.) * beta
76                                x[k] = (90. - latitude[k]) * alpha * cos(theta)
77                                y[k] = (90. - latitude[k]) * alpha * sin(theta)
78                        else: 
79                                if ((longitude[k] >= -180.) & (longitude[k] < 0.)): # 2eme et 3eme quadrants
80                                        theta = (270. + longitude[k]) * beta
81                                        x[k] = (90. - latitude[k]) * alpha * cos(theta)
82                                        y[k] = (90. - latitude[k]) * alpha * sin(theta)
83        #################################################
84        # counting each x and y value in new grid cells #
85        #################################################
86        ix = np.zeros([L], int)
87        i = 0
88        for i in range (0, L): # boucle sur les points M (abscisses)
89                if x[i] == x0:
90                        ix[i] = 0
91                else:
92                        ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x vakue to a cell number
93        i = 0
94        iy = np.zeros([L], int) 
95        for i in range (0, L):# boucle sur les points M (ordonnees)
96                if y[i] == y0:
97                        iy[i] = 0
98                else:
99                        iy[i] = math.ceil((y[i] - y0) / dy) - 1 # associates each y vakue to a cell number
100        #########################################################################################
101        # calculation of distances between point M(x,y) and 4 grid points of its belonging cell #
102        #########################################################################################
103        close_point = np.zeros([L, 2], int)
104        k = 0
105        for k in range (0, L): # boucle sur les points M (x et y)
106                d1 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2))
107                d2 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2))
108                d3 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2))
109                d4 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2))
110                d_vec = np.array([d1, d2, d3, d4])
111                ind = np.where(d_vec == min(d_vec)) # finds the smallest distance between the 4 points of the grid
112                point_vec = np.array([(ix[k], iy[k]), (ix[k] + 1, iy[k]), (ix[k], iy[k] + 1), (ix[k] + 1, iy[k] + 1)]) 
113                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)
114        ################################################
115        # association of z value to closest grid point #
116        ################################################
117        zgrid = np.zeros([ny, nx], float) # z in new grid
118        ngrid = np.zeros([ny, nx], int) # nb of obs per new grid cell
119        z2grid = np.zeros([ny, nx], float) # z**2 in new grid
120        for k in range (0, L):
121                zgrid[close_point[k, 1], close_point[k, 0]] = zgrid[close_point[k, 1], close_point[k, 0]] + z[k]
122                ngrid[close_point[k, 1], close_point[k, 0]] = ngrid[close_point[k, 1], close_point[k, 0]] + 1
123                z2grid[close_point[k, 1], close_point[k, 0]] = z2grid[close_point[k, 1], close_point[k, 0]] + (z[k] * z[k])
124        ZGRID = zgrid / ngrid
125        ##############################
126        # variance in each grid cell #
127        ##############################
128        sigmagrid = np.zeros([ny, nx], float)
129        for j in range (0, nx):
130                for i in range (0, ny):
131                        if (ngrid[i, j] > 1): # take away the cells where no or one obs
132                                sigmagrid[i, j] = sqrt((1./(ngrid[i, j]-1.))*(z2grid[i, j] - ngrid[i, j]*ZGRID[i, j]*ZGRID[i, j]))
133                        else:
134                                if (ngrid[i, j] == 1):
135                                        sigmagrid[i, j] = 0.
136                                else:
137                                        sigmagrid[i, j] = NaN
138
139        zgrid_output[:, :, ijr] = ZGRID[:, :]
140        ngrid_output[:, :, ijr] = ngrid[:, :]
141        z2grid_output[:, :, ijr] = z2grid[:, :]
142        sigmagrid_output[:, :, ijr] = sigmagrid[:, :]
143
144
145    return zgrid_output, ngrid_output, z2grid_output, sigmagrid_output, xvec, yvec, xgrid_cart, ygrid_cart
146
147
148
149
150
Note: See TracBrowser for help on using the repository browser.