prepare_ship_xy_position_from_AIS_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocess AIS Data Set 2 and derieve XY positionings for vessels in Hamburg harbour \n",
    "\n",
    "### Purpose: Convert all geopositions of 6 different AIS ship tracks to cartesian coordinates to make transient (i.e. non-stationary) CFD simulations possible in which all ships move within the computational domain.\n",
    "\n",
    "### Challenge: The AIS tracks are partly longer than the CFD computational domain and each track enters or leaves the domain on different directions and times. Each vessel track characteristics needs to be considered. \n",
    "\n",
    "#### Requirements: <br>- build unified time frame (4 hours envisaged) and store all positions of each vessel in separete x,y positioning files (x & y separately). <br>- positions must be relative to each vessel's start coordinate\n",
    "\n",
    "#### 1st data sighting HGE MO 13 Aug 2018\n",
    "#### Productive version HGE THU 16 Aug 2018\n",
    "\n",
    "no visualisation in this notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import errno\n",
    "import csv\n",
    "import numpy as np\n",
    "import math\n",
    "import datetime\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "#import folium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Set the order\n",
    "\n",
    "### Aligns from the beginning of a vessel movement or 10:00:\n",
    "\n",
    "ID 1\n",
    "- 10:00:00 (- 11:08:30)\n",
    "\n",
    "ID 2\n",
    "- 10:00:00 (- 13:17:30)\n",
    "\n",
    "ID 3\n",
    "- 10:28:00 (- 11:13:30) \n",
    "\n",
    "ID 4\n",
    "- 10:51:30 (- 13:25:30)\n",
    "\n",
    "ID 5\n",
    "- 12:26:30 (- 13:59:30)\n",
    "\n",
    "ID 6\n",
    "- 13:04:30 (- 13:59:30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vessel_id = ['add_id','add_id','add_id',\\\n",
    "             'add_id','add_id','add_id']\n",
    "\n",
    "inpath = 'add_path_to/csv/'\n",
    "\n",
    "outpath = 'add_path_/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General defintions wrt mapping:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### _Define map projection to establish conversion from cartesian XX,YY grid:_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m = Basemap(projection='tmerc',resolution='i',\n",
    "            llcrnrlat=53., urcrnrlat=54.,\n",
    "            llcrnrlon=8.5, urcrnrlon=10.5,\n",
    "            lat_0=53.5,lon_0=10.) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### _set size of simulation domain (should be same size as region which is analysed):_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sim_domain_edge_length = 10000. # meter\n",
    "sim_domain_half_width  = sim_domain_edge_length/2.\n",
    "\n",
    "# Base coordinate of simulation domain (position of simulation domain centre point):\n",
    "# this is the coord. from which in each direction (N,S,W,E) the domain half-width of \n",
    "# 5km begins. Means, entering or leaving the domain is calculated relative the the \n",
    "# vertical or horizonatalaxis that is intersecting this center point. In this case \n",
    "# about the middle of Blohm+Voss shipyard:\n",
    "lon_sim_base = 9.961250\n",
    "lat_sim_base = 53.539373"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### _Map the geographic base coordinate of simulation domain to global (projection specific) cartesian coordinates:_\n",
    "This point is fixed and vessel independent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y: 98133.64755970874 58975.38045677248\n"
     ]
    }
   ],
   "source": [
    "xpt0_sim_base, ypt0_sim_base = m(lon_sim_base,lat_sim_base)\n",
    "print('Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y:',\\\n",
    "      xpt0_sim_base, ypt0_sim_base)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def vessel_read(filename):\n",
    "    if not os.path.isfile(filename):\n",
    "        raise FileNotFoundError(\n",
    "            errno.ENOENT, os.strerror(errno.ENOENT), filename)\n",
    "    else:\n",
    "        print('read file :',filename)\n",
    "        with open(filename) as csvfile:\n",
    "            reader = csv.DictReader(csvfile)\n",
    "            t_unix= [] # unix time stamp in ms\n",
    "            x_geo = []\n",
    "            y_geo = []\n",
    "            for row in reader:\n",
    "                t_unix.append(row['ts'])\n",
    "                x_geo.append(row['x'])\n",
    "                y_geo.append(row['y'])\n",
    "        return  t_unix, x_geo, y_geo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# just for control time intervalls in incoming data - this adds no operational value\n",
    "def inspect_time_intervals(t):\n",
    "    # loop until n-1 because upwind differences applied:\n",
    "    for idx,obj in enumerate(range(0,len(t)-1)):   \n",
    "        dt = int(t[idx+1]) - int(t[idx])\n",
    "        print('    ',idx,dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# converts date time string back into unix time since epoch in ms:\n",
    "def time_string_to_unixtime(obj):\n",
    "    dt_obj = datetime.datetime.strptime(obj,  '%Y-%m-%d %H:%M:%S') \n",
    "    millisec = dt_obj.timestamp() * 1000\n",
    "    return int(millisec)\n",
    "\n",
    "\n",
    "# converts datetime object (that is something different than above!!!) back \n",
    "# into unix time since epoch in ms:\n",
    "def time_datetime_to_unixtime(obj):\n",
    "    millisec = obj.timestamp() * 1000\n",
    "    return int(millisec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# not all direction filters implemented yet\n",
    "# further modularisation possible (rel. hardwired so far)\n",
    "# Potentially a python class is better suiting the purpose of a filter\n",
    "# Note that each AIS track shows very individual characteristics - hence each filter \n",
    "# applies very different methodologies. \n",
    "# Only track 1 & 4 have enough similarities that they can be processed by the same method.\n",
    "\n",
    "# idea: loop through time \n",
    "#       each dt has a unique x & y position\n",
    "#       compare & calc. difference vessel pos to overall base-pos. (=centre) of simulation \n",
    "#       domain partly look 2min ahead or backward for modifying track positions\n",
    "\n",
    "\n",
    "\n",
    "# ID 6  enters westward (coming from east), hence in negative x direction\n",
    "#       northbound when entering domain (this is not necessarily the overall/average \n",
    "#       direction!)\n",
    "#       Here:\n",
    "#       - the initial vessel postion must NOT be reset its close enough (-4819m off \n",
    "#         sim-base)\n",
    "#       - hence, not requiring recalculation of all dispalcements relative to reset \n",
    "#         initial position\n",
    "#\n",
    "def enter_from_east_northbound(t_unmod, x_pos_unmod, y_pos_unmod,\\\n",
    "                               dx_unmod, dy_unmod, x_ini, y_ini):    # ID 6 (green)\n",
    "    \n",
    "    print('    Track filtering enter_from_east_northbound')\n",
    "        \n",
    "    # loop through reversed X direction \n",
    "    \n",
    "    x_mod_x_scan = x_pos_unmod.copy()   \n",
    "    y_mod_y_scan = y_pos_unmod.copy()  \n",
    "    \n",
    "    dx_mod_x_scan = dx_unmod.copy()      \n",
    "    dy_mod_y_scan = dy_unmod.copy()   \n",
    "    \n",
    "    \n",
    "    # calc difference vessel start pos to overall base-pos. (=centre) of simulation domain:\n",
    "    \n",
    "    dx_vessel_start_to_sim_centre = xpt0_sim_base - x_ini\n",
    "    \n",
    "    dy_vessel_start_to_sim_centre = ypt0_sim_base - y_ini\n",
    "    \n",
    "     \n",
    "    #!!!!!!!!!!!!! DEBUGGING !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n",
    "    print('xpt0_sim_base: ',xpt0_sim_base,'  x_ini: ',x_ini)\n",
    "    print('ypt0_sim_base: ',ypt0_sim_base,'  y_ini: ',y_ini)\n",
    "\n",
    "    print('dx_vessel_start_to_sim_centre: ',dx_vessel_start_to_sim_centre,\\\n",
    "          '  T/F: ',abs(dx_vessel_start_to_sim_centre) > sim_domain_half_width)\n",
    "    print('dy_vessel_start_to_sim_centre: ',dy_vessel_start_to_sim_centre,\\\n",
    "          '  T/F: ',abs(dy_vessel_start_to_sim_centre) > sim_domain_half_width)\n",
    "\n",
    "    \n",
    "    \n",
    "    # is this position outside simulation domain?\n",
    "    # if so, raise error because this shouldn't happen here (scan for leaving domain only)\n",
    "    \n",
    "    if abs(dx_vessel_start_to_sim_centre) > sim_domain_half_width:\n",
    "        raise ValueError('Vessel initial X position outside simulation domain: {}'.format(dx_vessel_start_to_sim_centre))\n",
    "    else:\n",
    "        print('X direction distance of vessel initial position:',dx_vessel_start_to_sim_centre)\n",
    "        \n",
    "    if abs(dy_vessel_start_to_sim_centre) > sim_domain_half_width:\n",
    "        # since in this case its only about 400m I assume its still within the \n",
    "        # computational domain which is slightly bigger (circular) than the simulation or\n",
    "        # analysis domain and we do not raise the error here\n",
    "        print('Vessel initial Y position outside simulation domain: ',\\\n",
    "              dy_vessel_start_to_sim_centre)\n",
    "    else:\n",
    "        print('Vessel initial Y position inside simulation domain:',\\\n",
    "              dy_vessel_start_to_sim_centre)\n",
    "\n",
    "    \n",
    "    \n",
    "    #- X scan\n",
    "    for idt, t in enumerate(t_unmod): \n",
    "\n",
    "        # for ID 6 that way around (positive distance from sim-base) \n",
    "        # which is getting smaller with time & vessel motion (at beginning of track)\n",
    "        \n",
    "        dx_distance_to_sim_base =  x_pos_unmod[idt] - xpt0_sim_base         \n",
    "                 \n",
    "            \n",
    "    #- Y scan\n",
    "    for idt, t in enumerate(t_unmod):\n",
    "        \n",
    "        dy_distance_to_sim_base =  y_pos_unmod[idt] - ypt0_sim_base\n",
    "        \n",
    "   \n",
    "    # correct just the very last displacements (otherwise it jumps back to inital pos):\n",
    "    \n",
    "    dx_mod_x_scan[ len(dx_mod_x_scan)-1 ] = dx_mod_x_scan[ len(dx_mod_x_scan)-2 ]\n",
    "    dy_mod_y_scan[ len(dy_mod_y_scan)-1 ] = dy_mod_y_scan[ len(dy_mod_y_scan)-2 ]\n",
    "    \n",
    "    return x_mod_x_scan, y_mod_y_scan, dx_mod_x_scan, dy_mod_y_scan\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# ID 5: enters eastward (coming from west), hence in positve x direction \n",
    "#       here:\n",
    "#       - the initial vessel postion must be set to a position near simulation domain\n",
    "#       - requires recalculation of all dispalcements relative to reset initial position\n",
    "#       - Y scan is necessary because original AIS vessel starting position's y > y_base+5km\n",
    "#       - the time stamp of the new intial position is determined from the latest occurence\n",
    "#         of entering the sim-domain (of the original track) in either X or Y direction\n",
    "#       -> actually in this case (ID5) in Y it reaches sim-domain in time step 306 whereas \n",
    "#         in X it reaches sim-domain at time step 353 -> hence 353 determines start pos \n",
    "#         and displacements \n",
    "#\n",
    "def enter_from_west_southbound(t_unmod, x_pos_unmod, y_pos_unmod,\\\n",
    "                                  dx_unmod, dy_unmod, x_ini, y_ini):  # ID 5  (light blue)\n",
    "    \n",
    "    print('    Track filtering enter_from_west_southbound')\n",
    "        \n",
    "    # loop in positive X direction \n",
    "    \n",
    "    x_mod_x_scan = x_pos_unmod.copy()    \n",
    "    y_mod_y_scan = y_pos_unmod.copy()   \n",
    "    \n",
    "    dx_mod_x_scan = dx_unmod.copy()     \n",
    "    dy_mod_y_scan = dy_unmod.copy()     \n",
    "    \n",
    "    \n",
    "    # calc difference vessel start pos to overall base-pos. (=centre) of simulation domain:\n",
    "    \n",
    "    dx_vessel_start_to_sim_centre = xpt0_sim_base - x_ini\n",
    "    \n",
    "    dy_vessel_start_to_sim_centre = ypt0_sim_base - y_ini\n",
    "    \n",
    "     \n",
    "    print('xpt0_sim_base: ',xpt0_sim_base,'  x_ini: ',x_ini)\n",
    "    print('ypt0_sim_base: ',ypt0_sim_base,'  y_ini: ',y_ini)\n",
    "\n",
    "    print('dx_vessel_start_to_sim_centre: ',dx_vessel_start_to_sim_centre,\\\n",
    "          '  T/F: ',abs(dx_vessel_start_to_sim_centre) > sim_domain_half_width)\n",
    "    print('dy_vessel_start_to_sim_centre: ',dy_vessel_start_to_sim_centre,\\\n",
    "          '  T/F: ',abs(dy_vessel_start_to_sim_centre) > sim_domain_half_width)\n",
    "\n",
    "    \n",
    "    # is this position outside simulation domain?\n",
    "    # if so, raise error because this shouldn't happen here (we scan for leaving domain only)\n",
    "    \n",
    "    if abs(dx_vessel_start_to_sim_centre) > sim_domain_half_width:\n",
    "        print('Vessel initial X position outside simulation domain: ',\\\n",
    "              dx_vessel_start_to_sim_centre)\n",
    "    else:\n",
    "        print('X direction distance of vessel initial position:',\\\n",
    "              dx_vessel_start_to_sim_centre)\n",
    "        \n",
    "    if abs(dy_vessel_start_to_sim_centre) > sim_domain_half_width:\n",
    "        print('Vessel initial Y position outside simulation domain: ',\\\n",
    "              dy_vessel_start_to_sim_centre,\\\n",
    "              'SPECIAL FILTERING NECESSARY !!!!!!!!!')\n",
    "    else:\n",
    "        print('Vessel initial Y position inside simulation domain:',\\\n",
    "              dy_vessel_start_to_sim_centre)\n",
    "\n",
    "    \n",
    "    # scan for positions ahead:\n",
    "    \n",
    "    t_offset = 4 # for ID 2min are likely 4 timesteps\n",
    "    \n",
    "    \n",
    "    #- X scan\n",
    "    for idt, t in enumerate(t_unmod[:-t_offset]):               \n",
    "        \n",
    "        dt_4_pos_backward = t_unmod[idt + t_offset]\n",
    "        \n",
    "        \n",
    "        \n",
    "        # X modification case ID 5:\n",
    "        # comes from smaller x values going to bigger x values in positive \n",
    "        # direction!\n",
    "        \n",
    "        # define vessel stopping criterium: if vessel reaches sim-domain edge, \n",
    "        # all values up to the end of overall time frame are replicated from \n",
    "        # 2 min backward from the edge position (and loop is interrupted):\n",
    "        \n",
    "        \n",
    "        dx_distance_to_sim_base =  x_pos_unmod[idt] - xpt0_sim_base        \n",
    "         \n",
    "        \n",
    "        if dx_distance_to_sim_base > (-sim_domain_half_width):     \n",
    "            \n",
    "            \n",
    "            # should interrupt at idx 355 - save value of 4 steps backward \n",
    "            # as new intial position: (2nd criterium here is that last time \n",
    "            # step's value is not fullfilling the 1st criterium)\n",
    "            \n",
    "            if (x_pos_unmod[idt-1] - xpt0_sim_base) < (-sim_domain_half_width):\n",
    "                \n",
    "                xpt0_vessel_reset = x_mod_x_scan[idt - t_offset]   \n",
    "                \n",
    "            \n",
    "            # modify range from 2min backward of this position up the \n",
    "            # beginning of list (!!!):\n",
    "        \n",
    "            x_mod_x_scan[:idt - t_offset]  = x_mod_x_scan[idt - t_offset]  \n",
    "            dx_mod_x_scan[:idt - t_offset] = dx_mod_x_scan[idt - t_offset]\n",
    "            \n",
    "            \n",
    "            # keep index with offset (!!) to compare with Y scan\n",
    "            \n",
    "            idt_x = idt - t_offset\n",
    "           \n",
    "             \n",
    "            # to ensure that the position 4 time steps away is not too far away \n",
    "            # (sometimes spacing is not 30sec but up to 120sec or even more), \n",
    "            # calculate the time difference between the two points (in seconds):\n",
    "            \n",
    "            t_difference_this_scan = (t_unmod[idt - t_offset] - t ) /1000.\n",
    "            \n",
    "            print('   delta time is : ',t_difference_this_scan,' sec')\n",
    "            \n",
    "            if abs(t_difference_this_scan) > 120.:\n",
    "                raise ValueError('Vessel stop time diff > 120sec: {}'.format(t_difference_this_scan))\n",
    "    \n",
    "            break\n",
    "        \n",
    "    #- Y scan - very special treatment here because of coarse & irrregular AIS positioning \n",
    "    for idt, t in enumerate(t_unmod):\n",
    "       \n",
    "        # Y modification case ID 5:\n",
    "        # comes from bigger y values going to smaller y values in negative direction!\n",
    "        \n",
    "        \n",
    "        dy_distance_to_sim_base =  y_pos_unmod[idt] - ypt0_sim_base\n",
    "        \n",
    "    \n",
    "        if dy_distance_to_sim_base < sim_domain_half_width: \n",
    "            \n",
    "            if idt < idt_x:\n",
    "                print('ID 5 Y scan: found ini pos at {} earlier than X: {} '.format(idt,idt_x))\n",
    "                print('-> Hence X index determines new intial position!! Interupt Y scan.')\n",
    "                print('y_mod_y_scan[idt]   : ',y_mod_y_scan[idt])\n",
    "                print('y_mod_y_scan[idt_x] : ',y_mod_y_scan[idt_x])\n",
    "                            \n",
    "                ypt0_vessel_reset = y_mod_y_scan[idt_x]       \n",
    "            \n",
    "                y_mod_y_scan[:idt_x]  = y_mod_y_scan[idt_x]   \n",
    "                dy_mod_y_scan[:idt_x] = dy_mod_y_scan[idt_x] \n",
    "                \n",
    "            else:\n",
    "                raise ValueError('ID 5: Error resetting start pos in Y direction.')\n",
    "\n",
    "            \n",
    "            print('Carefully initialise plume/engine as f(t) here (in simulation)!')\n",
    "            \n",
    "            break\n",
    "    \n",
    "    \n",
    "    print('Resetted initial X position {} at time step {}'.format(xpt0_vessel_reset,idt_x))\n",
    "    print('Resetted initial Y position {} at time step {}'.format(ypt0_vessel_reset,idt_x))\n",
    "    \n",
    "    #- Recalculation of all displacements relative to reset initial position:\n",
    "    \n",
    "    dx_reset = [0.]  \n",
    "    dy_reset = [0.]\n",
    "\n",
    "    for idx, obj in enumerate(t_unmod):\n",
    "        if idx >0:\n",
    "        \n",
    "            dxi = x_mod_x_scan[idx] - xpt0_vessel_reset\n",
    "            dyi = y_mod_y_scan[idx] - ypt0_vessel_reset\n",
    "                \n",
    "            dx_reset.append(dxi)\n",
    "            dy_reset.append(dyi)\n",
    "            \n",
    "            \n",
    "     #- round values and format to ensure 4 decimals are returned\n",
    "    \n",
    "    dx_reset_rnd = ['%.4f' % np.around(val,4) for val in dx_reset]\n",
    "    dy_reset_rnd = ['%.4f' % np.around(val,4) for val in dy_reset] \n",
    "   \n",
    "    \n",
    "   \n",
    "    return x_mod_x_scan, y_mod_y_scan, dx_reset_rnd, dy_reset_rnd\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# ID 3: moves out eastward (\"right side\"), hence in positve x direction\n",
    "def leave_on_east_side_southbound(t_unmod, x_pos_unmod, y_pos_unmod,\\\n",
    "                                  dx_unmod, dy_unmod, x_ini, y_ini):  # ID 3 (purple)\n",
    "    \n",
    "    print('    Track filtering leave_on_east_side_southbound')\n",
    "        \n",
    "    # loop in positive X direction \n",
    "    \n",
    "    x_mod_x_scan = x_pos_unmod.copy()   \n",
    "    y_mod_y_scan = y_pos_unmod.copy()  \n",
    "    \n",
    "    dx_mod_x_scan = dx_unmod.copy()    \n",
    "    dy_mod_y_scan = dy_unmod.copy()     \n",
    "    \n",
    "    \n",
    "    # calc difference vessel start pos to overall base-pos. (=centre) of simulation domain:\n",
    "   \n",
    "    dx_vessel_start_to_sim_centre = x_ini - xpt0_sim_base   \n",
    "    \n",
    "    dy_vessel_start_to_sim_centre = y_ini - ypt0_sim_base\n",
    "    \n",
    "     \n",
    "    print('xpt0_sim_base: ',xpt0_sim_base,'  x_ini: ',x_ini)\n",
    "    print('ypt0_sim_base: ',ypt0_sim_base,'  y_ini: ',y_ini)\n",
    "\n",
    "    print('dx_vessel_start_to_sim_centre: ',dx_vessel_start_to_sim_centre,\\\n",
    "          '  T/F: ',abs(dx_vessel_start_to_sim_centre) > sim_domain_half_width)\n",
    "    print('dy_vessel_start_to_sim_centre: ',dy_vessel_start_to_sim_centre,\\\n",
    "          '  T/F: ',abs(dy_vessel_start_to_sim_centre) > sim_domain_half_width)\n",
    "    \n",
    "    \n",
    "    if abs(dx_vessel_start_to_sim_centre) > sim_domain_half_width:\n",
    "        raise ValueError('Vessel initial X position outside simulation domain: {}'.format(dx_vessel_start_to_sim_centre))\n",
    "    else:\n",
    "        print('X direction distance of vessel initial position:',dx_vessel_start_to_sim_centre)\n",
    "        \n",
    "    if abs(dy_vessel_start_to_sim_centre) > sim_domain_half_width:\n",
    "        print('Vessel initial Y position outside simulation domain: ',\\\n",
    "              dy_vessel_start_to_sim_centre)\n",
    "    else:\n",
    "        print('Vessel initial Y position inside simulation domain:',\\\n",
    "              dy_vessel_start_to_sim_centre)\n",
    "    \n",
    "    \n",
    "    \n",
    "    # scan for positions ahead:\n",
    "    \n",
    "    t_offset = 4 # 2min are likely 4 timesteps a 30sec\n",
    "    \n",
    "    \n",
    "    # loop until 4th last time step\n",
    "    for idt, t in enumerate(t_unmod[:-t_offset]):                      \n",
    "        \n",
    "        dt_4_pos_ahead = t_unmod[idt + t_offset]\n",
    "        \n",
    "        \n",
    "        # in this case ID 3:\n",
    "        # if vessel moves out of domain it comes from smaller x values going to bigger x \n",
    "        # values (in positive direction!)\n",
    "        \n",
    "        # define vessel stopping criterium: if vessel reaches sim-domain edge, all values up \n",
    "        # to end of overall time frame are replicated from 2 min away from the edge position \n",
    "        # (and loop is interrupted):\n",
    "        \n",
    "        \n",
    "        dx_distance_to_sim_base =  x_pos_unmod[idt] - xpt0_sim_base \n",
    "        \n",
    "      \n",
    "        if dx_distance_to_sim_base > sim_domain_half_width:\n",
    "        # should interrupt at idx 143:\n",
    "            \n",
    "            # modify range from 2min away of this position until end of list:\n",
    "            \n",
    "            x_mod_x_scan[idt + t_offset:]  = x_mod_x_scan[idt + t_offset]\n",
    "            dx_mod_x_scan[idt + t_offset:] = dx_mod_x_scan[idt + t_offset]\n",
    "            \n",
    "            # also set Y displacement (otherwise vessel would \"jump\" back in Y during simulation):\n",
    "            \n",
    "            y_mod_y_scan[idt + t_offset:]  = y_mod_y_scan[idt + t_offset]\n",
    "            dy_mod_y_scan[idt + t_offset:] = dy_mod_y_scan[idt + t_offset]\n",
    "             \n",
    "            # to ensure that the position 4 time steps away is not too far away (sometimes \n",
    "            # spacing is not 30sec but up to 120sec or even more), calculate the time \n",
    "            # difference between the two points (in seconds):\n",
    "            \n",
    "            t_difference_this_scan = (t_unmod[idt + t_offset] - t ) /1000.\n",
    "            \n",
    "            print('   delta time is : ',t_difference_this_scan,' sec')\n",
    "            \n",
    "            if t_difference_this_scan > 120.:\n",
    "                raise ValueError('Vessel stop time diff > 120sec: {}'.format(t_difference_this_scan))\n",
    "            \n",
    "            break\n",
    "    \n",
    "    return x_mod_x_scan,y_mod_y_scan, dx_mod_x_scan, dy_mod_y_scan\n",
    "\n",
    "\n",
    "\n",
    "# ID 1 & 4:  moves out westward (\"left side\"), in negative x direction \n",
    "# (and positive in Y)\n",
    "def leave_on_west_side_northbound(t_unmod, x_pos_unmod, y_pos_unmod,\\\n",
    "                                  dx_unmod, dy_unmod, x_ini, y_ini):  # ID 1 & 4 (red)\n",
    "    \n",
    "    print('    Track filtering leave_on_west_side_northbound')\n",
    "    \n",
    "    # loop through reversed X direction \n",
    "    \n",
    "    x_mod_x_scan = x_pos_unmod.copy()   \n",
    "    y_mod_y_scan = y_pos_unmod.copy()    \n",
    "    \n",
    "    dx_mod_x_scan = dx_unmod.copy()     \n",
    "    dy_mod_y_scan = dy_unmod.copy()      \n",
    "    \n",
    "    \n",
    "    # calc difference vessel start pos to overall base-pos. (=centre) of simulation domain:\n",
    "    \n",
    "    dx_vessel_start_to_sim_centre = xpt0_sim_base - x_ini\n",
    "    \n",
    "    dy_vessel_start_to_sim_centre = ypt0_sim_base - y_ini\n",
    "    \n",
    "    \n",
    "    # is this position outside simulation domain?\n",
    "    # if so, raise error because this shouldn't happen here (we scan for leaving domain only)\n",
    "    \n",
    "    if abs(dx_vessel_start_to_sim_centre) > sim_domain_half_width:\n",
    "        raise ValueError('Vessel initial X position outside simulation domain: {}'.format(dx_vessel_start_to_sim_centre))\n",
    "    else:\n",
    "        print('X direction distance of vessel initial position:',\\\n",
    "              dx_vessel_start_to_sim_centre)\n",
    "        \n",
    "    if abs(dy_vessel_start_to_sim_centre) > sim_domain_half_width:\n",
    "        print('Vessel initial Y position outside simulation domain: ',\\\n",
    "              dy_vessel_start_to_sim_centre)\n",
    "    else:\n",
    "        print('Vessel initial Y position inside simulation domain:',\\\n",
    "              dy_vessel_start_to_sim_centre)\n",
    "\n",
    "    \n",
    "    \n",
    "    # scan for positions ahead:\n",
    "    \n",
    "    t_offset = 4 # 2min are likely 4 timesteps a 30sec\n",
    "    \n",
    "    \n",
    "    \n",
    "    for idt, t in enumerate(t_unmod[:-t_offset]):                   \n",
    "        \n",
    "        dt_4_pos_ahead = t_unmod[idt + t_offset]\n",
    "        \n",
    "              \n",
    "        \n",
    "        dx_distance_to_sim_base =  xpt0_sim_base - x_pos_unmod[idt]\n",
    "        \n",
    "        \n",
    "        if dx_distance_to_sim_base > sim_domain_half_width: \n",
    "            \n",
    "            # modify range from 2min away of this position until end of list:\n",
    "            \n",
    "            x_mod_x_scan[idt + t_offset:]  = x_mod_x_scan[idt + t_offset]\n",
    "            dx_mod_x_scan[idt + t_offset:] = dx_mod_x_scan[idt + t_offset]\n",
    "            \n",
    "            # also set Y displacement (otherwise vessel would \"jump\" back in Y during \n",
    "            # simulation):\n",
    "            \n",
    "            y_mod_y_scan[idt + t_offset:]  = y_mod_y_scan[idt + t_offset]\n",
    "            dy_mod_y_scan[idt + t_offset:] = dy_mod_y_scan[idt + t_offset]\n",
    "          \n",
    "            t_difference_this_scan = (t_unmod[idt + t_offset] - t ) /1000.\n",
    "                        \n",
    "            if t_difference_this_scan > 120.:\n",
    "                raise ValueError('Vessel stop time diff > 120sec: {}'.format(t_difference_this_scan))\n",
    "            \n",
    "            break\n",
    "    \n",
    "    return x_mod_x_scan,y_mod_y_scan, dx_mod_x_scan, dy_mod_y_scan\n",
    "\n",
    "\n",
    "    \n",
    "# ID 2:  not entering or leaving sim-domain - but stop positions need to be corrected:\n",
    "def set_last_pos_dxdy(t_unmod, x_pos_unmod, y_pos_unmod,\\\n",
    "                      dx_unmod, dy_unmod, x_ini, y_ini):        # ID 2\n",
    "    print('    Setting stop positions.')\n",
    "    \n",
    "    \n",
    "    dx_mod_x_scan = dx_unmod.copy()     \n",
    "    dy_mod_y_scan = dy_unmod.copy()      \n",
    "    \n",
    "   \n",
    "    \n",
    "    # scan for detecting stop position:\n",
    "    for idt, t in enumerate(t_unmod):  #[:-1]\n",
    "        \n",
    "        # stop pos is characterised as no deviation from mean from to stop to \n",
    "        # end of time series!\n",
    "        \n",
    "        avg_x_eots = np.mean(x_pos_unmod[idt:])\n",
    "        avg_y_eots = np.mean(y_pos_unmod[idt:])\n",
    "                \n",
    "        \n",
    "        # compare truncated values to avoid numerical precision problems:\n",
    "\n",
    "        if '%.6f' % x_pos_unmod[idt] == '%.6f' % avg_x_eots and \\\n",
    "           '%.6f' % y_pos_unmod[idt] == '%.6f' % avg_y_eots:  \n",
    "       \n",
    "            print('ID 2: stop position detected at :',\\\n",
    "                  idt,t_unmod[idt],x_pos_unmod[idt],y_pos_unmod[idt])\n",
    "            \n",
    "            dx_mod_x_scan[idt:] = dx_mod_x_scan[idt] \n",
    "            dy_mod_y_scan[idt:] = dy_mod_y_scan[idt]\n",
    "            \n",
    "            break\n",
    "        else:   \n",
    "            if idt == len(t_unmod)-1:\n",
    "                raise ValueError('ID 2: No vessel final stop position detected')\n",
    "\n",
    "    \n",
    "    return x_pos_unmod, y_pos_unmod, dx_mod_x_scan, dy_mod_y_scan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def write_files(filelist, varlist):\n",
    "\n",
    "    for item,obj in enumerate(filelist):\n",
    "        with open(outpath + obj,'w') as f:\n",
    "            writer = csv.writer(f)\n",
    "            for val in varlist[item]:\n",
    "                writer.writerow(val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General definitions wrt time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial time : 2018-06-28 09:55:00    unix time stamp: 1530172500000\n",
      "Overall stop time:  2018-06-28 14:00:00 <class 'datetime.datetime'> 2018-06-28 14:00:00\n",
      "                                      unix time stamp: 1530187200000.0\n"
     ]
    }
   ],
   "source": [
    "# start 5 min to 10:00 in order to establish 1st vessels exhaust plume\n",
    "\n",
    "t_initial = datetime.datetime(2018, 6, 28, 9, 55)\n",
    "print ('Initial time :',t_initial.isoformat(' '),'   unix time stamp:',\\\n",
    "       time_datetime_to_unixtime(t_initial))\n",
    "\n",
    "\n",
    "# define fixed stop time for vessel motion:\n",
    "\n",
    "t_stop_overall = datetime.datetime(2018, 6, 28, 14, 00)\n",
    "print ('Overall stop time: ',t_stop_overall,type(t_stop_overall),|\n",
    "       t_stop_overall.strftime('%Y-%m-%d %H:%M:%S'))\n",
    "print ('                                      unix time stamp:',|\n",
    "       t_stop_overall.timestamp() * 1000)\n",
    "\n",
    "\n",
    "# define time window for searching vessel motion adjustment in order \n",
    "# not allowing a vessel to leave the computational domain no more than \n",
    "# two minutes away (maybe 2min is too much)\n",
    "track_correction_time_window = 2.*60.*1000.  # milliseconds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Process all files and tracks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read file : /Users/rene/Documents/Environmod_iMac/Projekte2/Bullemer_Harbour_Project/Data_Akio/Set2/Python/csv/vesselmove_id211443180_dt20180629_v0.csv\n",
      "id211443180 :    1530172800000 9.955818 53.490818\n",
      "Vessel ID:  id211443180 :   time to add prior to track:  300000 ms ( 5.0  min)\n",
      "Vessel ID:  id211443180 :   time to add subsequently to track:  10290000 ms ( 171.5  min)\n",
      "Initial time                 (t_initial)  : 1530172500000\n",
      "Vessels start time           (base_time)  : 1530172800000\n",
      "...\n",
      "Vessels last time steps time (t_unix last): 1530176910000\n",
      "Overall stop time            (t_initial)  : 1530187200000.0\n",
      "* * * Cartesian mapping of all track positions:\n",
      "Global cartesian coords:\n",
      "  Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y: 98133.64755970874 58975.38045677248\n",
      "  1st pos of vessel X:  97770.16883059918   last pos X:  82806.93333683792\n",
      "  1st pos of vessel Y:  53571.616008122386   last pos Y:  61824.08800349565\n",
      "  min/max deviations:  in X: -14963.235493761269 -40.973636532755336   in Y: 71.81118883905583 8252.471995373264\n",
      "\n",
      " * * * Running Order:  1\n",
      "    Track filtering leave_on_west_side_northbound\n",
      "X direction distance of vessel initial position: 363.4787291095563\n",
      "Vessel initial Y position outside simulation domain:  5403.764448650094\n",
      " \n",
      "read file : /Users/rene/Documents/Environmod_iMac/Projekte2/Bullemer_Harbour_Project/Data_Akio/Set2/Python/csv/vesselmove_id211590990_dt20180629_v0.csv\n",
      "id211590990 :    1530172800000 10.003243 53.537183\n",
      "Vessel ID:  id211590990 :   time to add prior to track:  300000 ms ( 5.0  min)\n",
      "Vessel ID:  id211590990 :   time to add subsequently to track:  2550000 ms ( 42.5  min)\n",
      "Initial time                 (t_initial)  : 1530172500000\n",
      "Vessels start time           (base_time)  : 1530172800000\n",
      "...\n",
      "Vessels last time steps time (t_unix last): 1530184650000\n",
      "Overall stop time            (t_initial)  : 1530187200000.0\n",
      "* * * Cartesian mapping of all track positions:\n",
      "Global cartesian coords:\n",
      "  Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y: 98133.64755970874 58975.38045677248\n",
      "  1st pos of vessel X:  100917.6919531231   last pos X:  98794.40446299927\n",
      "  1st pos of vessel Y:  58730.9475883466   last pos Y:  57257.09807905132\n",
      "  min/max deviations:  in X: -4614.351483123464 -313.93416075318237   in Y: -3520.940116241858 817.9177589868705\n",
      "\n",
      " * * * Running Order:  2\n",
      "    Setting stop positions.\n",
      "ID 2: stop position detected at : 132 1530184650000 98794.404463 57257.0980791\n",
      " \n",
      "read file : /Users/rene/Documents/Environmod_iMac/Projekte2/Bullemer_Harbour_Project/Data_Akio/Set2/Python/csv/vesselmove_id211177630_dt20180629_v0.csv\n",
      "id211177630 :    1530174480000 9.943532 53.515927\n",
      "Vessel ID:  id211177630 :   time to add prior to track:  1980000 ms ( 33.0  min)\n",
      "Vessel ID:  id211177630 :   time to add subsequently to track:  9990000 ms ( 166.5  min)\n",
      "Initial time                 (t_initial)  : 1530172500000\n",
      "Vessels start time           (base_time)  : 1530174480000\n",
      "...\n",
      "Vessels last time steps time (t_unix last): 1530177210000\n",
      "Overall stop time            (t_initial)  : 1530187200000.0\n",
      "* * * Cartesian mapping of all track positions:\n",
      "Global cartesian coords:\n",
      "  Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y: 98133.64755970874 58975.38045677248\n",
      "  1st pos of vessel X:  96956.91948880129   last pos X:  104405.2115375147\n",
      "  1st pos of vessel Y:  56366.71563871107   last pos Y:  56673.302168695765\n",
      "  min/max deviations:  in X: -659.7027396964259 7448.292048713411   in Y: 154.79126906157035 3111.094198034698\n",
      "\n",
      " * * * Running Order:  3\n",
      "    Track filtering leave_on_east_side_southbound\n",
      "xpt0_sim_base:  98133.64755970874   x_ini:  96956.91948880129\n",
      "ypt0_sim_base:  58975.38045677248   y_ini:  56366.71563871107\n",
      "dx_vessel_start_to_sim_centre:  -1176.7280709074548   T/F:  False\n",
      "dy_vessel_start_to_sim_centre:  -2608.6648180614065   T/F:  False\n",
      "X direction distance of vessel initial position: -1176.7280709074548\n",
      "Vessel initial Y position inside simulation domain: -2608.6648180614065\n",
      "   delta time is :  120.0  sec\n",
      " \n",
      "read file : /Users/rene/Documents/Environmod_iMac/Projekte2/Bullemer_Harbour_Project/Data_Akio/Set2/Python/csv/vesselmove_id211202770_dt20180629_v0.csv\n",
      "id211202770 :    1530175890000 10.031312 53.520415\n",
      "Vessel ID:  id211202770 :   time to add prior to track:  3390000 ms ( 56.5  min)\n",
      "Vessel ID:  id211202770 :   time to add subsequently to track:  2070000 ms ( 34.5  min)\n",
      "Initial time                 (t_initial)  : 1530172500000\n",
      "Vessels start time           (base_time)  : 1530175890000\n",
      "...\n",
      "Vessels last time steps time (t_unix last): 1530185130000\n",
      "Overall stop time            (t_initial)  : 1530187200000.0\n",
      "* * * Cartesian mapping of all track positions:\n",
      "Global cartesian coords:\n",
      "  Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y: 98133.64755970874 58975.38045677248\n",
      "  1st pos of vessel X:  102779.51425716786   last pos X:  69499.74555774448\n",
      "  1st pos of vessel Y:  56865.18431387271   last pos Y:  70199.96262355089\n",
      "  min/max deviations:  in X: -33279.76869942338 -559.5205833170185   in Y: 1397.444551739245 13334.778309678179\n",
      "\n",
      " * * * Running Order:  4\n",
      "    Track filtering leave_on_west_side_northbound\n",
      "X direction distance of vessel initial position: -4645.866697459118\n",
      "Vessel initial Y position inside simulation domain: 2110.1961428997674\n",
      " \n",
      "read file : /Users/rene/Documents/Environmod_iMac/Projekte2/Bullemer_Harbour_Project/Data_Akio/Set2/Python/csv/vesselmove_id244650903_dt20180629_v0.csv\n",
      "id244650903 :    1530181590000 9.59428 53.596177\n",
      "Vessel ID:  id244650903 :   time to add prior to track:  9090000 ms ( 151.5  min)\n",
      "Vessel ID:  id244650903 :   time to add subsequently to track:  30000 ms ( 0.5  min)\n",
      "Initial time                 (t_initial)  : 1530172500000\n",
      "Vessels start time           (base_time)  : 1530181590000\n",
      "...\n",
      "Vessels last time steps time (t_unix last): 1530187170000\n",
      "Overall stop time            (t_initial)  : 1530187200000.0\n",
      "* * * Cartesian mapping of all track positions:\n",
      "Global cartesian coords:\n",
      "  Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y: 98133.64755970874 58975.38045677248\n",
      "  1st pos of vessel X:  73840.51717902496   last pos X:  96758.29040971972\n",
      "  1st pos of vessel Y:  65373.340781673316   last pos Y:  54985.472461581674\n",
      "  min/max deviations:  in X: 1591.3950064704113 22917.77323069476   in Y: -10387.868320091642 -1637.339909970724\n",
      "\n",
      " * * * Running Order:  5\n",
      "    Track filtering enter_from_west_southbound\n",
      "xpt0_sim_base:  98133.64755970874   x_ini:  73840.51717902496\n",
      "ypt0_sim_base:  58975.38045677248   y_ini:  65373.340781673316\n",
      "dx_vessel_start_to_sim_centre:  24293.130380683782   T/F:  True\n",
      "dy_vessel_start_to_sim_centre:  -6397.960324900836   T/F:  True\n",
      "Vessel initial X position outside simulation domain:  24293.130380683782\n",
      "Vessel initial Y position outside simulation domain:  -6397.960324900836 SPECIAL FILTERING NECESSARY !!!!!!!!!\n",
      "   delta time is :  -120.0  sec\n",
      ". . . . . . . . . . . Y  modification next . . . . . . . . . . . . . . . . . . . . .\n",
      "ID 5 Y scan: found intial position at 304 earlier than in X: 351 \n",
      "-> Hence X index determines new intial position!! Interupt Y scan.\n",
      "y_mod_y_scan[idt]   :  63736.0008717\n",
      "y_mod_y_scan[idt_x] :  59262.0594659\n",
      "Carefully initialise plume/engine as f(t) here (in simulation)! ! ! ! ! ! ! ! ! !\n",
      "Resetted initial X position 92558.7070472971 at time step 351\n",
      "Resetted initial Y position 59262.05946585979 at time step 351\n",
      " \n",
      "read file : /Users/rene/Documents/Environmod_iMac/Projekte2/Bullemer_Harbour_Project/Data_Akio/Set2/Python/csv/vesselmove_id211177670_dt20180629_v0.csv\n",
      "id211177670 :    1530183870000 10.03393 53.529852\n",
      "Vessel ID:  id211177670 :   time to add prior to track:  11370000 ms ( 189.5  min)\n",
      "Vessel ID:  id211177670 :   time to add subsequently to track:  30000 ms ( 0.5  min)\n",
      "Initial time                 (t_initial)  : 1530172500000\n",
      "Vessels start time           (base_time)  : 1530183870000\n",
      "...\n",
      "Vessels last time steps time (t_unix last): 1530187170000\n",
      "Overall stop time            (t_initial)  : 1530187200000.0\n",
      "* * * Cartesian mapping of all track positions:\n",
      "Global cartesian coords:\n",
      "  Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y: 98133.64755970874 58975.38045677248\n",
      "  1st pos of vessel X:  102952.65864147633   last pos X:  97108.89033536252\n",
      "  1st pos of vessel Y:  57915.56534595281   last pos Y:  54445.074858300446\n",
      "  min/max deviations:  in X: -6695.726810052132 -219.08549613402283   in Y: -3470.490487652365 1650.1042204881887\n",
      "\n",
      " * * * Running Order:  6\n",
      "    Track filtering enter_from_east_northbound\n",
      "xpt0_sim_base:  98133.64755970874   x_ini:  102952.65864147633\n",
      "ypt0_sim_base:  58975.38045677248   y_ini:  57915.56534595281\n",
      "dx_vessel_start_to_sim_centre:  -4819.011081767589   T/F:  False\n",
      "dy_vessel_start_to_sim_centre:  1059.815110819669   T/F:  False\n",
      "X direction distance of vessel initial position: -4819.011081767589\n",
      "Vessel initial Y position inside simulation domain: 1059.815110819669\n",
      " \n"
     ]
    }
   ],
   "source": [
    "running_order = 1          #... human readable\n",
    "\n",
    "for id in vessel_id:\n",
    "    \n",
    "    ## 1) read each AIS file\n",
    "    ##    Code change (compared to prepare_xy_positionings_displacements_03.ipynb): renamed \n",
    "    ##    x_geo/y_geo into lon_ais, lat_ais (in def above old names still used)\n",
    "    \n",
    "    t_unix, lon_ais, lat_ais = vessel_read(inpath+'vesselmove_'+ id +'_dt20180629_v0.csv')  \n",
    "    \n",
    "    \n",
    "    ## 2) set base coords and time:  (little code change to simplify handling)\n",
    "\n",
    "    lon_ais_0 = lon_ais[0]\n",
    "    lat_ais_0 = lat_ais[0]\n",
    "    base_time = t_unix[0]\n",
    "\n",
    "    print(id,':   ',t_unix[0],lon_ais[0],lat_ais[0])\n",
    "    \n",
    "    \n",
    "    \n",
    "    ## 3) Convert time to store in seconds since begin (list comprehension):\n",
    "    ##    minor code change: datetime.datetime now\n",
    "    \n",
    "    t_hr=[datetime.datetime.fromtimestamp(int(t)/1000).strftime('%Y-%m-%d %H:%M:%S') for t in t_unix]\n",
    "        \n",
    "\n",
    "    ## 4) check if there are time intervalls != 30sec (important for temporal 'mapping'):\n",
    "    #inspect_time_intervals(t_unix)\n",
    "    \n",
    "    \n",
    "    \n",
    "    ## 5) temporal mapping of track's time frame into the overall time frame of simulation:\n",
    "    ##    this is not a true mapping - here just the period/timeseries up to t_initial is added \n",
    "    ##    at the beginning, prior to starting time of vessel and, respectively, and at the end \n",
    "    ##    we extend the time to t_stop_overall\n",
    "    \n",
    "    #   5.1) calc difference between start time of vessel and initial time (in ms):\n",
    "\n",
    "    t_delta_begin_ms = int(base_time) - time_datetime_to_unixtime(t_initial)\n",
    "        \n",
    "    print('Vessel ID: ',id,':   time to add prior to track: ',t_delta_begin_ms,\\\n",
    "          'ms (',t_delta_begin_ms/1000/60,' min)')\n",
    "\n",
    "    \n",
    "    #   5.2) calc difference between last time stamp in data set and overall stop time (in ms):\n",
    "    \n",
    "\n",
    "    ## NOT USED LATER (direct defintiuon of time array below)\n",
    "    t_delta_subsequently_ms = time_datetime_to_unixtime(t_stop_overall) - int(t_unix[len(t_unix)-1])\n",
    "    \n",
    "    \n",
    "    print('Vessel ID: ',id,':   time to add subsequently to track: ',t_delta_subsequently_ms,\\\n",
    "          'ms (',t_delta_subsequently_ms/1000/60,' min)')\n",
    "\n",
    "    \n",
    "    #   5.3) set time frames to be added in 30sec intervalls:\n",
    "        \n",
    "    time_frame_prior_unixtime = np.arange(\n",
    "                                        time_datetime_to_unixtime(t_initial), \n",
    "                                        time_datetime_to_unixtime(t_initial) + t_delta_begin_ms, \n",
    "                                        30000)\n",
    "    \n",
    "    print('Initial time                 (t_initial)  :',time_datetime_to_unixtime(t_initial))\n",
    "    print('Vessels start time           (base_time)  :',base_time)\n",
    "\n",
    "       \n",
    "    # here one need to add 30 sec offset to last time step in data file AND 30sec to end of array \n",
    "    time_frame_subsequent_unixtime = np.arange(\n",
    "                                        int(t_unix[len(t_unix)-1]) + 30000,        # ms !!!\n",
    "                                        int(t_stop_overall.timestamp() * 1000 + 30000), \n",
    "                                        30000)\n",
    "    \n",
    "    print('Vessels last time steps time (t_unix last):',int(t_unix[len(t_unix)-1]))\n",
    "    print('Overall stop time            (t_initial)  :',t_stop_overall.timestamp() * 1000)\n",
    "\n",
    "   \n",
    "    \n",
    "    ## 6) mapping of each vessel track's initial position: geo-coords --> global cartesian \n",
    "    #    (better: projection region cartesian)\n",
    "    ##   Note that lon,lat can be scalars, lists or numpy arrays\n",
    "    \n",
    "    xpt0_vessel,ypt0_vessel = m(lon_ais_0,lat_ais_0)\n",
    "    \n",
    "       \n",
    "    \n",
    "    ## 7) convert all positions to global cartesian coords (cartesian mapping of all \n",
    "    #   track positions):\n",
    "    ##    for better debuggin no list comprehension here\n",
    "    ##    (code changed using lon_ais, lat_ais)\n",
    "    \n",
    "    x_glob_cartesian = []\n",
    "    y_glob_cartesian = []\n",
    "\n",
    "    \n",
    "    print('* * * Cartesian mapping of all track positions:')\n",
    "    \n",
    "    for idx, obj in enumerate(lon_ais):\n",
    "\n",
    "        xpt,ypt = m(lon_ais[idx],lat_ais[idx])\n",
    "    \n",
    "        x_glob_cartesian.append(xpt)\n",
    "        y_glob_cartesian.append(ypt)\n",
    "    \n",
    "   \n",
    "    # store last position to fill up overall time frame should the vesser stop before \n",
    "    # (which is the case for all vessels)\n",
    "    # Note: this is a value corresponding to last entry of vessel track file!\n",
    "    \n",
    "    xpt_last_vessel = x_glob_cartesian[len(x_glob_cartesian)-1]\n",
    "    ypt_last_vessel = y_glob_cartesian[len(y_glob_cartesian)-1]\n",
    "    \n",
    "    print('Global cartesian coords:')\n",
    "    print('  Cartesian base coordinate of simulation domain (9.961250/53.539373)  x/y:',\\\n",
    "          xpt0_sim_base, ypt0_sim_base)  \n",
    "    print('  1st pos of vessel X: ',xpt0_vessel,'  last pos X: ',xpt_last_vessel)\n",
    "    print('  1st pos of vessel Y: ',ypt0_vessel,'  last pos Y: ',ypt_last_vessel)\n",
    "\n",
    "    \n",
    "    \n",
    "    ## 8) calculate change from each position to REFERENCE position as delta_x, delta_y:\n",
    "    ##    note 1: reference is vessels initial position! \n",
    "    ##    note 2: resulting vectors do not have entries for start position! --> first position's \n",
    "    #             displacement is zero!!! this nereds to be taken into account when later fill\n",
    "    #             up the positions for the whole time frame!\n",
    "\n",
    "    dx = []\n",
    "    dy = []\n",
    "\n",
    "    distances = [] # from origin !!!\n",
    "\n",
    "    for idx, obj in enumerate(lon_ais):\n",
    "        if idx >0:\n",
    "        \n",
    "            dxi = x_glob_cartesian[idx] - xpt0_vessel \n",
    "            dyi = y_glob_cartesian[idx] - ypt0_vessel \n",
    "                \n",
    "            dx.append(dxi)\n",
    "            dy.append(dyi)\n",
    "        \n",
    "            # distance according to pythagoras:\n",
    "            distances.append(math.sqrt(dxi**2 + dyi**2))\n",
    "        \n",
    "            #print(idx,dxi,'  ',dyi,'  ',math.sqrt(dxi**2 + dyi**2))   \n",
    "    \n",
    "    \n",
    "    ## 9) calc. max deviations from origin to later diefine size of simulation domain:\n",
    "    ##    Note: when a ship is moving too far, it will be stopped when it leaves a core \n",
    "    ##          simulation domain - that is the same domain that is later used in \n",
    "    #           postprocessing\n",
    "    ##\n",
    "    ##    These deviations are stored in separete file\n",
    "    ##\n",
    "    \n",
    "    min_x_displacement_total = min(dx)\n",
    "    max_x_displacement_total = max(dx)\n",
    "    min_y_displacement_total = min(dy)\n",
    "    max_y_displacement_total = max(dy)\n",
    "    \n",
    "    print('  min/max deviations:  in X:',min(dx),max(dx),'  in Y:',min(dy),max(dy))\n",
    "    \n",
    "    \n",
    "    \n",
    "    ## 10) round numbers\n",
    "    \n",
    "    #dx_rnd = [round(val,4) for val in dx]\n",
    "    #dy_rnd = [round(val,4) for val in dy]\n",
    "     \n",
    "    #- round values and format to ensure 4 decimals are returned\n",
    "    \n",
    "    dx_rnd = ['%.4f' % np.around(val,4) for val in dx]\n",
    "    dy_rnd = ['%.4f' % np.around(val,4) for val in dy] \n",
    "    \n",
    "    \n",
    "    \n",
    "    ## 11) concatenate prior, real, subsequent time frames:\n",
    "    \n",
    "    debug_this = False\n",
    "    \n",
    "    if debug_this:\n",
    "        print('Times to be merged: \\n',\n",
    "            time_frame_prior_unixtime[len(time_frame_prior_unixtime)-1],'\\n',\n",
    "            t_unix[0],'\\n',\n",
    "            int(t_unix[len(t_unix)-1]),'\\n',\n",
    "            time_frame_subsequent_unixtime[0])\n",
    "    \n",
    "        print('Times to be merged: \\n',\n",
    "            type(time_frame_prior_unixtime[len(time_frame_prior_unixtime)-1]),'\\n',\n",
    "            type(int(t_unix[0])),'\\n',\n",
    "            type(  np.int64(int(t_unix[len(t_unix)-1]))  ),'\\n',\n",
    "            type(time_frame_subsequent_unixtime[0]))\n",
    "   \n",
    "        print('Times to be merged: \\n',\n",
    "            type(time_frame_prior_unixtime),'\\n',\n",
    "            type(t_unix),'\\n',\n",
    "            type(np.int64(t_unix)),'\\n',       # should be converted to numpy int array\n",
    "            type(time_frame_subsequent_unixtime))\n",
    "    \n",
    "    \n",
    "    \n",
    "    #   11.1) convert the native python list into ndarray containing np.int64\n",
    "    #        (not converting vice versa because this ways is more convenient)\n",
    "    \n",
    "    tmp = [ np.int64(item)  for item in t_unix]\n",
    "    \n",
    "    time_frame_vesselmotion_unixtime = np.array(tmp)\n",
    "    \n",
    "    \n",
    "    #   11.2) concatenate:\n",
    "    \n",
    "    t_track_merged = np.concatenate([time_frame_prior_unixtime, \\\n",
    "                                     time_frame_vesselmotion_unixtime, \\\n",
    "                                     time_frame_subsequent_unixtime])\n",
    "  \n",
    "    \n",
    "    \n",
    "    \n",
    "        \n",
    "    ## 12) position the vessel in frames prior and subsequently to given track:\n",
    "    \n",
    "    #  12.1) positions in cartesian coords - needed to later correct for simulation  \n",
    "    #        if vessel goes to far or enters the simualtion domain: \n",
    "    \n",
    "    #  12.1.1) set starting point-filled arrays in length of the two artifical \n",
    "    #          time frames, apply to both directions    \n",
    "    \n",
    "    x_pos_prior      = np.empty(len(time_frame_prior_unixtime))\n",
    "    x_pos_prior.fill(xpt0_vessel)\n",
    "\n",
    "    y_pos_prior      = np.empty(len(time_frame_prior_unixtime))\n",
    "    y_pos_prior.fill(ypt0_vessel)\n",
    "    \n",
    "    x_pos_subsequent = np.empty(len(time_frame_subsequent_unixtime))\n",
    "    x_pos_subsequent.fill(xpt_last_vessel)\n",
    "    \n",
    "    y_pos_subsequent = np.empty(len(time_frame_subsequent_unixtime))\n",
    "    y_pos_subsequent.fill(ypt_last_vessel)\n",
    "    \n",
    "\n",
    "   \n",
    "\n",
    "    #  12.1.2) concatenate\n",
    "    \n",
    "    x_pos_merged = np.concatenate([ x_pos_prior, x_glob_cartesian , x_pos_subsequent ])\n",
    "    y_pos_merged = np.concatenate([ y_pos_prior, y_glob_cartesian , y_pos_subsequent ])   \n",
    "\n",
    "    \n",
    "    \n",
    "    #  12.2) displacements dx & dy: \n",
    "    \n",
    "    #  12.2.1) set 0-filled arrays in length of the two artifical time frames \n",
    "    #         (0== no displacement), apply to both directions\n",
    "    \n",
    "    xy_pos_prior      = np.zeros(len(time_frame_prior_unixtime))\n",
    "    xy_pos_subsequent = np.zeros(len(time_frame_subsequent_unixtime))\n",
    "    \n",
    "   \n",
    "    #  12.2.2) conversion list to ndarray\n",
    "    tmpx = [ '%.4f' % np.around(np.float64(item),4)  for item in dx_rnd]\n",
    "    tmpy = [ '%.4f' % np.around(np.float64(item),4)  for item in dy_rnd]\n",
    "   \n",
    "    dx_rnd_np = np.array(tmpx)\n",
    "    dy_rnd_np = np.array(tmpy)\n",
    "    \n",
    "    \n",
    "    # debugging:\n",
    "    \n",
    "    debug_this = False\n",
    "    \n",
    "    if debug_this:\n",
    "        print('global cartesian X:')\n",
    "        for idx, obj in enumerate(x_glob_cartesian[:-1]):\n",
    "            if idx > 0:\n",
    "                print('idx/x_glob_cartesian/dx : ',idx,obj,dx_rnd_np[idx])\n",
    "        print('- - - - - -')\n",
    "\n",
    "    \n",
    "    \n",
    "    #  12.2.3) concatenate\n",
    "    #          note, the 0 displacement for vessels starting position (note that they are calculated \n",
    "    #          only for positions > intial!)\n",
    "    #                                           \n",
    "    dx_merged = np.concatenate([ xy_pos_prior, np.zeros(1), dx_rnd_np , xy_pos_subsequent ])\n",
    "    dy_merged = np.concatenate([ xy_pos_prior, np.zeros(1), dy_rnd_np , xy_pos_subsequent ])\n",
    "    \n",
    "\n",
    "    \n",
    "    ## 13) check if corrections on the vessel tracks are necessary \n",
    "    ##     e.g. stop the vessel when moving to far away or shift forward the position when entering the domain\n",
    "\n",
    "    \n",
    "    # track filter assignments (ID 2 not filtered): \n",
    "    filter_assignment = { 1 : leave_on_west_side_northbound,\n",
    "                          4 : leave_on_west_side_northbound,\n",
    "                          2 : set_last_pos_dxdy,\n",
    "                          3 : leave_on_east_side_southbound,\n",
    "                          5 : enter_from_west_southbound,\n",
    "                          6 : enter_from_east_northbound,\n",
    "                        }\n",
    "    \n",
    "    print('\\n * * * Running Order: ',running_order)\n",
    "    \n",
    "    # pythonic case switch:\n",
    "    x_pos_filtered, y_pos_filtered, dx_filtered , dy_filtered = filter_assignment[running_order](t_track_merged,\\\n",
    "                                                                   x_pos_merged, y_pos_merged,\\\n",
    "                                                                   dx_merged,    dy_merged,\\\n",
    "                                                                   xpt0_vessel,  ypt0_vessel)\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    ## 14) zip relative time in sec and each modified list\n",
    "    ##     this is done for both types of output: \n",
    "    ##     - the unmodified tracks & times\n",
    "    ##     - and the simulation domain adapted tracks & times\n",
    "    \n",
    "    \n",
    "    #  14.1) tracks:\n",
    "    \n",
    "    # unmodified:  x_pos_merged   & y_pos_merged\n",
    "    # modified:    x_pos_filtered & y_pos_filtered\n",
    "    \n",
    "    t_x_pos_unmod , t_y_pos_unmod = [] , []\n",
    "    t_x_pos_modif , t_y_pos_modif = [] , []\n",
    "\n",
    "   \n",
    "    # unmodified:\n",
    "    for val in zip( t_track_merged ,x_pos_merged ):    # X\n",
    "       #print(val)\n",
    "        t_x_pos_unmod.append(val)\n",
    "    \n",
    "    for val in zip( t_track_merged ,y_pos_merged):     # Y\n",
    "        t_y_pos_unmod.append(val)\n",
    "    \n",
    "    # modified:\n",
    "    for val in zip( t_track_merged ,x_pos_filtered ):  # X\n",
    "        #print(val)\n",
    "        t_x_pos_modif.append(val)\n",
    "    \n",
    "    for val in zip( t_track_merged ,y_pos_filtered):   # Y\n",
    "        t_y_pos_modif.append(val)\n",
    "    \n",
    "    \n",
    "    \n",
    "    #  14.2) displacements:\n",
    "    \n",
    "    # unmodified:  dx_merged   & dy_merged\n",
    "    # modified:    dx_filtered & dy_filtered\n",
    "    \n",
    "    t_dx_unmod , t_dy_unmod = [] , []\n",
    "    t_dx_modif , t_dy_modif = [] , []\n",
    "    \n",
    "    \n",
    "    # unmodified:\n",
    "    for val in zip( t_track_merged ,dx_merged ):    # X\n",
    "        #print(val)\n",
    "        t_dx_unmod.append(val)\n",
    "    \n",
    "    for val in zip( t_track_merged ,dy_merged):     # Y\n",
    "        t_dy_unmod.append(val)\n",
    "    \n",
    "    # modified:\n",
    "    for val in zip( t_track_merged ,dx_filtered ):  # X\n",
    "        #print(val)\n",
    "        t_dx_modif.append(val)\n",
    "    \n",
    "    for val in zip( t_track_merged ,dy_filtered):   # Y\n",
    "        t_dy_modif.append(val)\n",
    "    \n",
    "\n",
    "    \n",
    "    ## 15) write files\n",
    "    \n",
    "    fID = str(running_order)\n",
    "    \n",
    "    #  15.1) tracks:\n",
    "    \n",
    "    s1x = 't_x_pos_unmod_ID'\n",
    "    s1y = 't_y_pos_unmod_ID'\n",
    "    s2x = 't_x_pos_modif_ID'\n",
    "    s2y = 't_y_pos_modif_ID'\n",
    "   \n",
    "    filelist = [ s1x+fID+'.dat' , s1y+fID+'.dat' , s2x+fID+'.dat' , s2y+fID+'.dat' ]\n",
    "    \n",
    "    varlist  = [ t_x_pos_unmod , t_y_pos_unmod , t_x_pos_modif , t_y_pos_modif ]\n",
    "    \n",
    "    write_files(filelist, varlist)\n",
    "\n",
    "    \n",
    "    #  15.2) displacements:\n",
    "    \n",
    "    s1x = 't_dx_unmod_ID'\n",
    "    s1y = 't_dy_unmod_ID'\n",
    "    s2x = 't_dx_modif_ID'\n",
    "    s2y = 't_dy_modif_ID'\n",
    "        \n",
    "   \n",
    "    filelist = [ s1x+fID+'.dat' , s1y+fID+'.dat' , s2x+fID+'.dat' , s2y+fID+'.dat' ]\n",
    "    \n",
    "    varlist  = [ t_dx_unmod , t_dy_unmod , t_dx_modif , t_dy_modif ]\n",
    "    \n",
    "    write_files(filelist, varlist)\n",
    "\n",
    "\n",
    "    running_order +=1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}