source: trunk/src/scripts_Laura/ARCTIC/arctic_map.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
from netCDF4 import Dataset



def arctic_map_lat50():


    ## read z_cote from netcfd file 'z_cote_global_final.nc' ##
    cote_bin = Dataset ('/net/argos/data/parvati/lahlod/ARCTIC/z_cote_global_final_z0_lat50.nc','r', format = 'NETCDF4')
    #cote_bin.variables
    #lon_cote = cote_bin.variables['longitude'][:]
    #lat_cote = cote_bin.variables['latitude'][:]
    lim_cote = cote_bin.variables['cote'][:]
    cote_bin.close()

    coast_bin = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/z_cote_global_final_cartes_z0_lat50.nc', 'r', format = 'NETCDF4')
    #coast_bin.variables
    #lon_cote = coast_bin.variables['longitude'][:]
    #lat_cote = coast_bin.variables['latitude'][:]
    x_cote = coast_bin.variables['x'][:]
    y_cote = coast_bin.variables['y'][:]
    coast_bin.close()

    X = reshape(x_cote, size(x_cote))
    Y = reshape(y_cote, size(y_cote))
    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]

    volume = z_ind/50.
    #fig = plt.figure()
    #fig = plt.scatter(x_ind, y_ind, c = z_ind, s = volume)

    return(x_ind, y_ind, z_ind, volume)