{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Tutorial on how to combine different Fields into a `NestedField` object" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some applications, you may have access to different fields that each cover only part of the region of interest. Then, you would like to combine them all together. You may also have a field covering the entire region and another one only covering part of it, but with a higher resolution. The set of those fields form what we call nested fields.\n", "\n", "It is possible to combine all those fields with kernels, either with different if/else statements depending on particle position, or using recovery kernels (if only two levels of nested fields).\n", "\n", "However, an easier way to work with nested fields in Parcels is to combine all those fields into one `NestedField` object. The Parcels code will then try to successively interpolate the different fields. \n", "\n", "For each Particle, the algorithm is the following:\n", "\n", "1. Interpolate the particle onto the first `Field` in the `NestedFields` list.\n", "\n", "2. If the interpolation succeeds or if an error other than `ErrorOutOfBounds` is thrown, the function is stopped.\n", "\n", "3. If an `ErrorOutOfBounds` is thrown, try step 1) again with the next `Field` in the `NestedFields` list \n", "\n", "4. If interpolation on the last `Field` in the `NestedFields` list also returns an `ErrorOutOfBounds`, then the Particle is flagged as OutOfBounds.\n", "\n", "This algorithm means that **the order of the fields in the `NestedField` matters**. In particular, the smallest/finest resolution fields have to be listed _before_ the larger/coarser resolution fields.\n", "\n", "This tutorial shows how to use these `NestedField` with a very idealised example." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "from parcels import Field, NestedField, FieldSet, ParticleSet, JITParticle, plotTrajectoriesFile, AdvectionRK4\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First define a zonal and meridional velocity field defined on a high resolution (dx = 100m) 2kmx2km grid with a flat mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to 0.5 * cos(lon / 200 * pi / 2) m/s." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "dim = 21\n", "lon = np.linspace(0., 2e3, dim, dtype=np.float32)\n", "lat = np.linspace(0., 2e3, dim, dtype=np.float32)\n", "lon_g, lat_g = np.meshgrid(lon, lat)\n", "V1_data = np.cos(lon_g / 200 * np.pi/2)\n", "U1 = Field('U1', np.ones((dim, dim), dtype=np.float32), lon=lon, lat=lat)\n", "V1 = Field('V1', V1_data, grid=U1.grid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now define the same velocity field on a low resolution (dx = 2km) 20kmx4km grid." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "xdim = 11\n", "ydim = 3\n", "lon = np.linspace(-2e3, 18e3, xdim, dtype=np.float32)\n", "lat = np.linspace(-1e3, 3e3, ydim, dtype=np.float32)\n", "lon_g, lat_g = np.meshgrid(lon, lat)\n", "V2_data = np.cos(lon_g / 200 * np.pi/2)\n", "U2 = Field('U2', np.ones((ydim, xdim), dtype=np.float32), lon=lon, lat=lat)\n", "V2 = Field('V2', V2_data, grid=U2.grid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now combine those fields into a `NestedField` and create the fieldset" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "U = NestedField('U', [U1, U2])\n", "V = NestedField('V', [V1, V2])\n", "fieldset = FieldSet(U, V)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/h0/01fvrmn11qb62yjw7v1kn62r0000gq/T/parcels-503/ff29215bba0bd5244da08955369801d5.so\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEICAYAAACqMQjAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4HNWZ7/Hva3lhsUFiE8Y74CwGBscS2BOyKENYbzIGwuqZxNnGZAbuHZ7JvQ8kZOFCSMjcySSTSQIxgWASjHFiFg9xAoZBEBZvchy8sQhj2cLGBtwCK4CxpPf+UdVSqd3VUktd3S3593kePeo6VdX1nmqp3q5Tp+qYuyMiIpLNkFIHICIi5UtJQkREYilJiIhILCUJERGJpSQhIiKxlCRERCSWkoQMOmb2dTP7RS+Wu8PMvlOMmGK2/3szm12q7Ufi6NX+kv2T6T4JKTYz2wxUA+3AX4AlwP9099Y+vFcd8Gt3H9uHde8Amt39G31Y14HJ7t6Y77qFZGbXAce7+9+XMg4ZvHQmIaXyaXcfCUwDTgH6cqAeWvCoCqScY4saKHFK6ShJSEm5+yvA74ETAczsC2a20cx2m9kmM7s8vayZ1ZlZs5ldbWavAneH6x5jZq3hzzFmdp2Z/Tqy3kfM7GkzazGzrWb2+WyxmNmnzGxNuNzTZvZXMcs9Eb78c7jNS7LE9kszqzKzB83sNTNLha/HRt6n3sy+HJn+Ylj3lJk9ZGYTIvNOMLOlZrbLzHaETURnA18HLgnj+HO47DFmtjhcttHM/iHyPteZ2W/N7Ndm9hbw+Sz7a0Zkf/05PFtLz/t8+LnsNrOXzezv4j9dGQyUJKSkzGwccC7wp7BoJ/Ap4BDgC8APzWxaZJWjgcOACcDngHOAbe4+MvzZlvH+4wkSyX8CRwJTgTVZ4pgG3A5cDhwO/BxYbGYjMpd194+FL08Ot3lPltjmEPx//TKcHg+8A/wkZj+cR3DAvyCM848ESRAzGwU8AvwBOAY4HnjU3f8AfBe4J4zj5PDt7gaaw2UvBL5rZqdHNjcT+C1QCdyVEccY4HfAd8K6/G9gkZkdaWYHAz8GznH3UcCHs+1LGVyUJKRU7jezFuBJ4HGCgx3u/jt3f8kDjwMPAx+NrNcBfNvd97j7O73Yzt8Bj7j73e6+193fcPdsB7Z/AH7u7svdvd3d5wF7gBl51KlbbOG2Frn72+6+G7gR+HjMupcD33P3je7eRrA/poZnE58CXnX3H7j7u+6+292XZ3uTMOl+BLg6XHYN8Avgs5HFnnH3+929I8s+/HtgibsvCecvBVYRJPJ0HU80swPdfbu7r89j/8gApCQhpXKeu1e6+wR3/6f0wcrMzjGzZWFTSQvBwemIyHqvufu7eWxnHPBSL5abAHw1bGJpCbc9juDbeG91i83MDjKzn5tZU9i08wRQaWYVMdv/j8i2dwEGjMmjDoTx7gqTUlpT+D5pW3OsPwG4KGM/fAQY7e5/AS4BvgJsN7PfmdkHehmXDFBKElI2wqadRcC/AdXuXknQ88kii2V2x+upe95W4LhebH4rcGOYuNI/B7n73b0MP1ssXwXeD0x390OAdDOVsa+twOUZ2z/Q3Z/uoQ6Z29wGHBY2UaWNB17JsU5mHL/KiONgd78JwN0fcvczgNHAc8CtOd5LBgElCSknw4ERwGtAm5mdA5zZwzo7gMPN7NCY+XcBnzSzi81sqJkdbmZTsyx3K/AVM5tugYPN7H9kHGwzt3tsD7GNIrgO0WJmhwHfzrHsLcDXzOwEADM71MwuCuc9CBxtZleZ2QgzG2Vm0yNxTDSzIQDuvhV4GviemR0QXnz/EhnXHnL4NfBpMzvLzCrC96gzs7FmVm1mfxtem9gDtBJ0Y5ZBTElCykbYRPK/gIVACpgFLO5hnecILtRuCptHjsmYv4WgyeqrBE04a4CTs7zPKoLrEj8Jt90IfD7Hpq8D5oXbvDhmmR8BBwKvA8sILjzH1eM+4PvAgrBpah3BRfn0fjkD+DTwKvAi8Ilw1d+Ev98ws9Xh68uAiQRnFfcRXCdZmqMu0Ti2ElzY/jpBst4K/B+CY8UQgv24jWBffhz4p968rwxcuplOpETCrrS/cPc7Sx2LSBydSYiUgJkdRNBc9XKpYxHJRUlCpMjM7CiCZqPHCboAi5QtNTeJiEgsnUmIiEisAf9wryOOOMInTpyYdd5f/vIXDj744OIG1AuKKz+KK3/lGpviyk+ScTU0NLzu7kf2uKC7D+ifmpoaj/PYY4/FzislxZUfxZW/co1NceUnybiAVd6LY6yam0REJJaShIiIxFKSEBGRWEoSIiISS0lCRERi9TtJhE+JXBEOc7jezP5vWD7JzJab2Ytmdo+ZDQ/LR4TTjeH8iZH3+lpY/ryZndXf2EREpH8KcSaxB/gbD4ZOnAqcbWYzCJ5o+UN3n0zwVM0vhct/CUi5+/HAD8PlMLMpwKXACcDZwM9iBmcREZEi6ffNdGF/29Zwclj448DfEDzqGWAewaOVbyZ4DPF1YflvgZ+YmYXlC9x9D/CymTUCpwLP9DfGbOY2zGX+2vlJvHWPWlpaqNxcWZJt56K48lOucUH5xqa48pMrrlknzWJOzZzEYyjIHdfhN/4GggHaf0ow1GKLB2P1QjAoe3r4xDGEwye6e5uZvUkw8PwYgmfuk2WdzO3NIRhonurqaurr67PG1draGjvvZ2t+RmNrI8ePPL5XdSyk9vZ2Wlpair7dniiu/JRrXFC+sSmu/MTF1djaSEtLC+/b/b7EYyhIknD3doJB2ysJBjn5YLbFwt/Zhm70HOXZtjcXmAtQW1vrdXV1WeOqr68nbl7l5kpqK2up/3x91vlJyhVXKSmu/JRrXFC+sSmu/MTFVXdHUFaMmAvau8ndW4B6YAbBgO/pJDSWYDQrCM4QxgGE8w8lGOWqszzLOiIiUgKF6N10ZHgGgZkdCHwS2Ag8BlwYLjYbeCB8vTicJpz/3+F1jcXApWHvp0nAZGBFf+MTEZG+K0Rz02iCsX4rCJLOQnd/0Mw2EIzX+x3gT8Bt4fK3Ab8KL0zvIujRhLuvN7OFwAagDbgibMYSEZESKUTvpmeBD2Up30TQOymz/F3gopj3uhG4sb8xiYhIYeiOaxERiaUkISIisZQkREQklpKEiIjEUpIQEZFYShIiIhJLSUJERGIpSYiISCwlCRERiaUkISIisZQkREQklpKEiIjEUpIQEZFYShIiIhJLSUJERGIpSYiISKxCDF86zsweM7ONZrbezP45LL/OzF4xszXhz7mRdb5mZo1m9ryZnRUpPzssazSza/obm4iI9E8hhi9tA77q7qvNbBTQYGZLw3k/dPd/iy5sZlMIhiw9ATgGeMTM3hfO/ilwBtAMrDSzxe6+oQAxiohIHxRi+NLtwPbw9W4z2wiMybHKTGCBu+8BXg7Huk4Pc9oYDnuKmS0Il1WSEBEpkUKcSXQys4kE410vB04DrjSzzwGrCM42UgQJZFlktWa6ksrWjPLpMduZA8wBqK6upr6+Pms8ra2tsfNaWloAYucnKVdcpaS48lOucUH5xqa48hMXVzGPXwVLEmY2ElgEXOXub5nZzcANgIe/fwB8EbAsqzvZr494tm25+1xgLkBtba3X1dVljam+vp64eZWbKwFi5ycpV1ylpLjyU65xQfnGprjyExdXMY9fBUkSZjaMIEHc5e73Arj7jsj8W4EHw8lmYFxk9bHAtvB1XLmIiJRAIXo3GXAbsNHd/z1SPjqy2PnAuvD1YuBSMxthZpOAycAKYCUw2cwmmdlwgovbi/sbn4iI9F0hziROAz4LrDWzNWHZ14HLzGwqQZPRZuByAHdfb2YLCS5ItwFXuHs7gJldCTwEVAC3u/v6AsQnIiJ9VIjeTU+S/TrDkhzr3AjcmKV8Sa71RESkuHTHtYiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYmlJCEiIrGUJEREJJaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYmlJCEiIrGUJEREJFYhhi8dZ2aPmdlGM1tvZv8clh9mZkvN7MXwd1VYbmb2YzNrNLNnzWxa5L1mh8u/aGaz+xubiIj0TyHOJNqAr7r7B4EZwBVmNgW4BnjU3ScDj4bTAOcQjGs9GZgD3AxBUgG+DUwHTgW+nU4sIiJSGoUYvnQ7sD18vdvMNgJjgJlAXbjYPKAeuDosv9PdHVhmZpVmNjpcdqm77wIws6XA2cDd/Y2xrDTOZerrP4NHKkuz/Ymz4Pg5pdm2iAw4/U4SUWY2EfgQsByoDhMI7r7dzI4KFxsDbI2s1hyWxZVn284cgrMQqqurqa+vzxpPa2tr7LyWlhaA2PlJmfr6zzj4vUZaWo4v6nYBRu5tpLWlhTXN78s6P9f+KiXFlb9yjU1x5ScurmIevwqWJMxsJLAIuMrd3zKz2EWzlHmO8n0L3ecCcwFqa2u9rq4u64bq6+uJm1e5OfgmHzc/MY9U0tJyPJUXrinudgEeqaOS+Drn2l+lpLjyV66xKa78xMVVzONXQXo3mdkwggRxl7vfGxbvCJuRCH/vDMubgXGR1ccC23KUi4hIiRSid5MBtwEb3f3fI7MWA+keSrOBByLlnwt7Oc0A3gybpR4CzjSzqvCC9ZlhmYiIlEghmptOAz4LrDWzdBvK14GbgIVm9iVgC3BROG8JcC7QCLwNfAHA3XeZ2Q3AynC569MXsUVEpDQK0bvpSbJfTwA4PcvyDlwR8163A7f3NyYRESkM3XEtIiKxlCRERCSWkoSIiMRSkhARkVhKEiIiEktJQkREYilJiIhILCUJERGJpSQhIiKxlCRERCSWkoSIiMRSkhARkVhKEiIiEktJQkREYilJiIhIrEINX3q7me00s3WRsuvM7BUzWxP+nBuZ9zUzazSz583srEj52WFZo5ldU4jYRESk7woxMh3AHcBPgDszyn/o7v8WLTCzKcClwAnAMcAjZva+cPZPgTMIxrteaWaL3X1DgWKUAWj+8i384um3mb9lFZd//DhqJlSVOqTEzV++hXtWbqH6kAP2izo3NKVYtLqZ13fv4chRI7hg2thBX+eBpCBJwt2fMLOJvVx8JrDA3fcAL5tZI3BqOK/R3TcBmNmCcFklif1M+qDR0LSL519tBWDThh08vGEHxx81ki+eNolZ08eXOMrCStd51eY3eGHHX8LSNwd9nW95/CUe2bADj5TftXwLp06s4upzPqhkUQYsGE20AG8UJIkH3f3EcPo64PPAW8Aq4KvunjKznwDL3P3X4XK3Ab8P3+Zsd/9yWP5ZYLq7X5llW3OAOQDV1dU1CxYsyBpTa2srI0eOzDrvqjVXAfCjqT/Ku679MfX1q2hvb2dt9X8WdbvpbQOsOSJ7nXPtr2JpTLXzveXv0t7DcudOHMrFHxhRlJjiFGp/JVHncvgss0nHVb9lL/M2vEeuo48Bs6cMp278sKLFVW7i4irE8esTn/hEg7vX9rRcoZqbsrkZuAHw8PcPgC+SfTxsJ/v1kax/Q+4+F5gLUFtb63V1dVkDqK+vJ25e5eZKgNj5iXmkkpaWluJvN9w2xNc51/4qhoamFL9asabHgyXAks1ttB10eEmbYwqxvxqaUtyx/E8Fr3OpP8s49fX1bDvwWOZtWJszQUDwzz9vw3tsoyrxz7mc91e2uIp5/Eqsd5O773D3dnfvAG6lq0mpGRgXWXQssC1HuewH5i/fwkW3PE3Trrd7vc7DG3Zw8c+fYf7yLQlGlpx0nZtT7/R6nYc37OCSuc/Q0JRKMLLk1G/Zy7X39Zwg0pyBX+eBLrEkYWajI5PnA+meT4uBS81shJlNAiYDK4CVwGQzm2Rmwwkubi9OKj4pHw1NKb55/1o6shw5zpxSzTemH8Cif/wwp0zc95tke4fzzfvXDrgDyPzlW7g2R50X/eOHY+vc1u7c8vhLRYiysOYv35K1iWlM1YH83fTxLPrHD/Pd80/i+CMP3qe5YaDWeTAoSHOTmd0N1AFHmFkz8G2gzsymEnwZ2AxcDuDu681sIcEF6TbgCndvD9/nSuAhoAK43d3XFyI+KW/3NjTTnnHkqBhi3DDzRGZNH099fT01E6r4zVc+zPzlW/jm/Wu7Ld/usGh184C5yNnQlOJbD6wj83JgtM5p6Tpnfvv+7+d20tCUGnh1zigfWmH8+NIPddajZkIVs6aPHxR1HiwK1bvpsizFt+VY/kbgxizlS4AlhYhJBo5nX+l+FnD8USP5/mf+KuvBYNb08bz/6FFcvehZGne2dpa/vntP4nEWyrJNb9AWOYUw4Iwp1bHt7umkce39azsTS0eHD6jE2Oc6RxLFQKvzYKE7rqWkbvvjJta+srtzemiFxSaItJoJVXz/M3/F0CFdjRKPbtw5YK5NbN31l27Tl3/sWOZ+rjZnnWdNH8+N551ERVhlB+5ZsXXA1Hnj9je7Tfe6zuefxJBonVcOnDoPFkoSUjLzl2/hht9t7Jw24OLacb36plgzoYpLTunq59DuzjfuX1v2B5D5y7ewYGVz57QBow7sXRfPWdPHc+mpXU1R7T4wrsf8elkTDz77aud0vnW+LFrnAXoNaiBTkpCSaGhK8c0H1nUrqxhifGba2F6/xwXTxnY7m+hw+NYD68r2ABJcoN+3zjOOPbzX73HBtLGdZxPQdT2mXKWvRUQN9joPNkoSUhLLNr1Be6SNeojB9TNPzKu9uWZCFdfPPBGLHkA6nGWb3ihkqAVz7+pm2r3/dT79g9Xdysr5esxvG7Z268FlDP46DzZKElISr6S63w8x56PH9umxE7Omj+fyjx7bOe3A7nf29je8gmtoSrFgRVdTWMUQ4zvnndSnOl/+8eMGxPWYhqYUC1du7ZyuGGLMnjK8MHV+rjzrPBgpSUjRNTSluDty8MinjTqbzHV/8eTLZdfkFJxFdE2f/oGj+vwspmzXY8qxmW1Rw7517usjNvapc0d51nkwUpKQort7RVO3ewTybaPONOPYw7t9y2wPu0qWk82vd+/RdMSo/j13KvN6TDnWeeOrb3Wb3h/qPBgpSUhRNTSluG/1K53TFUMs7zbqTOlrE+k/Zgd+29BcNt8yG5pSPBO5TjK0Ir8L9Nl0Xo8Jp8uxzmu2tHROF7LOaeVW58FKSUKKKnp3tQGXnDKuII/Azuweureto2y+ZS5cuaXz4m0+3Xx7Mmv6eC6q7TrwtrV1lM1F+189s7nzJrhC13nmycd0TpdTnQcrJQkpmoamFPes6roWMWzokH5/u4z6TM3YbjeblcO3zIamFL+NJKtC13nquK4DbwdQddDwgr13XzU0pfivP3c9m7PQdZ4eaZoslzoPZkoSUjTRRzMYcGFNYUcgq5lQxd+ePKZzuhy+Zd67upn2juB1EnVOvf1et4fhrdv2ZuyyxbIo42xxf6jzYKYkUUINTSl++lhjXt92812nL9tIykHDKzpfO3DiMYcWfBunTDqs83Wpv2U2NKVYmOCZEwQX7YdF7jQr9dnT/ljnwS7JQYckh4amFJf+/Bn2djhDjF71mf+PR1/gR0tfxIHhFcbdc/465ze0uU+8xE2/f44Oh6FDjKdq91B9SOlGc3t4fdejGYYQfCMstPS3zHR7eCm/ZS7b9AZ725M7c4Lg7Omi2nHcFd4zkL4Wc2aJnoF37+rmRM8WIb7OevBfMnQmUSL/+oeN7A3/mTo8eNplrpuD5j7xEj8MEwTAez08X//Xz2zmu0ue67xg2tbhvPx6K7vfbStQDfLT0JRi2aZdndNDhw7pV7fXOOX0LbMqcv9GUmdOEHQNTdc5fS2mMdWbse4Kq6EpxW8SPotIi3aHLZfrT4OVkkQJvNPmLH+5+x+0E//coYamFN/7/XP7lKefr59t+W8tzj4Ux2utpXmcwaMbd3Tr7ZLEN0zo+paZ1t5euusSG1/terptUmdOENT54ow6P7er+EmiGGdOadnqXOrrT4OVkkQJ7Ho3++CNcTcHLVi5ZZ8BaqDr+fqZ7m1ozjriGcDe9FXUItvWEgzRaQYjhiX3DRP2/WZdqusSTeEjwYcYDB+WzJlTWrcbzcwYOSzbUPLJmnHs4Z3P0RpWkexnDGFvtkid1cspGQVJEmZ2u5ntNLN1kbLDzGypmb0Y/q4Ky83MfmxmjWb2rJlNi6wzO1z+RTObXYjYys3ud9vY/V73wVfSsp02NzSluHd190dL57qBqqEpxT0NXaf8QwwiN6nS8vbeop+WNzSlWBx2iaww41ufOiHR9uOaCVV861NTgKAp7/oH1xe/zpt38ccXXgdgSJHq/H/Oej8QfNmY/9x7Ra/ze23tkS8nvR3Fuu9qJlTxlY9NAoIvTKX4nPcHhTqTuAM4O6PsGuBRd58MPBpOA5xDMK71ZGAOcDMESYVg2NPpwKnAt9OJZTB5LfL0SgMumz6ev4vcBJbZbTN4Wmr35aMXuDNvGlu26Q3aIqf8l57a/Xn87sV/lMGCFV03k7l7Ys0uUW9Frr3s2Vv8G+tuf6rrZrJi1Tk68tvejuI/TvuWxzd1vi7W03gPGhH0vXHgvb1qckpCQZKEuz8B7MoongnMC1/PA86LlN/pgWVApZmNBs4Clrr7LndPAUvZN/EMaA1NqW7XBNIX9k4Y03VBM7PbZvQibPri5wXTxlIxpKssejZRmeVi6QXTxjK0RBdzG5pS3PunroNVRUWyzS5pM449PHYfJa2hKcUf1m3vnC5qnUt0M2FDU4o/vvBa53Sx6lx1UFdvvVJ3eR6skuwCW+3u2wHcfbuZHRWWjwG2RpZrDsviyvdhZnMIzkKorq6mvr4+awCtra2x81pagufKxM1PwoMvvceF3vUt/7Sjjd0v/5nVL3Xvtrl05QaOeSf4VvbbZe90rm/A6nXPccw7w/nr0RU8+UpwcXJvWwd3P7KS3ccNZ9Hqd7Mu/5HRXd8H9rZ18JP/Ws7sEw/oFl+u/dVX89bt6XYmlK5z/cu9f4++xnXa6AqeyLKPCiUurnnr9nS7mawvde6rU4+u4JntydU5zi/XvUv6aleuOhf6b2z1S11naF1/75viV4iRxN9+IcTFVczjVynuk8h2Rc1zlO9b6D4XmAtQW1vrdXV1WTdUX19P3LzKzZUAsfOTMGpSivaHjQ53RgwbwpWfnk7NhCpGTUqxeNMzvBceWZ561bny0ycD8OIfnu5cf9jQIVz2yVOomVDFtgO38OR9a4FgJ0078QOMOnoUax7KvvyoSSn2PkTn8ultRNvJc+2vvmhoSvHU0me6xZOucz76Gte2A7fwRMY+qivAc6JyxdXQlOLJAtS5r7YduIVnEqxzNg1NKZ56uHd1LvTf2KhJKR7cvIy9bR3d/t7zVei4CiUurmIev5JMEjvMbHR4FjEa2BmWNwPjIsuNBbaF5XUZ5fUJxld0NROq2D36EHamdnPXl2d0/jHH3Rz0zp622G6j2W4aW/tKS7cHyUWXr5lQxYsHj+CNvwTNXekug0kevIrZJTKb6D4ykuuCGpV5TagUdU4bYkWsc8I30MWpmVDFXV+ewbJNbzDj2MN1Q10CkkwSi4HZwE3h7wci5Vea2QKCi9RvhonkIeC7kYvVZwJfSzC+khh1wFDaDzSOy/hjvmDaWBau2sredseBhau20hEZsSXzxqT0TWPps4+Fq7biHfHLAxx96AGdSaK/Yzj0RrpLpHtxukRm2/6IoUN4t60DrDjt1dF9muTNZLm2P2LoEPa0BY0/xajzoQd0HUaS7t6cTc2EKiWHBBWqC+zdwDPA+82s2cy+RJAczjCzF4EzwmmAJcAmoBG4FfgnAHffBdwArAx/rg/L9guZNwe1tXu3Nt7Mb2eZN421tXsvH6rWdf9A0t5+r62oXSIz1Uyo4lufPiE4myhSV9jtb77TVdNsN7ckrGZCFd/+9AlAcbr/NjSluP7BjUBw5pJ0V18pvkL1brrM3Ue7+zB3H+vut7n7G+5+urtPDn/vCpd1d7/C3Y9z95PcfVXkfW539+PDn18WIraBJNprKSruG2n0prGouAFe3np3L+mD9d725LvC3lzf2Pm6WF0iM0WbW4rRFXZuCbqBZipmnZdteoP3IjdoFqN5S4pLd1yXkZoJVVmf7xN3VpB59pE2ZfQhWZc/5IBhmBWnK2xDU4pnXuo6ESxWl8hM0Wc5Jd0ttGHzLp59peuBgqWsc7G6wo46oOvJvh2uLqiDkZJEmbnklO49UYb1MOxj8DiG3O+RNuqAoRwZGWc4ydHb7lrWVJRnNfUks1kuyTrf+seus4hS1/ljY7quEyRZ54fX7+x8neTzqaR0lCTKzKzp4/nu+Sdx8thDOXNKNQt6eBx4zYQq7rn8w5wxpZqTxx7Kd8/P/cjxI0eOSPybdUNTigfWdI1jXYoLuFHFeGJoQ1OKhzfs6JwudZ1PGzM08ZsJG5pSPNX4eud0Uk/2ldLSeBJlaFbGozd6UjOhils/V9urZUcdMJSLE34Wf+Y41qX6Rp2WbpabvyKoc/rRJ4WMKXMc61LX+fiqCi6YNpbfrArOIJKoc+Y41qWusyRDZxL7oVyP9eivpMex7qsTczz6pL+SHse6rz6U4PjXSY9jLeVDSWI/VDOhiotquv6hCzkWdClvrMolyXGRl216PdFxrPsqyTr/ZtXWsjpblOQoSeyn/mps92+Zu9/ZW5D3fSX1dufrJEdjy1fmiHX3rNiacyTAfLywo7Xz9f5Q52KOQCelpySxn8rshfKLJ1/u95CXDU0p7l7ZdfAo1qMwemOfEevcY0cCzEdDU4rFa7qaXfaHOkevOQHUve9InUUMYkoS+6kZxx7eNZIZwY1fT73Sv7OJ+Suaut1kXIxHf+Sj2+htxI8EmI9fPvVyt3vJ94c6b9zxVrfpIyLdqmXwUZLYT9VMqOL6mSd2jlrnwOOvtPe5OWL+8i3c29DV7bViiHH9zBPL6htmZ53DaQfuWdn3Jpj6LXt58NnIuBFlXOfoaIYLV23t89nE/GVNrG5q6ZyOu7tfBg8lif3YrOndR63rcPrUHNHQlOKbD6zr1h3yklPG5dWNt1gy69ze0bcmmIamFHdu7D6WQTnX+fQPHtk53dbHR7I0NKX4xgOdIxRjwMW148oqKUomL35NAAAQzElEQVThKUns54LusN2bI/Lt6RQMsdrV6FIxpLy/XV5Q0/8639vQHHl4YfnXufqQA7tNvx4ZRre3Fq7aMqDqLIWhJLGfq5lQxT98ZFLntJN/T6fnX+3eRv3lj0wq62+X/a1zQ1OKBSu7mqjKsZkpU+a1iUc37sirma2hKcVvVkaGoR0AdZbCUJIQRh04rFt/+rl/3NTrA8gdT73M4j93tctb+H7lLjPGfOr840df7Na75/QPHFWWzUxRNROquOSUaE8n+Mb9a3td5x889DwdkemBUGcpDCUJCZ4aGvmW2dHLA8j8ZU1c918bupWVW++eOJm9u3pb518/s5nHX3itW9lA6d2TeTbR4fDN+9f2eD3m9ic38XRGc9xAqbP0n5KE7NMDBoIDyLX3xR805y/fwtfvX9etbIgxYJogMnt3Qc91/t6SjXzjgfXdyiqMAdMun+1zbne4etGzsYniOw+u7xxUKG0g1Vn6L/EkYWabzWytma0xs1Vh2WFmttTMXgx/V4XlZmY/NrNGM3vWzKYlHZ8EZk0fz+wpw7sdQBz4+n1rmXPnqs6DSENTin+Yt5Kv37e22/pm8J3zcj+BttzMmj6e75x3Uq/qfOHNT/HzJzZ1W9+AG847aUAkxbRZ08dzxpTqbmWNO1u58JanuWlJVzJoaErxmZ89yS+e3Nxt2SE28Oos/VOsp8B+wt1fj0xfAzzq7jeZ2TXh9NXAOcDk8Gc6cHP4W4qgbvwwtlHV7ZHXAA9v2MHSDTt4/9GjeO7V3VnX/eQHqwdUgkibNX089c/vjK3z+MMOomnX2/usZ8DsKcMHZJ0v//hxPPrczm490tzhlic28fgLr9GB8/yrrfusZwy8LwLSf6VqbpoJzAtfzwPOi5TfGQ5xugyoNLPRpQhwf3X5x49jaJYhUR1iE8TQCuMrHz8u4ciSk6vO2RIEwOUfO5a68eV/gT6bmglV3DDzRLJUmY2v7s6aICCosxLE/sc84cHazexlIEXwP/dzd59rZi3uXhlZJuXuVWb2IHCTuz8Zlj8KXB0dBzssnwPMAaiurq5ZsGBB1m23trYycuTIrPOuWnMVAD+a+qN+1jA/U1+/ivb2dtZW/2dRt5veNsCaI7LXOb2/GlPtLNn0Hqtf68i6XNS0oyo4d9Iwjq+q6HHZvsr1ORZKb+tswDkTh3LxB0YUJa6+6k1sjal2Fj6/hxdach8DonUuRlylMNDiKsTx6xOf+ESDu/c4EE0xmptOc/dtZnYUsNTMnsuxbJbvNuzzF+zuc4G5ALW1tV5XV5f1zerr64mbV7k5yFFx8xPzSCUtLS3F3264bYivc3p/1QFfJrg4/Y3713a7gSptSBGvQeT6HAuljp7rfOrEKq4+54Od7fHFiKuvehNbHUGdb1qykZ8/sWnffzT2rXMx4iqFgRZXMY9fiScJd98W/t5pZvcBpwI7zGy0u28Pm5PSA+U2A+Miq48FtiElMWv6eN5/9CgWrW7udofukaNGcMG0wTl+QGadW95+jz1tHVxySn6jBQ4k15z7Qc444ej9qs7Se4kmCTM7GBji7rvD12cC1wOLgdnATeHvB8JVFgNXmtkCggvWb7r79n3fWYqlZkLVoEwGuajOIl2SPpOoBu4zs/S25rv7H8xsJbDQzL4EbAEuCpdfApwLNAJvA19IOD4REckh0STh7puAk7OUvwGcnqXcgSuSjElERHpPd1yLiEgsJQkREYmlJCEiIrGUJEREJJaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYmlJCEiIrGUJEREJJaShIiIxFKSEBGRWEoSIiISq+yShJmdbWbPm1mjmV1T6nhERPZnZZUkzKwC+ClwDjAFuMzMppQ2KhGR/VfSw5fm61SgMRzRjnCs65nAhpJGVUg7H6cS4JG64m87tQaqphZ/uyIyYJVbkhgDbI1MNwPTMxcysznAHIDq6mrq6+uzvllra2vsvJaWFoDY+Umpy9h+UdlEduypYXsf9lcpKa78lWtsiis/cXEV8/hVbknCspT5PgXuc4G5ALW1tV5XV5f1zerr64mbV7m5EiB2fnI8Z1xJqwTeHzOvlHHlorjyV66xKa78xMVVzONXWV2TIDhzGBeZHgtsK1EsIiL7vXJLEiuByWY2ycyGA5cCi0sck4jIfqusmpvcvc3MrgQeAiqA2919fYnDEhHZb5VVkgBw9yXAklLHISIi5dfcJCIiZURJQkREYilJiIhILCUJERGJpSQhIiKxlCRERCSWkoSIiMRSkhARkVhKEiIiEktJQkREYilJiIhILCUJERGJpSQhIiKxlCRERCSWkoSIiMRSkhARkViJJQkzu87MXjGzNeHPuZF5XzOzRjN73szOipSfHZY1mtk1ScUmIiK9k/TIdD9093+LFpjZFIKxq08AjgEeMbP3hbN/CpwBNAMrzWyxu29IOEYREYlRiuFLZwIL3H0P8LKZNQKnhvMa3X0TgJktCJdVkhARKZGkk8SVZvY5YBXwVXdPAWOAZZFlmsMygK0Z5dOzvamZzQHmAFRXV1NfX591462trbHzWlpaAGLnJylXXKWkuPJTrnFB+camuPITF1cxj1/9ShJm9ghwdJZZ1wI3AzcAHv7+AfBFwLIs72S/PuLZtuvuc4G5ALW1tV5XV5c1vvr6euLmVW6uBIidn6RccZWS4spPucYF5Rub4spPXFzFPH71K0m4+yd7s5yZ3Qo8GE42A+Mis8cC28LXceUiIlICSfZuGh2ZPB9YF75eDFxqZiPMbBIwGVgBrAQmm9kkMxtOcHF7cVLxiYhIz5K8JvGvZjaVoMloM3A5gLuvN7OFBBek24Ar3L0dwMyuBB4CKoDb3X19gvGJiEgPEksS7v7ZHPNuBG7MUr4EWJJUTCIikh/dcS0iIrGUJEREJJaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYmlJCEiIrGUJEREJJaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYnVryRhZheZ2Xoz6zCz2ox5XzOzRjN73szOipSfHZY1mtk1kfJJZrbczF40s3vCIUxFRKSE+nsmsQ64AHgiWmhmUwjGqD4BOBv4mZlVmFkF8FPgHGAKcFm4LMD3gR+6+2QgBXypn7GJiEg/9Wv4UnffCGBmmbNmAgvcfQ/wspk1AqeG8xrdfVO43gJgppltBP4GmBUuMw+4Dri5P/Hl8njT4wDU3VGX1CZitbS0ULm5sujb7Yniyk+5xgXlG5viyk9cXGteXcPUo6cWJYakxrgeAyyLTDeHZQBbM8qnA4cDLe7elmX5fZjZHGAOQHV1NfX19VmXa21tjZ2X1tLSknN+Etrb20uy3Z4orvyUa1xQvrEprvzExTXxgInUDK/p8fhWCD0mCTN7BDg6y6xr3f2BuNWylDnZm7c8x/JZuftcYC5AbW2t19XVZV2uvr6euHleF/v2icsVVykprvyUa1xQvrEprvyUQ1w9Jgl3/2Qf3rcZGBeZHgtsC19nK38dqDSzoeHZRHR5EREpkaS6wC4GLjWzEWY2CZgMrABWApPDnkzDCS5uL3Z3Bx4DLgzXnw3EnaWIiEiR9LcL7Plm1gz8NfA7M3sIwN3XAwuBDcAfgCvcvT08S7gSeAjYCCwMlwW4GviX8CL34cBt/YlNRET6r7+9m+4D7ouZdyNwY5byJcCSLOWb6OoBJSIiZUB3XIuISCwlCRERiaUkISIisZQkREQklgW9TwcuM3sNaIqZfQTBPRjlRnHlR3Hlr1xjU1z5STKuCe5+ZE8LDfgkkYuZrXL32p6XLC7FlR/Flb9yjU1x5acc4lJzk4iIxFKSEBGRWIM9ScwtdQAxFFd+FFf+yjU2xZWfksc1qK9JiIhI/wz2MwkREekHJQkREYk1oJOEmf0/M3vOzJ41s/vMrDIsn2hm75jZmvDnlsg6NWa21swazezHFo69amaHmdlSM3sx/F2VUMxnm9nz4favSWIbGdsbZ2aPmdlGM1tvZv8cll9nZq9E9tG5kXW+Fsb3vJmdlVTsZrY5/CzWmNmqsCzr52CBH4fbftbMpkXeZ3a4/ItmNrufMb0/sk/WmNlbZnZVKfaXmd1uZjvNbF2krGD7J+5/oY9xlfx/MSaugn1uFgxxsDyM6x4Lhjvoa1z3RGLabGZrir2/es3dB+wPcCYwNHz9feD74euJwLqYdVYQPNrcgN8D54Tl/wpcE76+Jv1eBY63AngJOBYYDvwZmJLwPhoNTAtfjwJeAKYQjCH+v7MsPyWMawQwKYy3IonYgc3AERllWT8H4Nzw8zJgBrA8LD8M2BT+rgpfVxXw83oVmFCK/QV8DJgW/Vsu5P6J+1/oY1wl/1+MiatgnxvB8AeXhq9vAf6xr3FlzP8B8K1i76/e/gzoMwl3f9i7xsVeRjCiXSwzGw0c4u7PeLBH7wTOC2fPBOaFr+dFygvpVKDR3Te5+3vAgnC7iXH37e6+Ony9m2Acj9jxw8N4Frj7Hnd/GWgM4y5W7HGfw0zgTg8sIxjJcDRwFrDU3Xe5ewpYCpxdoFhOB15y97g7+tNxJbK/3P0JYFeW7fV7//Twv5B3XOXwvxizv+Lk9bmF39r/BvhtIeMK3/di4O5c71HKY9eAThIZvkiQXdMmmdmfzOxxM/toWDaGYGjVtGa6DpjV7r4dggMrcFQCMY4BtsZsP3FmNhH4ELA8LLoybB64PXKKGhdjErE78LCZNZjZnLAs7nMoZlxpl9L9n7fU+wsKt39y/S/0V7n9LxbiczscaIkkwkLtr48CO9z9xUhZqfdXN2WfJMzsETNbl+VnZmSZa4E24K6waDsw3t0/BPwLMN/MDiE4TctUzD7AJdu+mY0EFgFXuftbwM3AccBUgv31gx5iTCL209x9GnAOcIWZfSzHssWMi7C9+W+B34RF5bC/csk3jqT2W7n9Lxbqc0sq3svo/kWk1PtrH/0ama4Y3P2TueZbcCHuU8Dp4WkY7r4H2BO+bjCzl4D3EWTf6GnwWGBb+HqHmY129+3hqd3OwtYEwu2Pi9l+YsxsGEGCuMvd7wVw9x2R+bcCD/YixoLG7u7bwt87zew+glP9uM8hLq5moC6jvL4/cYXOAVan91M57K9QofZPrv+FPinH/8UCfm6vEzThDQ3PJgqxv4YCFwA1kXjL7thV9mcSuZjZ2QRjY/+tu78dKT/SzCrC18cCk4FN4anYbjObEbYFfg54IFxtMZDu+TE7Ul5IK4HJYS+J4QTNGYsT2E6nsJ63ARvd/d8j5aMji50PpHteLAYuNbMRZjaJYN+tKHTsZnawmY1Kvya48LmO+M9hMfA5C8wA3gw/z4eAM82sKmxKODMs669u3/BKvb8iCrJ/evhfyFu5/i8W6nMLk95jwIWFiCv0SeA5d+9sRir1/sqqkFfBi/1DcLFpK7Am/LklLP8MsJ6gZ8Jq4NORdWoJ/lBeAn5C113nhwOPAi+Gvw9LKOZzCXoYvQRcW4R99BGC09JnI/vpXOBXwNqwfDEwOrLOtWF8zxPp8VLI2Al6j/w5/Fmffr+4z4HgdPun4bbXArWR9/pi+LfQCHyhAPvsIOAN4NBIWdH3F0GS2g7sJfgm+aVC7p+4/4U+xlXy/8WYuAr2uYV/syvCuv4GGNHXuMLyO4CvZCxbdscuPZZDRERiDejmJhERSZaShIiIxFKSEBGRWEoSIiISS0lCRERiKUmIiEgsJQkREYn1/wFIbEZYMCCTCwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[1000])\n", "output_file = pset.ParticleFile(name='NestedFieldParticle.nc', outputdt=50)\n", "pset.execute(AdvectionRK4, runtime=14000, dt=10, output_file=output_file)\n", "plt = plotTrajectoriesFile('NestedFieldParticle.nc', show_plt=False)\n", "plt.plot([0,2e3,2e3,0,0],[0,0,2e3,2e3,0], c='orange')\n", "plt.plot([-2e3,18e3,18e3,-2e3,-2e3],[-1e3,-1e3,3e3,3e3,-1e3], c='green');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we observe, there is a change of dynamic at lon=2000, which corresponds to the change of grid.\n", "\n", "The analytical solution to the problem:\n", "\\begin{align}\n", "dx/dt &= 1;\\\\\n", "dy/dt &= \\cos(x \\pi/400);\\\\\n", "\\text{with } x(0) &= 0, y(0) = 1000\n", "\\end{align}\n", "\n", "is\n", "\\begin{align}\n", "x(t) &= t;\\\\\n", "y(t) &= 1000 + 400/\\pi \\sin(t \\pi / 400)\n", "\\end{align}\n", "which is captured by the High Resolution field (orange area) but not the Low Resolution one (green area)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Keep track of the field interpolated\n", "\n", "For different reasons, you may want to keep track of the field you have interpolated. You can do that easily by creating another field that share the grid with original fields.\n", "Watch out that this operation has a cost of a full interpolation operation." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "F1 = Field('F1', np.ones((U1.grid.ydim, U1.grid.xdim), dtype=np.float32), grid=U1.grid)\n", "F2 = Field('F2', 2*np.ones((U2.grid.ydim, U2.grid.xdim), dtype=np.float32), grid=U2.grid)\n", "F = NestedField('F', [F1, F2])\n", "fieldset.add_field(F)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO: Compiled SampleParticleSampleNestedFieldIndex ==> /var/folders/h0/01fvrmn11qb62yjw7v1kn62r0000gq/T/parcels-503/57e7f8f8b9c55074a13384b46fceb0eb.so\n", "WARNING: dt or runtime are zero, or endtime is equal to Particle.time. The kernels will be executed once, without incrementing time\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Particle (1000, 500) interpolates Field #1\n", "Particle (10000, 500) interpolates Field #2\n" ] } ], "source": [ "from parcels import Variable\n", "\n", "def SampleNestedFieldIndex(particle, fieldset, time):\n", " particle.f = fieldset.F[time, particle.depth, particle.lat, particle.lon]\n", "\n", "class SampleParticle(JITParticle):\n", " f = Variable('f', dtype=np.int32)\n", " \n", "pset = ParticleSet(fieldset, pclass= SampleParticle, lon=[1000], lat=[500])\n", "pset.execute(SampleNestedFieldIndex, runtime=0, dt=0)\n", "print('Particle (%g, %g) interpolates Field #%d' % (pset[0].lon, pset[0].lat, pset[0].f))\n", "\n", "pset[0].lon = 10000\n", "pset.execute(SampleNestedFieldIndex, runtime=0, dt=0)\n", "print('Particle (%g, %g) interpolates Field #%d' % (pset[0].lon, pset[0].lat, pset[0].f))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "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.15" } }, "nbformat": 4, "nbformat_minor": 2 }