<jupyter_start><jupyter_text>Scenario: Velocity heat/vector map in Little Egg Inlet, NJ Get a dataset (In this case, ROMS) that covers the area<jupyter_code>from netCDF4 import Dataset
import numpy as np
import pysgrid

url = ('http://geoport-dev.whoi.edu/thredds/dodsC/clay/usgs/users/zdefne/run076/his/00_dir_roms_display.ncml')

nc = Dataset(url)
sgrid = pysgrid.load_grid(nc)
sgrid  # We need a better __repr__ and __str__ !!!<jupyter_output><empty_output><jupyter_text>Generate a mesh grid of points within the lat/lon bounds<jupyter_code>%%time
lons, lats = np.mgrid[-74.38:-74.26:600j, 39.45:39.56:600j]

points = np.stack((lons, lats), axis = -1)
                
print points.shape<jupyter_output>(600L, 600L, 2L)
Wall time: 12 ms<jupyter_text>After determining the time, depth, and points you are interested in, just interpolate!<jupyter_code>%%time
time_idx = 0
depth_idx = -1

interp_u = sgrid.interpolate_var_to_points(points, sgrid.u, slices=[time_idx, depth_idx])
interp_v = sgrid.interpolate_var_to_points(points, sgrid.v, slices=[time_idx, depth_idx])

#other function signatures:
interp_u2 = sgrid.interpolate_var_to_points(points, sgrid.u[time_idx, depth_idx])
interp_u3 = sgrid.interpolate_var_to_points(points, nc.variables['u'], slices=[time_idx, depth_idx])
print np.all(interp_u == interp_u2)
print np.all(interp_u == interp_u3)
#rotation is still ugly...
from pysgrid.processing_2d import rotate_vectors, vector_sum
from pysgrid.processing_2d import vector_sum

ind = sgrid.locate_faces(points)
ang_ind = ind + [1, 1]
angles = sgrid.angles[:][ang_ind[:, 0], ang_ind[:, 1]]
u_rot, v_rot = rotate_vectors(interp_u, interp_v, angles)
u_rot = u_rot.reshape(600,-1) # reshape for pcolormesh
v_rot = v_rot.reshape(600,-1)

uv_vector_sum = vector_sum(u_rot, v_rot)
%matplotlib inline

#Credit to Filipe for the plotting cells below

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=(12, 12)):
    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
mscale = 1
vscale = 15
scale = 0.04
lon_data = lons
lat_data = lats

fig, ax = make_map()

kw = dict(scale=1.0/scale, pivot='middle', width=0.003, color='black')
q = plt.quiver(lon_data[::vscale, ::vscale], lat_data[::vscale, ::vscale],
               u_rot[::vscale, ::vscale], v_rot[::vscale, ::vscale], zorder=2, **kw)

cs = plt.pcolormesh(lon_data[::mscale, ::mscale],
                    lat_data[::mscale, ::mscale],
                    uv_vector_sum[::mscale, ::mscale], zorder=1, cmap=plt.cm.rainbow)
ax.coastlines('10m');
plt.title('SGrid interpolation')
plt.show()<jupyter_output>C:\Users\Jay.Hennen\AppData\Local\Continuum\Anaconda2\lib\site-packages\matplotlib\artist.py:221: MatplotlibDeprecationWarning: This has been deprecated in mpl 1.5, please use the
axes property.  A removal date has not been set.
  warnings.warn(_get_axes_msg, mplDeprecation, stacklevel=1)
{'fraction_non_alphanumeric': 0.12294081757169005, 'fraction_numerical': 0.02593044539353264, 'mean_word_length': 5.072222222222222, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '398445', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1314, 'n_tokens_neox': 1245, 'n_words': 269, 'path': 'notebooks/SGrid_velocity_interpolation_demo.ipynb', 'repo_name': 'ocefpaf/pysgrid', 'size': '3278'}