DepthOfUpwellingArianePlan.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Depth of upwelling Ariane experiment \n",
    "\n",
    "Planning for where and when will particles be released. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import matplotlib.colors as mcolors\n",
    "\n",
    "from math import *\n",
    "\n",
    "import scipy.io\n",
    "\n",
    "import scipy as spy\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "from netCDF4 import Dataset\n",
    "\n",
    "import pylab as pl\n",
    "\n",
    "import os \n",
    "\n",
    "import sys\n",
    "\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lib_path = os.path.abspath('../../Building_canyon/BuildCanyon/PythonModulesMITgcm') # Add absolute path to my python scripts\n",
    "sys.path.append(lib_path)\n",
    "\n",
    "import ReadOutTools_MITgcm as rout \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sns.set()\n",
    "sns.set_style('white')\n",
    "sns.set_context('paper')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "CanyonGrid='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run03/gridGlob.nc'\n",
    "CanyonGridOut = Dataset(CanyonGrid)\n",
    "#for dimobj in CanyonGridOut.variables.values():\n",
    "#    print dimobj\n",
    "CanyonState='/ocean/kramosmu/MITgcm/TracerExperiments/CNTDIFF/run03/stateGlob.nc'\n",
    "CanyonStateOut = Dataset(CanyonState)\n",
    "#for dimobj in CanyonStateOut.variables.values():\n",
    "#    print dimobj\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'HFacS'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-5-9e2d9bd71fed>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     24\u001b[0m \u001b[0mMaskW\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrout\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetMask\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mCanyonGrid\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'HFacW'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     25\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 26\u001b[1;33m \u001b[0mhFacS\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrout\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetField\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mCanyonGrid\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'HFacS'\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# (Z,Yp1,X) for meridional fluxes\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     27\u001b[0m \u001b[0mMaskS\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrout\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetMask\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mCanyonGrid\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'HFacS'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/ocean/kramosmu/Building_canyon/BuildCanyon/PythonModulesMITgcm/ReadOutTools_MITgcm.pyc\u001b[0m in \u001b[0;36mgetField\u001b[1;34m(statefile, fieldname)\u001b[0m\n\u001b[0;32m     20\u001b[0m     \u001b[0mStateOut\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mDataset\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstatefile\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     21\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 22\u001b[1;33m     \u001b[0mFld\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mStateOut\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvariables\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mfieldname\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     23\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     24\u001b[0m     \u001b[0mshFld\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mFld\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyError\u001b[0m: 'HFacS'"
     ]
    }
   ],
   "source": [
    "# General input\n",
    "\n",
    "nx = 360\n",
    "ny = 360\n",
    "nz = 90\n",
    "\n",
    "nt = 19 # t dimension size \n",
    "\n",
    "Z = CanyonStateOut.variables['Z']\n",
    "#print(z[10])\n",
    "Time = CanyonStateOut.variables['T']\n",
    "#print(Time[:])\n",
    "\n",
    "xc = rout.getField(CanyonGrid, 'XC') # x coords tracer cells\n",
    "rc = CanyonGridOut.variables['RC']\n",
    "yc = rout.getField(CanyonGrid, 'YC') # y coords tracer cells\n",
    "\n",
    "bathy = rout.getField(CanyonGrid, 'Depth')\n",
    "\n",
    "hFacC = rout.getField(CanyonGrid, 'HFacC')\n",
    "MaskC = rout.getMask(CanyonGrid, 'HFacC')\n",
    "\n",
    "hFacW = rout.getField(CanyonGrid, 'HFacW') # (Z,Y,Xp1) for zonal fluxes\n",
    "MaskW = rout.getMask(CanyonGrid,'HFacW') \n",
    "\n",
    "hFacS = rout.getField(CanyonGrid, 'HFacS') # (Z,Yp1,X) for meridional fluxes\n",
    "MaskS = rout.getMask(CanyonGrid,'HFacS') \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6,6))\n",
    "\n",
    "CS = plt.plot(range(360),-bathy[:,120],'ko-',markersize=5,linewidth=0.5)\n",
    "plt.plot(226,-bathy[226,120],'ro',markersize=8)\n",
    "plt.axhline(y=Z[29],xmin=0,xmax=360,linestyle='-', color='0.75')\n",
    "\n",
    "plt.ylabel('Depth (m)')\n",
    "plt.xlabel('y index')\n",
    "print('The shelf-break isobath (sb=%.2f m) upstream of the canyon corresponds to the y-index 198' % bathy[198,120])\n",
    "print('The corresponding index for z is   (z(29)=%.2f m)' % -Z[29])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.set(context='paper', style='white', font='sans-serif', font_scale=2)\n",
    "\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "\n",
    "#plt.plot(xc[150:300:5,15],yc[150:300:5,15], 'o',color = 'orange')\n",
    "#plt.plot(xc[150:267:10,45],yc[150:267:10,45], 'o',color = 'orange')\n",
    "#plt.plot(xc[150:300:5,80],yc[150:300:5,80], 'o',color = 'orange')\n",
    "plt.plot(xc[150:267:10,150],yc[150:267:10,150], 'o',color = 'orange')\n",
    "\n",
    "#for ii in range(0,360,10):\n",
    "#    plt.axhline(y=yc[ii,1],xmin=0,xmax=360,linestyle='-', color='0.75')\n",
    "\n",
    "levels = [20,50,100,148.17,200,300,400,500,600,700,800,900,1000,1100,1200]\n",
    "\n",
    "cmapPurple = sns.cubehelix_palette(light=1, as_cmap=True)\n",
    "cmapBlue = sns.cubehelix_palette(start=.5, rot=-.75, as_cmap=True)\n",
    "\n",
    "CS = plt.contourf(xc,yc,bathy,levels,cmap=cmapPurple)\n",
    "                  \n",
    "                  #linewidths=[0.75,0.75,0.75,2,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75] )\n",
    "plt.colorbar(CS)\n",
    "#plt.clabel(CS, \n",
    "#          inline=1,\n",
    "#          fmt='%1.1f',\n",
    "#           fontsize=12)\n",
    "\n",
    "    \n",
    "plt.xticks(np.arange(0,120000,10000), [' ','10','20','30','40','50','60','70','80','90','100','110','120'])\n",
    "plt.yticks(np.arange(0,90000,10000), [' ','10','20','30','40','50','60','70','80','90'])\n",
    "\n",
    "plt.xlabel('km',fontsize=14)\n",
    "plt.ylabel('km',fontsize=14)\n",
    "\n",
    "plt.title('Depth (m)',fontsize=16)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.set(context='paper', style='white', font='sans-serif', font_scale=2)\n",
    "\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "\n",
    "y,z = np.meshgrid(yc[:,150],rc[:])\n",
    "\n",
    "y=np.ma.masked_array(y,mask=MaskC[:,:,150])\n",
    "z=np.ma.masked_array(z,mask=MaskC[:,:,150])\n",
    "\n",
    "plt.plot(y[15:56:1,150:267:10],z[15:56:1,150:267:10], 'o',color = 'orange')\n",
    "\n",
    "plt.plot(yc[:,150],-bathy[:,150],'k-',linewidth=1)\n",
    "\n",
    "plt.axhline(-30,color='0.50',linestyle='--')\n",
    "plt.axhline(-500,color='0.50',linestyle='--')\n",
    "#plt.xticks(np.arange(0,120000,10000), [' ','10','20','30','40','50','60','70','80','90','100','110','120'])\n",
    "#plt.yticks(np.arange(0,90000,10000), [' ','10','20','30','40','50','60','70','80','90'])\n",
    "\n",
    "plt.xlabel('km',fontsize=14)\n",
    "plt.ylabel('m',fontsize=14)\n",
    "\n",
    "plt.title('Cross-shelf slice at upstream side of canyon',fontsize=16)\n",
    "\n",
    "print(np.shape(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zoom in plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.set(context='paper', style='white', font='sans-serif', font_scale=2)\n",
    "\n",
    "fig,(ax1,ax2) = plt.subplots(1,2 ,figsize=(15,6))\n",
    "\n",
    "#ax1.plot(xc[150:267:10,15],yc[150:267:10,15], 'o',color = 'red')\n",
    "#ax1.plot(xc[150:267:10,45],yc[150:267:10,45], 'o',color = 'red')\n",
    "#ax1.plot(xc[150:267:10,80],yc[150:267:10,80], 'o',color = 'red')\n",
    "ax1.plot(xc[150:267:10,150],yc[150:267:10,150], 'o',color = 'red')\n",
    "\n",
    "#for ii in range(0,360,10):\n",
    "#    plt.axhline(y=yc[ii,1],xmin=0,xmax=360,linestyle='-', color='0.75')\n",
    "\n",
    "levels = [20,50,100,148.17,200,267,400,500,600,700,800,900,1000,1100,1200]\n",
    "\n",
    "cmapPurple = sns.cubehelix_palette(light=1, as_cmap=True)\n",
    "cmapBlue = sns.cubehelix_palette(start=.5, rot=-.75, as_cmap=True)\n",
    "\n",
    "CS = ax1.contourf(xc[148:301,14:240],yc[148:301,14:240],bathy[148:301,14:240],levels,cmap=cmapPurple)\n",
    "                  \n",
    "plt.colorbar(CS,ax=ax1)\n",
    "    \n",
    "ax1.set_xticks(np.arange(10000,70000,10000), ['10','20','30','40','50','60','70'])\n",
    "ax1.set_yticks(np.arange(40000,70000,10000), ['40','50','60','70'])\n",
    "\n",
    "ax1.set_xlabel('km',fontsize=14)\n",
    "ax1.set_ylabel('km',fontsize=14)\n",
    "\n",
    "ax1.set_title('Depth (m)',fontsize=16)\n",
    "\n",
    "# Ax2\n",
    "\n",
    "plt.plot(y[15:35:2,150:267:10],z[15:35:2,150:267:10], 'o',color = 'orange')\n",
    "plt.plot(y[35:56:1,150:267:10],z[35:56:1,150:267:10], 'o',color = 'red')\n",
    "\n",
    "plt.plot(yc[148:301,150],-bathy[148:301,150],'k-',linewidth=2)\n",
    "\n",
    "\n",
    "plt.axhline(-30,color='0.50',linestyle='--')\n",
    "plt.axhline(-500,color='0.50',linestyle='--')\n",
    "#plt.xticks(np.arange(0,120000,10000), [' ','10','20','30','40','50','60','70','80','90','100','110','120'])\n",
    "#plt.yticks(np.arange(0,90000,10000), [' ','10','20','30','40','50','60','70','80','90'])\n",
    "\n",
    "plt.xlabel('km',fontsize=14)\n",
    "plt.ylabel('m',fontsize=14)\n",
    "\n",
    "plt.title('Cross-shelf slice at x=10km form west bdy',fontsize=16)\n",
    "\n",
    "print('Each red dot is a release site')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Horizontal and vertical spacing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.set(context='paper', style='white', font='sans-serif', font_scale=2)\n",
    "\n",
    "fig,(ax1,ax2) = plt.subplots(1,2 ,figsize=(15,6))\n",
    "\n",
    "ax1.plot(yc[150:262:10,150]/1000.0,yc[150:262:10,150]-yc[155:267:10,150], 'o',color = 'red')\n",
    "\n",
    "\n",
    "ax1.set_xlabel('Cross-shelf position (km from southern bdy) ',fontsize=14)\n",
    "ax1.set_ylabel('Spacing (m)',fontsize=14)\n",
    "\n",
    "ax1.set_title('Horizontal spacing',fontsize=16)\n",
    "\n",
    "# Ax2\n",
    "\n",
    "plt.plot(rc[6:32:2],rc[6:32:2]-rc[8:34:2], 'o',color = 'orange')\n",
    "plt.plot(rc[34:55],rc[34:55]-rc[35:56], 'o',color = 'red')\n",
    "\n",
    "\n",
    "plt.xlabel(' Depth starting position (m)',fontsize=14)\n",
    "plt.ylabel('Vertical spacing (m from surface)',fontsize=14)\n",
    "\n",
    "plt.title('Cross-shelf slice at x=10km form west bdy',fontsize=16)\n",
    "\n",
    "#print('Each red dot is a release site')\n",
    "print('Horizontal spacing 550 m near shelf break and grows to up to 1400 m at inner shelf ')\n",
    "print('Vertical spacing 10 m in the first 250 m and 20 meters from 250-500m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare initial_positions.txt file\n",
    "5 columns:\n",
    "\n",
    "| 1       | 2       |  3                                      | 4          | 5\n",
    "|:-------:|:-------:|:---------------------------------------:|:----------:|:-----------------------------------:\n",
    "| x index | y index | z index (negative to stay in one depth)   | time index | always 1 for qualitative    \n",
    "\n",
    "\n",
    "The idea is to create a dataframe with these 5 columns and as many rows as particles I want to include and then save as a csv file with delimiter set to 'tab'.\n",
    "\n",
    "Also, remember indices start at 1 in Ariane, not in 0 as in python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Columns 1,2,3\n",
    "\n",
    "#yy = np.arange(0,360)\n",
    "#xx = np.arange(0,360)\n",
    "#zi = np.arange(0,90)\n",
    "\n",
    "#xi,yi = np.meshgrid(xx,yy)  # These are like yc and xc but index version\n",
    "\n",
    "#xind = [150] # 45,110 x index for lines along y \n",
    "\n",
    "#times = [10] # 3,10 time indices\n",
    "\n",
    "#for time in times:\n",
    "    \n",
    "#    for ii in range(len(xind)):\n",
    "    \n",
    "#        ypmesh,zpmesh = np.meshgrid(yi[:,xind[ii]],zi[:])\n",
    "\n",
    "#        ypmask = np.ma.masked_array(ypmesh,mask=MaskC[:,:,xind[ii]])\n",
    "#        zpmask = np.ma.masked_array(zpmesh,mask=MaskC[:,:,xind[ii]])\n",
    "\n",
    "        # Select every 2 cells for the first 35 z-levles and then every 1; select every 5 cells along y-direction\n",
    "        \n",
    "#        ypslice1 = ypmask[15:35:2,150:267:10] # before 15:35:2 it was 6:35:2 \n",
    "#        ypslice2 = ypmask[35:56:1,150:267:10]\n",
    "#        yp = np.ma.append(ypslice1,ypslice2,axis=0)\n",
    "\n",
    "#        zpslice1 = zpmask[15:35:2,150:267:10] # before 15:35:2 it was 6:35:2 \n",
    "#        zpslice2 = zpmask[35:56:1,150:267:10]\n",
    "#        zp = np.ma.append(zpslice1,zpslice2,axis=0)\n",
    "\n",
    "#        yind = np.ma.reshape(yp,(1,np.shape(yp)[0]*np.shape(yp)[1]))\n",
    "#        zind = np.ma.reshape(zp,(1,np.shape(zp)[0]*np.shape(zp)[1]))\n",
    "#        xindcol = xind[ii]*np.ones(np.shape(yind))\n",
    "\n",
    "#        if ii == 0:\n",
    "#            col1part = xindcol\n",
    "#            col2part = yind\n",
    "#            col3part = zind\n",
    "             \n",
    "    \n",
    "#        else:\n",
    "#            col1part = np.ma.append(col1part,xindcol)\n",
    "#            col2part = np.ma.append(col2part,yind)\n",
    "#            col3part = np.ma.append(col3part,zind)\n",
    "    \n",
    "#    if time == times[0]:\n",
    "#        col1 = col1part\n",
    "#        col2 = col2part\n",
    "#        col3 = col3part\n",
    "#        col4 = time*np.ones(len(col1part))\n",
    "#        col5 = np.ones(len(col1part))\n",
    "    \n",
    "#    else:\n",
    "#        col1 = np.ma.append(col1,col1part)\n",
    "#        col2 = np.ma.append(col2,col2part)\n",
    "#        col3 = np.ma.append(col3,col3part)\n",
    "#        col4 = np.append(col4,time*np.ones(len(col1part)))\n",
    "#        col5 = np.append(col5,np.ones(len(col1part)))\n",
    "        \n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "#Use this cell if there is only one slice in x\n",
    "# Columns 1,2,3\n",
    "\n",
    "yy = np.arange(0,360)\n",
    "xx = np.arange(0,360)\n",
    "zi = np.arange(0,90)\n",
    "\n",
    "xi,yi = np.meshgrid(xx,yy)  # These are like yc and xc but index version\n",
    "\n",
    "xind = 150# 45,110 x index for lines along y \n",
    "\n",
    "times = 10 # 3,10 time indices\n",
    "\n",
    "time = 0.5\n",
    "    \n",
    "ii=0\n",
    "    \n",
    "ypmesh,zpmesh = np.meshgrid(yi[:,xind],zi[:])\n",
    "\n",
    "ypmask = np.ma.masked_array(ypmesh,mask=MaskC[:,:,xind])\n",
    "zpmask = np.ma.masked_array(zpmesh,mask=MaskC[:,:,xind])\n",
    "\n",
    "# Select every 2 cells for the first 35 z-levles and then every 1; select every 5 cells along y-direction\n",
    "        \n",
    "ypslice1 = ypmask[15:35:2,150:267:10] # before 15:35:2 it was 6:35:2 \n",
    "ypslice2 = ypmask[35:56:1,150:267:10]\n",
    "yp = np.ma.append(ypslice1,ypslice2,axis=0)\n",
    "\n",
    "zpslice1 = zpmask[15:35:2,150:267:10] # before 15:35:2 it was 6:35:2 \n",
    "zpslice2 = zpmask[35:56:1,150:267:10]\n",
    "zp = np.ma.append(zpslice1,zpslice2,axis=0)\n",
    "\n",
    "yind = np.ma.reshape(yp,(np.shape(yp)[0]*np.shape(yp)[1]))\n",
    "zind = np.ma.reshape(zp,(np.shape(zp)[0]*np.shape(zp)[1]))\n",
    "xindcol = xind*np.ones(np.shape(yind))\n",
    "\n",
    "\n",
    "col1part = xindcol\n",
    "col2part = yind\n",
    "col3part = zind\n",
    "             \n",
    "    \n",
    "    \n",
    "\n",
    "col1 = col1part\n",
    "col2 = col2part\n",
    "col3 = col3part\n",
    "col4 = time*np.ones(np.shape(col1part))\n",
    "col5 = np.ones(np.shape(col1part))\n",
    "    \n",
    "       \n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.shape(col4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "raw_data = { 'nx': col1, 'ny':col2, 'nz':col3, 'nt':col4,'param':col5}\n",
    "df = pd.DataFrame(raw_data, columns = ['nx','ny','nz','nt','param'])\n",
    "\n",
    "#remove rows with NaN's. NaN's correspond to masked values.\n",
    "df = df[np.isfinite(df['nz'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(len(col1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check that I am choosing the right positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.set(context='paper', style='white', font='sans-serif', font_scale=2)\n",
    "\n",
    "fig,(ax2) = plt.subplots(1,1 ,figsize=(10,6))\n",
    "\n",
    "YY=np.array(df.ny) \n",
    "ZZ=np.array(df.nz)\n",
    "# Ax2\n",
    "for jj,kk in zip(YY[:600],ZZ[:600]):\n",
    "    \n",
    "    plt.plot(yc[jj,150],z[kk,jj], 'o',color = 'orange')\n",
    "\n",
    "    plt.plot(yc[148:301,150],-bathy[148:301,150],'k-',linewidth=2)\n",
    "     \n",
    "    plt.plot(yc[260,150],z[24,260],'ko' )\n",
    "\n",
    "    plt.axhline(-30,color='0.50',linestyle='--')\n",
    "    plt.axhline(-500,color='0.50',linestyle='--')\n",
    "    #plt.xticks(np.arange(0,120000,10000), [' ','10','20','30','40','50','60','70','80','90','100','110','120'])\n",
    "    #plt.yticks(np.arange(0,90000,10000), [' ','10','20','30','40','50','60','70','80','90'])\n",
    "\n",
    "    plt.xlabel('km',fontsize=14)\n",
    "    plt.ylabel('m',fontsize=14)\n",
    "\n",
    "    plt.title('Cross-shelf slice at upstream side of canyon',fontsize=16)\n",
    "    \n",
    "    #print('Each red dot is a release site')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# add 1 to all indices x,y,z,t to match fortran indexing\n",
    "df.nx = df.nx+1\n",
    "df.ny = df.ny+1\n",
    "df.nz = df.nz # ensure particles don't touch the bottom\n",
    "df.nt = df.nt+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df = df[df.ny != 261]\n",
    "\n",
    "\n",
    "df.head(20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(len(df.nx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#save data frame\n",
    "filename1 = ('ArianeInitialFiles/UpwellingDepthExp372_upstreamCanyon_tt1p5_particles.csv') \n",
    "df.to_csv(filename1,sep = '\\t', index=False, header=False) # index=false, header=false avoid writing the index column nad header row\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}