source: trunk/src/scripts_Laura/ARCTIC/cartesian_grid_shoreline.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(lon, lat), y(lon, lat)) #
########################################################################

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


latitude = lat_cote
longitude = lon_cote

Rt = 6371.
alpha = (pi*Rt)/180.
beta = pi/180.
x = np.zeros([len(latitude), len(longitude)], float)
y = np.zeros([len(latitude), len(longitude)], float)

for ilon in range (0, len(longitude)):
    for ilat in range (0, len(latitude)):
        if ((longitude[ilon] >= 0.) & (longitude[ilon] <= 90.)): # 4eme quadrant
            theta = (90. - longitude[ilon]) * beta
            x[ilat, ilon] = (90. - latitude[ilat]) * alpha * cos(theta)
            y[ilat, ilon] = (90. - latitude[ilat]) * alpha * sin(-theta)
        else:
            if ((longitude[ilon] > 90.) & (longitude[ilon] <= 180.)): # 1er quadrant
                theta = (longitude[ilon] - 90.) * beta
                x[ilat, ilon] = (90. - latitude[ilat]) * alpha * cos(theta)
                y[ilat, ilon] = (90. - latitude[ilat]) * alpha * sin(theta)
            else: 
                if ((longitude[ilon] >= -180.) & (longitude[ilon] < 0.)): # 2eme et 3eme quadrants
                    theta = (270. + longitude[ilon]) * beta
                    x[ilat, ilon] = (90. - latitude[ilat]) * alpha * cos(theta)
                    y[ilat, ilon] = (90. - latitude[ilat]) * alpha * sin(theta)



rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/z_cote_global_final_cartes_z0.nc', 'w', format = 'NETCDF4')
rootgrp.createDimension('longitude', len(lon_cote))
rootgrp.createDimension('latitude', len(lat_cote))
nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))
nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))
nc_x = rootgrp.createVariable('x', 'f', ('latitude', 'longitude',))
nc_y = rootgrp.createVariable('y', 'f', ('latitude', 'longitude',))
nc_lon[:] = lon_cote
nc_lat[:] = lat_cote
nc_x[:] = x
nc_y[:] = y
rootgrp.close()



X = reshape(x, size(x))
Y = reshape(y, size(y))
Z_LIM = reshape(lim_cote, size(lim_cote))

indices = np.where(Z_LIM == 1.)[0]
x_ind = X[indices]
y_ind = Y[indices]
z_ind = Z_LIM[indices]


#x0 = x_ind.min() - 50.
#x1 = x_ind.max() + 50.
#dx = 10.
#xvec = np.arange(x0, x1+dx, dx)

#y0 = y_ind.min() - 50.
#y1 = y_ind.max() + 50.
#dy = 10.
#yvec = np.arange(y0, y1+dy, dy)

#xgrid, ygrid = np.meshgrid(x_vec, y_vec)


plt.figure()
plt.scatter(x_ind, y_ind, c = z_ind/10.)

# 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] = sqrt((1/(ngrid[i, j]-1))*(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