Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
notebook/folium/tri_tide_movie.ipynb
{
"cells": [
{
"metadata": {},
"cell_type": "markdown",
"source": "# Plotting predicted tides\n\njupyter notebook --NotebookApp.iopub_data_rate_limit=10000000 tri_tide_movie-signell.ipynb"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import netCDF4\n\n\n# Versions >1.2.4 are failing with http://nbviewer.jupyter.org/gist/anonymous/410bf777e19ac8ca14bb06a24b03745f\nprint(netCDF4.__version__)",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "from netCDF4 import Dataset\n\n\nurl = 'http://geoport.whoi.edu/thredds/dodsC/usgs/vault0/models/tides/vdatum_gulf_of_maine/adcirc54_38_orig.nc'\n\nnc = Dataset(url)\n\nprint(nc.variables.keys())",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "units = {\n 'knots': 1.9438,\n 'm/s': 1.0\n}\nconsts = ['STEADY', 'M2', 'S2', 'N2', 'K1', 'O1', 'P1', 'M4', 'M6']\n\n# Vineyard sound 2.\nbbox = [-70.7234, -70.4532, 41.4258, 41.5643]\n\nhalo = 0.1\n\nax2 = [\n bbox[0] - halo * (bbox[1] - bbox[0]),\n bbox[1] + halo * (bbox[1] - bbox[0]),\n bbox[2] - halo * (bbox[3] - bbox[2]),\n bbox[3] + halo * (bbox[3] - bbox[2])\n]",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "import pytz\nfrom datetime import datetime\nfrom pandas import date_range\n\n\nfmt = '%d-%b-%Y %H:%M'\ntzinfo = pytz.utc\n\nstart = datetime.strptime('18-Sep-2015 05:00', fmt).replace(tzinfo=tzinfo)\nstop = datetime.strptime('19-Sep-2015 05:00', fmt).replace(tzinfo=tzinfo)\n\nglocals = date_range(start, stop, freq='1H').to_pydatetime()\n\nntimes = len(glocals)",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "def parse_string(name):\n try:\n const = ''.join(name).strip()\n except TypeError:\n const = b''.join(name).strip()\n const = const.decode()\n return const ",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "data = nc['tidenames'][:]\n\nnames = []\nfor name in data.tolist():\n names.append(parse_string(name))",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "from scipy.spatial import Delaunay\n\n\ndepth = nc['depth'][:]\nlatf = nc['lat'][:]\nlonf = nc['lon'][:]\nfrequency = nc['tidefreqs'][:]\n\n# Not sure why this is not working.\n# trif = cubes.extract_strict('Horizontal Element Incidence List').data\ntrif = Delaunay(list(zip(lonf, latf))).vertices",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "# Find indices in box.\nimport numpy as np\n\n\ninbox = np.logical_and(\n np.logical_and(lonf >= ax2[0], lonf <= ax2[1]),\n np.logical_and(latf >= ax2[2], latf <= ax2[3])\n)\n\nlon = lonf[inbox]\nlat = latf[inbox]",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "from scipy.io import loadmat\n\n\nmat = 't_constituents.mat'\ncon_info = loadmat(mat, squeeze_me=True)\n# ignoring shallow water and sat constants.\ncon_info = con_info['const']",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"scrolled": true,
"trusted": true
},
"cell_type": "code",
"source": "from utide import _ut_constants_fname, loadbunch\n\n\ncon_info = loadbunch(_ut_constants_fname)['const']",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "# Find the indices of the tidal constituents.\n\nk = 0\nind_nc, ind_ttide = [], []\n\nconst_name = [e.strip() for e in con_info['name'].tolist()]\n\nfor name in consts:\n try:\n if name == 'STEADY':\n indx = const_name.index('Z0')\n else:\n indx = const_name.index(name)\n k += 1\n ind_ttide.append(indx)\n ind_nc.append(names.index(name))\n except ValueError:\n pass # `const` not found.",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "uamp = nc['u_amp'][0, inbox, :][:, ind_nc]\nvamp = nc['v_amp'][0, inbox, :][:, ind_nc]\n\nupha = nc['u_phase'][0, inbox, :][:, ind_nc]\nvpha = nc['v_phase'][0, inbox, :][:, ind_nc]",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "freq_nc = frequency[ind_nc]\n\nfreq_ttide = con_info['freq'][ind_ttide]\n\nt_tide_names = np.array(const_name)[ind_ttide]",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "omega_ttide = 2*np.pi * freq_ttide # Convert from radians/s to radians/hour.\n\nomega = freq_nc * 3600\n\nrllat = 55 # Reference latitude for 3rd order satellites (degrees) (55 is fine always)\n\nvscale = 10 # smaller makes arrows bigger",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "from matplotlib.dates import date2num\n\n\n# Convert to Matlab datenum.\n# (Soon UTide will take python datetime objects.)\njd_start = date2num(start) + 366.1667",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "from utide.harmonics import FUV\n\n\n# NB: I am not a 100% sure if this is identical to what we had with t_tide.\n# ngflgs -> [NodsatLint NodsatNone GwchLint GwchNone]\nv, u, f = FUV(\n t=np.array([jd_start]),\n tref=np.array([0]),\n lind=np.array([ind_ttide]),\n lat=rllat,\n ngflgs=[0, 0, 0, 0]\n)",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "# Convert phase in radians.\nv, u, f = map(np.squeeze, (v, u, f))\nv = v * 2 * np.pi\nu = u * 2 * np.pi\n\nthours = np.array(\n [d.total_seconds() for d in (glocals - glocals[0])]\n) / 60 / 60.",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true,
"collapsed": true
},
"cell_type": "code",
"source": "import json\nimport mplleaflet\nimport matplotlib.pyplot as plt\n\n\nfig = plt.figure(figsize=(10,5))\nU = (f * uamp * np.cos(v + thours[k] * omega + u - upha * np.pi/180)).sum(axis=1)\nV = (f * vamp * np.cos(v + thours[k] * omega + u - vpha * np.pi/180)).sum(axis=1)\n\nw = units['knots'] * (U + 1j * V)\n\nwf = np.NaN * np.ones_like(lonf, dtype=w.dtype)\nwf[inbox] = w\n\n# FIXME: Cannot use masked arrays and tricontour!\n# wf = ma.masked_invalid(wf)\n# cs = ax.tricontour(lonf, latf, trif, np.abs(wf).filled(fill_value=0))\n# fig.colorbar(cs)\nq = plt.quiver(lon, lat, U, V, scale=vscale, color='white')\nplt.axis(bbox) # Vineyard sound 2.\n#q.set_title('{}'.format(glocals[k]))\nmplleaflet.display(fig=fig, tiles='esri_aerial')",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "js = mplleaflet.fig_to_geojson(fig=fig)\nwith open('test.json', 'w') as outfile:\n json.dump(js, outfile, indent=4)",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true,
"collapsed": true
},
"cell_type": "code",
"source": "%matplotlib inline\nimport matplotlib.pyplot as plt\nfrom JSAnimation import IPython_display\nfrom matplotlib.animation import FuncAnimation\n\n\ndef update_figure(k):\n global ax, fig\n ax.cla()\n \n U = (f * uamp * np.cos(v + thours[k] * omega + u - upha * np.pi/180)).sum(axis=1)\n V = (f * vamp * np.cos(v + thours[k] * omega + u - vpha * np.pi/180)).sum(axis=1)\n \n w = units['knots'] * (U + 1j * V)\n \n wf = np.NaN * np.ones_like(lonf, dtype=w.dtype)\n wf[inbox] = w\n\n # FIXME: Cannot use masked arrays and tricontour!\n # wf = ma.masked_invalid(wf)\n # cs = ax.tricontour(lonf, latf, trif, np.abs(wf).filled(fill_value=0))\n # fig.colorbar(cs)\n q = ax.quiver(lon, lat, U, V, scale=vscale)\n ax.axis(bbox) # Vineyard sound 2.\n ax.set_title('{}'.format(glocals[k]))\n\nfig, ax = plt.subplots(figsize=(10, 5))\nFuncAnimation(fig, update_figure, interval=100, frames=ntimes)",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"gist": {
"id": "",
"data": {
"description": "notebook/folium/tri_tide_movie.ipynb",
"public": true
}
},
"gist_id": "ab5133b7feecba1dfacd",
"kernelspec": {
"name": "python3",
"display_name": "Python [default]",
"language": "python"
},
"language_info": {
"name": "python",
"version": "3.6.1",
"mimetype": "text/x-python",
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"pygments_lexer": "ipython3",
"nbconvert_exporter": "python",
"file_extension": ".py"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.