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

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



###############
# test values #
###############
longi = array([-109.208, -117.39 , -117.928, -123.927, -120.325, -128.82 , -125.947, -122.596, -133.173, -124.778, -137.312, -134.984, -132.538, -129.883, -126.898, -136.826, -134.414, -131.833, -128.978, -141.05 , -138.707, -136.313, -133.787, -131.036, -127.931, -140.634, -138.243, -135.757, -133.087, -130.122, -126.704, -110.352, -142.613, -140.213, -137.751, -135.144, -132.291, -129.058, -125.246, -114.38 , -147.118, -144.65 , -142.229, -139.777, -137.216, -134.453, -131.372, -127.803, -118.005, -146.753])

lati = array([ 68.164,  69.205,  70.121,  70.087,  70.422,  69.994,  70.359, 70.709,  69.846,  70.983,  69.619,  70.061,  70.474,  70.867, 71.246,  70.259,  70.693,  71.103,  71.496,  69.936,  70.437, 70.896,  71.325,  71.735,  72.129,  70.593,  71.08 ,  71.532, 71.959,  72.369,  72.764,  73.708,  70.723,  71.244,  71.721, 72.168,  72.595,  73.007,  73.403,  74.075,  70.188,  70.824, 71.384,  71.891,  72.361,  72.807,  73.235,  73.648,  74.394, 70.892])

z = array([ 0.938 ,  0.899 ,  0.958 ,  0.944 ,  0.991 ,  0.867 ,  0.948 , 0.963 ,  0.911 ,  0.978 ,  0.935 ,  0.932 ,  0.957 ,  0.972 , 0.983 ,  0.945 ,  0.941 ,  0.965 ,  0.986 ,  0.88  ,  0.942 , 0.945 ,  0.966 ,  0.989 ,  0.9999,  0.914 ,  0.908 ,  0.93  , 0.96  ,  0.996 ,  0.9999,  0.952 ,  0.898 ,  0.88  ,  0.852 , 0.86  ,  0.906 ,  0.984 ,  0.976 ,  0.888 ,  0.794 ,  0.89  , 0.846 ,  0.849 ,  0.869 ,  0.812 ,  0.764 ,  0.904 ,  0.853 , 0.934 ])
z0 = min(z)
z1 = max(z)


####################################################
# 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) 
Rt = 6371.
alpha = (pi*Rt)/180.
beta = pi/180.
# definition of theta angle and coord x,y #
x = np.zeros([len(lati), len(longi)], float)
y = np.zeros([len(lati), len(longi)], float)
for ilat in range (0, len(lati)):
    for ilon in range (0, len(longi)):
        if ((longi[ilon] >= 0.) & (longi[ilon] <= 90.)): # 4eme quadrant
            theta = (90. - longi[ilon]) * beta
            x[ilat, ilon] = (90. - lati[ilat]) * alpha * cos(theta)
            y[ilat, ilon] = (90. - lati[ilat]) * alpha * sin(-theta)
        else:
            if ((longi[ilon] > 90.) & (longi[ilon] <= 180.)): # 1er quadrant
                theta = (longi[ilon] - 90.) * beta
                x[ilat, ilon] = (90. - lati[ilat]) * alpha * cos(theta)
                y[ilat, ilon] = (90. - lati[ilat]) * alpha * sin(theta)
            else: 
                if ((longi[ilon] >= -180.) & (longi[ilon] < 0.)): # 2eme et 3eme quadrants
                    theta = (270. + longi[ilon]) * beta
                    x[ilat, ilon] = (90. - lati[ilat]) * alpha * cos(theta)
                    y[ilat, ilon] = (90. - lati[ilat]) * alpha * sin(theta)


# definition of the new cartesian grid (x, y) #
Sx = size(x)
Sy = size(y)
xx = reshape(x, size(x))
yy = reshape(y, size(y))
mx = min(xx)
Mx = max(xx)
my = min(yy)
My = max(yy)

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


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

iy = np.zeros([Sy], int) 
for j in range (0, Sy): # boucle sur les points M (ordonnees)
    if yy[j] == y0:
        iy[j] = 0
    else:
        iy[j] = math.ceil((yy[j] - y0) / dy) - 1 


#inx = (ix >= 0) & (ix < nx)
#iny = (iy >= 0) & (iy < ny)
#inz = (z > z0) & (z < z1)
#inn = inx & iny & inz
#iix = ix[inn]
#iiy = iy[inn]
#zz = z[inn]




# calculation of distances between point M(x,y) and 4 grid points of its belonging cell #
dist1 = np.zeros([Sx], float)
dist2 = np.zeros([Sx], float)
dist3 = np.zeros([Sx], float)
dist4 = np.zeros([Sx], float)
min_dist = np.zeros([Sx], float)
for k in range (0, Sx): # boucle sur les points M (x et y)
    dist1[k] = sqrt(((xx[k] - xvec[ix[k]]) ** 2) + ((yy[k] - yvec[iy[k]]) ** 2))
    dist2[k] = sqrt(((xx[k] - xvec[ix[k] + 1]) ** 2) + ((yy[k] - yvec[iy[k]]) ** 2))
    dist3[k] = sqrt(((xx[k] - xvec[ix[k]])**2) + ((yy[k] - yvec[iy[k] + 1]) ** 2))
    dist4[k] = sqrt(((xx[k] - xvec[ix[k] + 1]) ** 2) + ((yy[k] - yvec[iy[k] + 1]) ** 2))
    dist_vec = np.array([dist1[k], dist2[k], dist3[k], dist4[k]])
    ind = np.where(dist_vec == min(dist_vec))
    print 'grid cell no : ' + str(ind[0][0]+1)
    

##### probleme: z pas la même dimension (plus petite)!!#######






















for i in range (0, 5):
    np.where((xx - xvec[i]) <= dx)
        
    for j in range (0, len(yy)):
        np.where ((xx[i] - xvec) <= dx):
        print 'ok'
        else:
            print 'not ok'
        np.where((yy[i] - yvec) <= dy)



M = np.array([xx, yy])
for k in range (0, len(xx)):
    xcoord[k] = M[:,k][0]
    ycoord[k] = M[:,k][1]




ix = np.zeros([len(xx)], int)
for i in range (0, len(xx)):
    if xx[i] == x0:
        ix[i] = 0
    else:
        ix[i] = math.ceil((xx[i] - x0) / dx) - 1

iy = np.zeros([len(yy)], int)
for j in range (0, len(yy)):
    if yy[j] == y0:
        iy[j] = 0
    else:
        iy[j] = math.ceil((yy[j] - y0) / dy) - 1 







figure()
m = Basemap(projection = 'npaeqd',boundinglat = 60, lon_0 = 0, resolution = 'l')
clevs=arange(0.5, 1., 0.01)
m.drawcoastlines(linewidth=1)
m.drawparallels(np.arange(60, 90, 20))
m.drawmeridians(np.arange(-180, 180, 20))
xii,yii = m(*np.meshgrid(xx ,yy))
cs=m.contourf(xii,yii, zgrid, clevs, cmap=cm.s3pcpn_l_r)
cbar =colorbar(cs)