removed dependencies on segtrax code

Author: Andy Barrett

dataset-diagram.png

@@ -24,12 +24,9 @@
     "import numpy as np\n",
     "import glob\n",
     "import sys\n",
+    "import pyproj\n",
     "\n",
-    "from scipy.interpolate import interp2d\n",
-    "\n",
-    "# Our modules\n",
-    "sys.path.append('/home/jovyan/segtrax/source')\n",
-    "from utilities import transform_coord\n"
+    "import scipy.interpolate as scinterp\n"
    ]
   },
   {
@@ -444,7 +441,7 @@
     {
      "data": {
       "text/plain": [
-       "<matplotlib.collections.QuadMesh at 0x7f68644f0c88>"
+       "<matplotlib.collections.QuadMesh at 0x7f6e74668b38>"
       ]
      },
      "execution_count": 13,
@@ -483,6 +480,33 @@
    "cell_type": "code",
    "execution_count": 14,
    "metadata": {},
+   "outputs": [],
+   "source": [
+    "def transform_coord(from_epsg, to_epsg, x, y):\n",
+    "    \"\"\"Transform coordinates from proj1 to proj2 (EPSG num).\n",
+    "    \n",
+    "    from_epsg - EPSG code for from_proj\n",
+    "    to_epsg - EPSG code for to_proj\n",
+    "    x - x-coordinate to convert\n",
+    "    y - y-coordinate to convert\n",
+    "    \n",
+    "    Useful EPSG:\n",
+    "    4326 - WGS84\n",
+    "    3408 - EASE North (https://spatialreference.org/ref/epsg/3408/)\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Set full EPSG projection strings\n",
+    "    from_proj = pyproj.Proj(\"+init=EPSG:\"+str(from_epsg))\n",
+    "    to_proj = pyproj.Proj(\"+init=EPSG:\"+str(to_epsg))\n",
+    "    \n",
+    "    # Convert coordinates\n",
+    "    return pyproj.transform(from_proj, to_proj, x, y)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
    "outputs": [
     {
      "name": "stdout",
@@ -508,7 +532,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
@@ -542,7 +566,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [
     {
@@ -578,129 +602,50 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "<xarray.DataArray 'u' (time: 1, i: 3)>\n",
-      "dask.array<shape=(1, 3), dtype=float32, chunksize=(1, 3)>\n",
-      "Coordinates:\n",
-      "  * time     (time) object 2018-11-05 12:00:00\n",
-      "    x        (i) float64 0.0 3.876e+04 5.792e+04\n",
-      "    y        (i) float64 -5.558e+05 -4.43e+05 -3.285e+05\n",
-      "Dimensions without coordinates: i\n",
-      "Attributes:\n",
-      "    short_name:     U\n",
-      "    long_name:      sea_ice_x_velocity\n",
-      "    standard_name:  sea_ice_x_velocity\n",
-      "    units:          cm / s\n",
-      "    valid_min:      -200.0\n",
-      "    valid_max:      200.0\n",
-      "    grid_mapping:   crs\n",
-      "    comment:        U is the along-x component of the ice motion.  It is *not... <xarray.DataArray 'u' (time: 1, i: 3)>\n",
-      "dask.array<shape=(1, 3), dtype=float32, chunksize=(1, 3)>\n",
-      "Coordinates:\n",
-      "  * time     (time) object 2018-11-05 12:00:00\n",
-      "    x        (i) float64 0.0 3.876e+04 5.792e+04\n",
-      "    y        (i) float64 -5.558e+05 -4.43e+05 -3.285e+05\n",
-      "Dimensions without coordinates: i\n",
-      "Attributes:\n",
-      "    short_name:     U\n",
-      "    long_name:      sea_ice_x_velocity\n",
-      "    standard_name:  sea_ice_x_velocity\n",
-      "    units:          cm / s\n",
-      "    valid_min:      -200.0\n",
-      "    valid_max:      200.0\n",
-      "    grid_mapping:   crs\n",
-      "    comment:        U is the along-x component of the ice motion.  It is *not...\n"
+      "[[-1.94814917 -4.22065543 -4.93430168]] [[-6.13456718 -5.56727877 -4.18583824]]\n"
      ]
     }
    ],
    "source": [
     "ui = ds.sel(time='2018-11-05').u.interp(x=ind_x, y=ind_y)\n",
-    "vi = ds.sel(time='2018-11-05').u.interp(x=ind_x, y=ind_y)\n",
-    "print (ui, vi)"
+    "vi = ds.sel(time='2018-11-05').v.interp(x=ind_x, y=ind_y)\n",
+    "print (ui.values, vi.values)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Johan's approach"
+    "### Johan's approach"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": 36,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "<xarray.DataArray 'u' (y: 361, x: 361)>\n",
-       "dask.array<shape=(361, 361), dtype=float32, chunksize=(361, 361)>\n",
-       "Coordinates:\n",
-       "  * x        (x) float64 -4.512e+06 -4.487e+06 ... 4.487e+06 4.512e+06\n",
-       "  * y        (y) float64 -4.512e+06 -4.487e+06 ... 4.487e+06 4.512e+06\n",
-       "    time     object 2018-11-05 12:00:00\n",
-       "Attributes:\n",
-       "    short_name:     U\n",
-       "    long_name:      sea_ice_x_velocity\n",
-       "    standard_name:  sea_ice_x_velocity\n",
-       "    units:          cm / s\n",
-       "    valid_min:      -200.0\n",
-       "    valid_max:      200.0\n",
-       "    grid_mapping:   crs\n",
-       "    comment:        U is the along-x component of the ice motion.  It is *not..."
-      ]
-     },
-     "execution_count": 36,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
    "source": [
-    "data = ds.sel(time='2018-11-05').u.squeeze()\n",
-    "data"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 57,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "1"
-      ]
-     },
-     "execution_count": 57,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "ds.x.values.ndim"
+    "I've modified Johan's code to allow coordinate arrays to be 1D or 2D.  \n",
+    "\n",
+    "I've also added some checks to make sure the input arrays are numpy arrays.  This is what the `assert` statements do.  I also check that the sizes of input coordinate and data arrays are compatable."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 146,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "array([-1.9481492, -4.2206554, -4.934302 ], dtype=float32)"
-      ]
-     },
-     "execution_count": 146,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[-1.9481492 -4.2206554 -4.934302 ] [-6.1345673 -5.567279  -4.185838 ]\n"
+     ]
     }
    ],
    "source": [
@@ -725,8 +670,14 @@
     "    assert isinstance(xq, np.ndarray), 'xp must be a numpy array'\n",
     "    assert isinstance(yq, np.ndarray), 'yp must be a numpy array'\n",
     "\n",
-    "    assert x.shape == y.shape == z.shape, \\\n",
-    "           \"input array must have same size and number of dimensions\"\n",
+    "    if (x.ndim == y.ndim == 2):\n",
+    "        assert x.shape == y.shape == z.shape, \\\n",
+    "               \"input array must have same size and number of dimensions\"\n",
+    "    elif (x.ndim == y.ndim == 1):\n",
+    "        assert (*x.shape, *y.shape) == z.shape, \\\n",
+    "               \"input coordinate dimensions must be same size as input data\"\n",
+    "    else:\n",
+    "        raise \n",
     "    assert xq.shape == yq.shape, 'Output coordinate must have same shape'\n",
     "    \n",
     "    # If input coordinates are 2D, get single components\n",
@@ -763,84 +714,37 @@
     "    return zq\n",
     "\n",
     "xx, yy = np.meshgrid(ds.x.values, ds.y.values)\n",
-    "zz = ds.sel(time='2018-11-05').u.squeeze().values\n",
-    "#zz = np.where(~np.isnan(zz), zz, -99.9)\n",
-    "uj = interp2d(xx, yy, zz, np.array(x), np.array(y), order=1)\n",
-    "uj\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 138,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "[ 2  1  0 -1 -2]\n",
-      "[[16 17 18 19]\n",
-      " [12 13 14 15]\n",
-      " [ 8  9 10 11]\n",
-      " [ 4  5  6  7]\n",
-      " [ 0  1  2  3]]\n",
-      "[-2 -1  0  1  2]\n",
-      "[[ 0  1  2  3]\n",
-      " [ 4  5  6  7]\n",
-      " [ 8  9 10 11]\n",
-      " [12 13 14 15]\n",
-      " [16 17 18 19]]\n"
-     ]
-    }
-   ],
-   "source": [
-    "reverse = True\n",
-    "x1 = np.arange(3)-1\n",
-    "y1 = np.arange(5)-2\n",
-    "\n",
-    "x2, y2 = np.meshgrid(x1, y1)\n",
-    "data = np.arange(4*5).reshape(5,4)\n",
-    "\n",
-    "if reverse:\n",
-    "    y1 = np.flip(y1)\n",
-    "    data = np.flipud(data)\n",
-    "    \n",
-    "print (y1)\n",
-    "print (data)\n",
-    "\n",
-    "if y1[-1] < y1[0]:\n",
-    "        y1 = np.flip(y1)\n",
-    "        data = np.flipud(data)\n",
-    "        \n",
-    "print (y1)\n",
-    "print (data)\n",
-    "\n"
+    "zu = ds.sel(time='2018-11-05').u.squeeze().values\n",
+    "zv = ds.sel(time='2018-11-05').v.squeeze().values\n",
+    "uj = interp2d(xx, yy, zu, np.array(x), np.array(y), order=1)\n",
+    "vj = interp2d(xx, yy, zv, np.array(x), np.array(y), order=1)\n",
+    "print (uj, vj)\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "My approach"
+    "## My approach"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [],
    "source": [
     "u = ds.sel(time='2018-11-05').u.values\n",
     "u[np.isnan(u)] = -99.9\n",
     "v = ds.sel(time='2018-11-05').v.values\n",
     "v[np.isnan(v)] = -99.9\n",
-    "fu = interp2d(ds.x, ds.y, u, kind='linear')\n",
-    "fv = interp2d(ds.x, ds.y, v, kind='linear')"
+    "fu = scinterp.interp2d(ds.x, ds.y, u, kind='linear')\n",
+    "fv = scinterp.interp2d(ds.x, ds.y, v, kind='linear')"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 21,
    "metadata": {},
    "outputs": [
     {
@@ -853,7 +757,7 @@
     }
    ],
    "source": [
-    "print ('u = ', np.diag(fu(x, y)))\n",
+    "print ('u = ', np.diag(fu(x, y)))  # This returns a 2D array so use np.diag to get the unique values\n",
     "print ('v = ', np.diag(fv(x, y)))"
    ]
   },
@@ -866,7 +770,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -877,8 +781,8 @@
     "    um = np.where(~np.isnan(u), u, -99.9)\n",
     "    vm = np.where(~np.isnan(v), v, -99.9)\n",
     "    \n",
-    "    fu = interp2d(x, y, um, kind='linear')\n",
-    "    fv = interp2d(x, y, vm, kind='linear')\n",
+    "    fu = scinterp.interp2d(x, y, um, kind='linear')\n",
+    "    fv = scinterp.interp2d(x, y, vm, kind='linear')\n",
     "    \n",
     "    ui = np.array([fu(x_i, y_i) for x_i, y_i in zip(xi, yi)])\n",
     "    vi = np.array([fv(x_i, y_i) for x_i, y_i in zip(xi, yi)])\n",
@@ -896,9 +800,18 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 23,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[-1.94814917 -4.22065543 -4.93430168]\n",
+      "[-6.13456718 -5.56727877 -4.18583824]\n"
+     ]
+    }
+   ],
    "source": [
     "ui, vi = interp_uv(ds.x, ds.y, \n",
     "                   ds.sel(time='2018-11-05').u.values, \n",