Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Cell Tree, Index translation, and Bilinear Interpolation #78

Merged
merged 12 commits into from
Feb 9, 2016
Merged

Cell Tree, Index translation, and Bilinear Interpolation #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
dcb67bb
Reset repo from clean fork.
jay-hennen Feb 1, 2016
85668ab
begin refactoring for new nodes format
jay-hennen Feb 1, 2016
d31bacf
major refactoring to deal with API changes
jay-hennen Feb 2, 2016
9058a83
Refinement of interpolation and beginning of renaming
jay-hennen Feb 3, 2016
91e7e06
added edge1/2 lon/lat to constructor
jay-hennen Feb 3, 2016
b2bf440
Bugfixes, renaming complete
jay-hennen Feb 3, 2016
5799d53
Factored variable mask into calculations
jay-hennen Feb 3, 2016
12b236f
Docstring formatting and renamed XtoL
jay-hennen Feb 5, 2016
6b77f02
Added demo notebook, fixed minor bugs in Sgrid interpolation
jay-hennen Feb 8, 2016
bc55e06
refactored translation and interpolation alphas function signature
jay-hennen Feb 8, 2016
f8f7bcb
added timing to interpolation demo
jay-hennen Feb 8, 2016
2a5fde1
added timing in the correct place
jay-hennen Feb 8, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
370 changes: 370 additions & 0 deletions notebooks/Sgrid Interpolation Demo.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,370 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pysgrid Interpolation Demo \n",
" \n",
"Lets start with the pysgrid object."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from netCDF4 import Dataset\n",
"import pysgrid\n",
"import numpy as np\n",
"\n",
"url = ('http://geoport.whoi.edu/thredds/dodsC/clay/usgs/users/'\n",
" 'jcwarner/Projects/Sandy/triple_nest/00_dir_NYB05.ncml')\n",
"\n",
"sgrid = pysgrid.load_grid(Dataset(url))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Pysgrid variable objects can now serve data, similar to pyugrid variables. It calls directly on the netCDF Variable object, and has it's lazy-loading behavior."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SGridVariable object: sea_water_potential_temperature, on the faces\n",
"(201, 16, 107, 345)\n"
]
}
],
"source": [
"print sgrid.temp\n",
"print sgrid.temp.shape"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 18.07930946, 18.14973068, 18.2201519 , 18.21634102,\n",
" 18.21252823],\n",
" [ 18.15918732, 18.22960854, 18.30002975, 18.29621696,\n",
" 18.29240417],\n",
" [ 18.23906708, 18.30948639, 18.37990761, 18.37609482,\n",
" 18.32413864],\n",
" [ 18.31894493, 18.38936615, 18.45978546, 18.40782738,\n",
" 18.35586929],\n",
" [ 18.45632744, 18.52674675, 18.50500488, 18.4530468 ,\n",
" 18.40108871]], dtype=float32)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sgrid.temp[0,0,50:55,100:105]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First lets establish some points on the grid we're interested in. For simplicity, we'll derive them from the grid nodes."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-73.74568073 40.32299155]\n",
" [-73.7391676 40.32583585]\n",
" [-73.73265448 40.32868015]\n",
" [-73.72614136 40.33152445]\n",
" [-73.71962823 40.33436874]\n",
" [-73.74901747 40.33079974]\n",
" [-73.7425043 40.33364399]\n",
" [-73.73599114 40.33648823]\n",
" [-73.72947797 40.33933248]\n",
" [-73.72296481 40.34217672]]\n"
]
}
],
"source": [
"points_lon = sgrid.node_lon[50:52,100:105] * 0.9999\n",
"points_lat = sgrid.node_lat[50:52,100:105] * 0.9999\n",
"points = np.stack((points_lon, points_lat), axis = -1).reshape(-1,2)\n",
"print points"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The SGrid object now provides a way to get the x,y indices of these points on the node grid. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 49 100]\n",
" [ 49 101]\n",
" [ 49 102]\n",
" [ 49 103]\n",
" [ 49 104]\n",
" [ 50 100]\n",
" [ 50 101]\n",
" [ 50 102]\n",
" [ 50 103]\n",
" [ 50 104]]\n"
]
}
],
"source": [
"node_ind = sgrid.locate_faces(points)\n",
"print node_ind"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Scenario: I want to interpolate the temperature and salinity at a list of points."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[23.40559527721179 23.391546600118488 23.41423669831296 23.436921401441918\n",
" 23.45776700834991 23.419357486410757 23.4053081701454 23.427996148979414\n",
" 23.44884449221559 23.45918945470487]\n",
"[43.41301646604848 43.40833575274081 43.4036551192405 43.39897456556031\n",
" 43.394294092083555 43.41151307230791 43.40683296250011 43.4021529311478\n",
" 43.39747298010166 43.39279310876989]\n",
"Wall time: 18 ms\n"
]
}
],
"source": [
"%%time\n",
"interp_temps = sgrid.interpolate_var_to_points(points, sgrid.temp[0,6])\n",
"interp_salinity = sgrid.interpolate_var_to_points(points, sgrid.salt[0,6])\n",
"\n",
"print interp_temps\n",
"print interp_salinity"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Done!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Interpolating for a vector quantity like velocity is a bit more complicated, but the same approach can be used piecewise."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.02588658 -0.19636827]\n",
" [ 0.02969118 -0.19708817]\n",
" [ 0.03296548 -0.19143968]\n",
" [ 0.03621145 -0.18579166]\n",
" [ 0.03945695 -0.1801441 ]\n",
" [ 0.03022276 -0.19490993]\n",
" [ 0.03349701 -0.19477434]\n",
" [ 0.03674293 -0.18912618]\n",
" [ 0.03998839 -0.18347849]\n",
" [ 0.04323336 -0.17783126]]\n",
"Wall time: 22 ms\n"
]
}
],
"source": [
"%%time\n",
"interp_u = sgrid.interpolate_var_to_points(points, sgrid.u[0,6])\n",
"interp_v = sgrid.interpolate_var_to_points(points, sgrid.v[0,6])\n",
"print np.column_stack((interp_u, interp_v))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Trying to interpolate to points outside the grid will return nans."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-1 -1]\n",
" [-1 -1]]\n"
]
}
],
"source": [
"bad_pts = np.array(([0,0],[1,1]))\n",
"bad_ind = sgrid.locate_faces(bad_pts)\n",
"print bad_ind"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[nan nan]\n"
]
}
],
"source": [
"interp_u = sgrid.interpolate_var_to_points(bad_pts, sgrid.u[0,6])\n",
"print interp_u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you are interpolating multiple variables at the same time or interpolating to the same points repeatedly, you can save considerable computation by pre-computing some information for interpolate_var_to_points."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[23.40559527721179 23.391546600118488 23.41423669831296 23.436921401441918\n",
" 23.45776700834991 23.419357486410757 23.4053081701454 23.427996148979414\n",
" 23.44884449221559 23.45918945470487]\n",
"[43.41301646604848 43.40833575274081 43.4036551192405 43.39897456556031\n",
" 43.394294092083555 43.41151307230791 43.40683296250011 43.4021529311478\n",
" 43.39747298010166 43.39279310876989]\n",
"Wall time: 13 ms\n"
]
}
],
"source": [
"%%time\n",
"indices = sgrid.locate_faces(points)\n",
"\n",
"#re-use these variables if you plan to interpolate to the same points repeatedly\n",
"center_inds = sgrid.translate_index(points, indices, 'center')\n",
"interp_alphas = sgrid.interpolation_alphas(points, center_inds, 'center')\n",
"\n",
"interp_temps = sgrid.interpolate_var_to_points(points, sgrid.temp[0,6], center_inds, interp_alphas)\n",
"interp_salinity = sgrid.interpolate_var_to_points(points, sgrid.salt[0,6], center_inds, interp_alphas)\n",
"\n",
"print interp_temps\n",
"print interp_salinity"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.11"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Loading