source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/regrid_distrib.py @ 56

Last change on this file since 56 was 48, checked in by lahlod, 10 years ago

new scripts
File size: 6.5 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
10import draw_map
11
12####################################################
13# from polar (lon, lat) to cartesian coords (x, y) #
14####################################################
15
16def new_cartesian_grid_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 = np.zeros([ny, nx], float)
49    ngrid = np.zeros([ny, nx], float)
50    z2grid = np.zeros([ny, nx], float)
51    sigmagrid = np.zeros([ny, nx], float)
52    bblat = nonzero(lati >= 65.) # only select high latitude values of z
53    longitude = longi[bblat]
54    latitude = lati[bblat]
55    z = z_ini[bblat]
56    #################################################################################
57    # associates a (x, y) couple to each point of (longitude, latitude) coordinates #
58    #################################################################################
59    # origin of the (x, y) grid : (x=0, y=0) <=> (lon=0, lat=0)
60    L = len(z)
61    Rt = 6371.
62    alpha = (pi*Rt)/180.
63    beta = pi/180.
64    x = np.zeros([L], float)
65    y = np.zeros([L], float)
66    for k in range (0, L):
67            if ((longitude[k] >= 0.) & (longitude[k] <= 90.)): # 4eme quadrant
68                    theta = (90. - longitude[k]) * beta
69                    x[k] = (90. - latitude[k]) * alpha * cos(theta)
70                    y[k] = (90. - latitude[k]) * alpha * sin(-theta)
71            else:
72                    if ((longitude[k] > 90.) & (longitude[k] <= 180.)): # 1er quadrant
73                            theta = (longitude[k] - 90.) * beta
74                            x[k] = (90. - latitude[k]) * alpha * cos(theta)
75                            y[k] = (90. - latitude[k]) * alpha * sin(theta)
76                    else: 
77                            if ((longitude[k] >= -180.) & (longitude[k] < 0.)): # 2eme et 3eme quadrants
78                                    theta = (270. + longitude[k]) * beta
79                                    x[k] = (90. - latitude[k]) * alpha * cos(theta)
80                                    y[k] = (90. - latitude[k]) * alpha * sin(theta)
81    #################################################
82    # counting each x and y value in new grid cells #
83    #################################################
84    ix = np.zeros([L], int)
85    i = 0
86    for i in range (0, L): # boucle sur les points M (abscisses)
87            if x[i] == x0:
88                    ix[i] = 0
89            else:
90                    ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x vakue to a cell number
91    i = 0
92    iy = np.zeros([L], int) 
93    for i in range (0, L):# boucle sur les points M (ordonnees)
94            if y[i] == y0:
95                    iy[i] = 0
96            else:
97                    iy[i] = math.ceil((y[i] - y0) / dy) - 1 # associates each y vakue to a cell number
98    #########################################################################################
99    # calculation of distances between point M(x,y) and 4 grid points of its belonging cell #
100    #########################################################################################
101    close_point = np.zeros([L, 2], int)
102    k = 0
103    for k in range (0, L): # boucle sur les points M (x et y)
104            d1 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2))
105            d2 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k]]) ** 2))
106            d3 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2))
107            d4 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2))
108            d_vec = np.array([d1, d2, d3, d4])
109            ind = np.where(d_vec == min(d_vec)) # finds the smallest distance between the 4 points of the grid
110            point_vec = np.array([(ix[k], iy[k]), (ix[k] + 1, iy[k]), (ix[k], iy[k] + 1), (ix[k] + 1, iy[k] + 1)]) 
111            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)
112        ################################################
113        # association of z value to closest grid point #
114        ################################################
115    zgrid = np.zeros([ny, nx], float) # z in new grid
116    ngrid = np.zeros([ny, nx], int) # nb of obs per new grid cell
117    z2grid = np.zeros([ny, nx], float) # z**2 in new grid
118    for k in range (0, L):
119            zgrid[close_point[k, 1], close_point[k, 0]] = zgrid[close_point[k, 1], close_point[k, 0]] + z[k]
120            ngrid[close_point[k, 1], close_point[k, 0]] = ngrid[close_point[k, 1], close_point[k, 0]] + 1
121            z2grid[close_point[k, 1], close_point[k, 0]] = z2grid[close_point[k, 1], close_point[k, 0]] + (z[k] * z[k])
122    ZGRID = zgrid / ngrid
123    ##############################
124    # variance in each grid cell #
125    ##############################
126    sigmagrid = np.zeros([ny, nx], float)
127    for j in range (0, nx):
128            for i in range (0, ny):
129                    if (ngrid[i, j] > 1): # take away the cells where no or one obs
130                            sigmagrid[i, j] = sqrt((1./(ngrid[i, j]-1.))*(z2grid[i, j] - ngrid[i, j]*ZGRID[i, j]*ZGRID[i, j]))
131                    else:
132                            if (ngrid[i, j] == 1):
133                                    sigmagrid[i, j] = 0.
134                            else:
135                                    sigmagrid[i, j] = NaN
136
137
138
139    return ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart
140
141
142
143
144
Note: See TracBrowser for help on using the repository browser.