source: trunk/src/scripts_Laura/ARCTIC/read_bathy.py @ 54

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


## read bathymetry data from netcdf file 'ETOPO1_Ice_g_gmt4.nc' ##
bathy = Dataset ('/net/argos/data/parvati/lahlod/ARCTIC/ETOPO1_Ice_g_gmt4.nc','r', format = 'NETCDF4')
#bathy.variables
lon_bath = bathy.variables['lon'][:]
lat_bath = bathy.variables['lat'][:]
z_bath = bathy.variables['z'][:]
bathy.close()


lat_limit = 50.


## read all points with bathymerty around 0, create z_cote 0 or 1 matrix ##
z_cote = np.zeros([len(lat_bath[np.where(lat_bath >= lat_limit)]), len(lon_bath)], float)
nn = min(np.where(lat_bath >= lat_limit)[0])
for ilon in range (0, len(lon_bath)):
    for ilat in range (0, len(lat_bath[np.where(lat_bath >= lat_limit)])):
        if ((z_bath[ilat+nn, ilon] == 0.)):
            z_cote[ilat, ilon] = 1.
        else: 
            z_cote[ilat, ilon] = NaN


## write z_cote in netcdf file 'z_cote_global_final.nc' ##
rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/z_cote_global_final_z0_lat50.nc', 'w', format = 'NETCDF4')
rootgrp.createDimension('longitude', len(lon_bath))
rootgrp.createDimension('latitude', len(lat_bath[np.where(lat_bath >= lat_limit)]))
nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))
nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))
nc_zcote = rootgrp.createVariable('cote', 'f', ('latitude', 'longitude',))
nc_lon[:] = lon_bath
nc_lat[:] = lat_bath[np.where(lat_bath >= lat_limit)]
nc_zcote[:] = z_cote
rootgrp.close()


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


## tests ##
plt.ion()
plt.figure()
m = Basemap(llcrnrlon=-180, urcrnrlon=180, llcrnrlat=-90., urcrnrlat=90, projection='cyl', resolution='c', fix_aspect=True)
m.drawcoastlines(linewidth = 1)
m.drawparallels(np.arange(-90., 90., 20.))
m.drawmeridians(np.arange(-180., 180., 20.))
xii,yii = m(*np.meshgrid(lon_cote, lat_cote))
cs = m.scatter(xii, yii, lim_cote)


## conversion to cartesian grid points ( lon, lat ) --> ( x[lon, lat], y[lon, lat] ) ##
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)


## write x and y equivalents of lon, lat points in netcdf file 'z_cote_global_final_cartes_z0.nc' ##
rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/z_cote_global_final_cartes_z0_lat50.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()


## read x and y equivalents of lon, lat points in netcdf file 'z_cote_global_final_cartes_z0.nc' ##
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.
plt.figure()
plt.scatter(x_ind, y_ind, c = z_ind, s = volume)