In [1]:
from netCDF4 import Dataset
url = ('http://geoport.whoi.edu/thredds/dodsC/clay/usgs/users/'
'jcwarner/Projects/Sandy/triple_nest/00_dir_NYB05.ncml')
#url = '00_dir_NYB05.nc'
nc = Dataset(url)
In [2]:
import pysgrid
sgrid = pysgrid.load_sgrid(nc)
sgrid # We need a better __repr__ and __str__ !!!
Out[2]:
In [3]:
sgrid.edge1_coordinates, sgrid.edge1_dimensions, sgrid.edge1_padding
Out[3]:
In [4]:
u_var = sgrid.u
u_var.center_axis, u_var.node_axis
Out[4]:
In [5]:
v_var = sgrid.v
v_var.center_axis, v_var.node_axis
Out[5]:
In [6]:
u_var.center_slicing
Out[6]:
In [7]:
v_var.center_slicing
Out[7]:
(Don't be scared, you do not need the sgrid object to get the variables. This just shows that there is a one-to-one mapping from the sgrid object to the netCDF4 object.)
In [8]:
u_velocity = nc.variables[u_var.variable]
v_velocity = nc.variables[v_var.variable]
In [9]:
from datetime import datetime, timedelta
from netCDF4 import date2index
t_var = nc.variables['ocean_time']
start = datetime(2012, 10, 30, 0, 0)
time_idx = date2index(start, t_var, select='nearest')
v_idx = 0
# Slice of the slice!
u_data = u_velocity[time_idx, v_idx, u_var.center_slicing[-2], u_var.center_slicing[-1]]
v_data = v_velocity[time_idx, v_idx, v_var.center_slicing[-2], v_var.center_slicing[-1]]
In [10]:
angle = sgrid.angle
angles = nc.variables[angle.variable][angle.center_slicing]
In [11]:
sgrid.topology_dimension
Out[11]:
In [12]:
from pysgrid.processing_2d import avg_to_cell_center
u_avg = avg_to_cell_center(u_data, u_var.center_axis)
v_avg = avg_to_cell_center(v_data, v_var.center_axis)
In [13]:
from pysgrid.processing_2d import rotate_vectors
u_rot, v_rot = rotate_vectors(u_avg, v_avg, angles)
In [14]:
from pysgrid.processing_2d import vector_sum
uv_vector_sum = vector_sum(u_rot, v_rot)
In [15]:
grid_cell_centers = sgrid.centers # Array of lon, lat pairs.
lon_var_name, lat_var_name = sgrid.face_coordinates
# use the longitude and latitude variable names to get their respective sgrid variables
sg_lon = getattr(sgrid, lon_var_name)
sg_lat = getattr(sgrid, lat_var_name)
lon_data = grid_cell_centers[..., 0][sg_lon.center_slicing]
lat_data = grid_cell_centers[..., 1][sg_lat.center_slicing]
In [16]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io import shapereader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
def make_map(projection=ccrs.PlateCarree(), figsize=(9, 9)):
fig, ax = plt.subplots(figsize=figsize,
subplot_kw=dict(projection=projection))
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
return fig, ax
In [17]:
sub = 10
scale = 0.06
fig, ax = make_map()
kw = dict(scale=1.0/scale, pivot='middle', width=0.003, color='black')
q = plt.quiver(lon_data[::sub, ::sub], lat_data[::sub, ::sub],
u_rot[::sub, ::sub], v_rot[::sub, ::sub], zorder=2, **kw)
cs = plt.pcolormesh(lon_data[::sub, ::sub],
lat_data[::sub, ::sub],
uv_vector_sum[::sub, ::sub], zorder=1, cmap=plt.cm.rainbow)
ax.coastlines('10m');