ocetrac_flowchart.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "cultural-coaching",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import s3fs\n",
    "import xarray as xr\n",
    "import scipy.ndimage\n",
    "from skimage.measure import regionprops \n",
    "from skimage.measure import label as label_np\n",
    "import dask.array as dsa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
   "metadata": {},
   "outputs": [],
   "source": [
    "/burg/abernathey/users/hillary/oisst_v21"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are on Pangeo Cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
=======
   "id": "olympic-radical",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "endpoint_url = 'https://ncsa.osn.xsede.org'\n",
    "fs_osn = s3fs.S3FileSystem(anon=True, client_kwargs={'endpoint_url': endpoint_url},)\n",
    "\n",
    "path = \"Pangeo/pangeo-forge/noaa_oisst/v2.1-avhrr.zarr\"\n",
    "ds = xr.open_zarr(fs_osn.get_mapper(path), consolidated=True)\n",
    "print(ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "vanilla-greek",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a mask\n",
    "mask = ds.sst.isel(time=0, zlev=0).fillna(-999)\n",
    "mask = mask.where(mask==-999, other = 1)\n",
    "mask = mask.where(mask==1, other = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "rough-brush",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,12))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "mask.plot(transform=ccrs.PlateCarree(), vmin=0, vmax=1, cmap='Greys_r', extend='max', add_colorbar=False, label=False)\n",
    "ax.set_title('')\n",
    "ax.background_patch.set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "universal-bracket",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "sst_anom = ds.anom.sel(time=('2014-03-01'), zlev=0).squeeze()\n",
    "# sst_anom = sst_anom.where(sst_anom >=0.5, other=np.nan)\n",
    "plt.figure(figsize=(12,12))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "sst_anom.plot(transform=ccrs.PlateCarree(), vmin=-3, vmax=3, cmap='RdBu_r', extend='both', add_colorbar=False, label=False)\n",
    "ax.set_title('')\n",
    "ax.coastlines(resolution='110m', color='black', linewidth=1)\n",
    "ax.add_feature(cfeature.LAND, facecolor='w')\n",
    "ax.background_patch.set_visible(False)\n",
    "# ax.outline_patch.set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "varying-focus",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "sst_hot = ds.anom.sel(time=('2014-03-01'), zlev=0).squeeze()\n",
    "sst_hot = sst_hot.where(sst_hot >=0.5, other=np.nan)\n",
    "\n",
    "plt.figure(figsize=(12,12))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "sst_hot.plot(transform=ccrs.PlateCarree(), vmin=-3, vmax=3, cmap='RdBu_r', extend='max', add_colorbar=False, label=False)\n",
    "ax.set_title('')\n",
    "ax.coastlines(resolution='110m', color='black', linewidth=1)\n",
    "ax.add_feature(cfeature.LAND, facecolor='w')\n",
    "ax.background_patch.set_visible(False)\n",
    "# ax.outline_patch.set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "incorporate-moses",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to Binary\n",
    "bitmap_binary = sst_hot.where(sst_hot>0, drop=False, other=0)\n",
    "bitmap_binary = bitmap_binary.where(bitmap_binary==0, drop=False, other=1)\n",
    "\n",
    "plt.figure(figsize=(12,12))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "bitmap_binary.plot(transform=ccrs.PlateCarree(), vmin=0, vmax=1, cmap='Greys', extend='max', add_colorbar=False, label=False)\n",
    "ax.set_title('')\n",
    "ax.coastlines(color='k', linewidth=.5)\n",
    "# ax.add_feature(cfeature.LAND, facecolor='w')\n",
    "ax.background_patch.set_visible(False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "institutional-interface",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define structuring element\n",
    "radius =8\n",
    "diameter = radius*2\n",
    "x = np.arange(-radius, radius+1)\n",
    "x, y = np.meshgrid(x, x)\n",
    "r = x**2+y**2 \n",
    "se = r<radius**2\n",
    "\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "fig, ax = plt.subplots(1, 1,figsize=(6,6))\n",
    "plt.pcolormesh(se, cmap='Greys', alpha=0.5)\n",
    "ax.grid(True, which='both', axis='both', linestyle='-', color='k')\n",
    "ax.set_xticks(np.arange(0,16,1), minor=True)\n",
    "ax.set_yticks(np.arange(0,16,1), minor=True)\n",
    "plt.xlim(0, 17)\n",
    "plt.ylim(0, 17)\n",
    "# plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.axis('square');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "continued-drama",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "def binary_open_close(bitmap_binary):\n",
    "    bitmap_binary_padded = np.pad(bitmap_binary,\n",
    "                                  ((diameter, diameter), (diameter, diameter)),\n",
    "                                  mode='wrap')\n",
    "    s1 = scipy.ndimage.binary_closing(bitmap_binary_padded, se, iterations=1)\n",
    "    s2 = scipy.ndimage.binary_opening(s1, se, iterations=1)\n",
    "    unpadded= s2[diameter:-diameter, diameter:-diameter]\n",
    "    return unpadded\n",
    "\n",
    "\n",
    "\n",
    "sst_hot = ds.anom.sel(zlev=0).squeeze()\n",
    "sst_hot = sst_hot.where(sst_hot >=0.5, other=np.nan)\n",
    "bitmap_binary = sst_hot.where(sst_hot>0, drop=False, other=0)\n",
    "bitmap_binary = bitmap_binary.where(bitmap_binary==0, drop=False, other=1)\n",
    "\n",
    "if bitmap_binary.chunks:\n",
    "#     bitmap_binary = bitmap_binary.chunk({'time': -1})\n",
    "    bitmap_binary = bitmap_binary.chunk({'time': bitmap_binary.sizes['time']})\n",
    "    \n",
    "mo_binary = xr.apply_ufunc(binary_open_close, bitmap_binary,\n",
    "                           input_core_dims=[['time','lat', 'lon']],\n",
    "                           output_core_dims=[['time','lat', 'lon']],\n",
    "                           output_dtypes=[bitmap_binary.dtype],\n",
    "                           vectorize=True,\n",
    "                           dask='parallelized')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "modern-cornell",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "def _apply_mask(binary_images, mask):\n",
    "    binary_images_with_mask = binary_images.where(mask==1, drop=False, other=0)\n",
    "    return binary_images_with_mask\n",
    "\n",
    "mo_binary = _apply_mask(mo_binary, mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "liquid-password",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,12))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "mo_binary.sel(time=('2014-03-01')).plot(transform=ccrs.PlateCarree(), vmin=0, vmax=1, cmap='Greys', extend='max', add_colorbar=False, label=False)\n",
    "ax.set_title('')\n",
    "ax.coastlines(color='k', linewidth=.5)\n",
    "# ax.add_feature(cfeature.LAND, facecolor='w')\n",
    "ax.background_patch.set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "advisory-unknown",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "# label\n",
    "def _label_either(data, **kwargs):\n",
    "        if isinstance(data, dsa.Array):\n",
    "            try:\n",
    "                from dask_image.ndmeasure import label as label_dask\n",
    "                def label_func(a, **kwargs):\n",
    "                    ids, num = label_dask(a, **kwargs)\n",
    "                    return ids\n",
    "            except ImportError:\n",
    "                raise ImportError(\n",
    "                    \"Dask_image is required to use this function on Dask arrays. \"\n",
    "                    \"Either install dask_image or else call .load() on your data.\"\n",
    "                )\n",
    "        else:\n",
    "            label_func = label_np\n",
    "        return label_func(data, **kwargs)\n",
    "    \n",
    "def get_labels(binary_images):\n",
    "        blobs_labels = _label_either(binary_images, background=0)\n",
    "        return blobs_labels\n",
    "    \n",
    "blobs_labels = get_labels(mo_binary)\n",
    "\n",
    "labels = xr.DataArray(blobs_labels, dims=mo_binary.dims, coords=mo_binary.coords)\n",
    "labels = labels.where(labels>0, drop=False, other=np.nan)  \n",
    "\n",
    "# # The labels are repeated each time step, therefore we relabel them to be consecutive\n",
    "# for i in range(1, labels.shape[0]):\n",
    "#     labels[i,:,:] = labels[i,:,:].values + labels[i-1,:,:].max().values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "arctic-michael",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap\n",
    "maxl = int(np.nanmax(labels.values))\n",
    "cm = ListedColormap(np.random.random(size=(maxl, 3)).tolist())\n",
    "\n",
    "plt.figure(figsize=(12,12))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "labels.sel(time=('2014-03-01')).plot(cmap=cm, transform=ccrs.PlateCarree(), add_colorbar=False, label=False)\n",
    "ax.set_title('')\n",
    "ax.coastlines(color='k', linewidth=.5)\n",
    "# ax.add_feature(cfeature.LAND, facecolor='w')\n",
    "ax.background_patch.set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "diverse-blair",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "# wrap labels\n",
    "def _wrap(labels):\n",
    "        ''' Impose periodic boundary and wrap labels'''\n",
    "        \n",
    "        first_column = labels[..., 0]\n",
    "        last_column = labels[..., -1]\n",
    "\n",
    "        unique_first = np.unique(first_column[first_column>0])\n",
    "\n",
    "        # This loop iterates over the unique values in the first column, finds the location of those values in \n",
    "        # the first columnm and then uses that index to replace the values in the last column with the first column value\n",
    "        for i in enumerate(unique_first):\n",
    "            print(i)\n",
    "            first = np.where(first_column == i[1])\n",
    "            last = last_column[first[0]]#, first[1]] # this works only if based on the assumption that the data contains multiple times\n",
    "            bad_labels = np.unique(last[last>0])\n",
    "            replace = np.isin(labels, bad_labels)\n",
    "            labels[replace] = i[1]\n",
    "        \n",
    "#         labels_wrapped = np.unique(labels, return_inverse=True)[1].reshape(labels.shape)\n",
    "\n",
    "        # recalculate the total number of labels \n",
    "        N = np.max(labels)\n",
    "\n",
    "        return labels, N\n",
    "\n",
    "labels_wrapped, N = _wrap(np.array(labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "american-being",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "labels_xr = xr.DataArray(labels_wrapped, dims=labels.dims, coords=labels.coords)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "incorrect-nelson",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "# from matplotlib.colors import ListedColormap\n",
    "# maxl = int(np.nanmax(labels_wrapped))\n",
    "# cm = ListedColormap(np.random.random(size=(maxl, 3)).tolist())\n",
    "\n",
    "plt.figure(figsize=(12,12))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "labels_xr.plot(cmap=cm, transform=ccrs.PlateCarree(), add_colorbar=False, label=False)\n",
    "ax.set_title('')\n",
    "ax.coastlines(color='k', linewidth=.5)\n",
    "# ax.add_feature(cfeature.LAND, facecolor='w')\n",
    "ax.background_patch.set_visible(False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "smaller-railway",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "min_size_quartile = 0.75\n",
    "\n",
    "# Calculate Area of each object and keep objects larger than threshold\n",
    "props = regionprops(labels_xr.astype('int'))\n",
    "labelprops = [p.label for p in props]\n",
    "labelprops = xr.DataArray(labelprops, dims=['label'], coords={'label': labelprops}) \n",
    "area = xr.DataArray([p.area for p in props], dims=['label'], coords={'label': labelprops})  # Number of pixels of the region.\n",
    "min_area = np.percentile(area, min_size_quartile*100)\n",
    "print('minimum area: ', min_area) \n",
    "keep_labels = labelprops.where(area>=min_area, drop=True)\n",
    "keep_where = np.isin(labels_xr, keep_labels)\n",
    "out_labels = xr.DataArray(np.where(keep_where==False, np.nan, labels_xr), dims=mo_binary.dims, coords=mo_binary.coords)\n",
    "\n",
    "# Convert images to binary. All positive values == 1, otherwise == 0\n",
    "binary_labels = out_labels.where(out_labels==0, drop=False, other=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "overhead-relaxation",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,12))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "out_labels.plot(cmap=cm, transform=ccrs.PlateCarree(), add_colorbar=False, label=False)\n",
    "ax.set_title('')\n",
    "ax.coastlines(color='k', linewidth=.5)\n",
    "# ax.add_feature(cfeature.LAND, facecolor='w')\n",
    "ax.background_patch.set_visible(False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
<<<<<<< HEAD
=======
   "id": "typical-glenn",
>>>>>>> 696c5c17b9519cfca3ea335833f5d2c46a7c2213
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ocetrac-scipy-2021]",
   "language": "python",
   "name": "conda-env-ocetrac-scipy-2021-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}