Flood_Estimate_from_Surge_Height.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flood Estimate Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code produces flood heights for given storm parameters. The code calculates surface volume functions that can relate the volume of water to the flood height at a given division. Using Mannning's equation, it takes topographic data and storm parameter data to find the velocity and subsequently the volume of the water. This volumes are then redistributed by propagating to the nearby divisions, and the final flood heights for each division are determined.\n",
    "\n",
    "\n",
    "Inputs: surface data (topography, roughness, slope, surface volume, and surface volume grouped), flood data for cold and warm storms (peak, surge, time), and wall data for cold and warm storms.\n",
    "\n",
    "Outputs: Flood heights for cold and warm storms with and without flood wall protection in csv format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports and Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-30T16:59:46.177495Z",
     "start_time": "2020-11-30T16:59:40.172283Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import time\n",
    "import csv\n",
    "import glob                                                                                                \n",
    "import datetime\n",
    "from scipy import stats\n",
    "from scipy.stats import genpareto\n",
    "from scipy.stats import lognorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-30T16:59:46.205381Z",
     "start_time": "2020-11-30T16:59:46.190431Z"
    }
   },
   "outputs": [],
   "source": [
    "# Surge height\n",
    "Sandy_h = 2.864\n",
    "ssh = Sandy_h\n",
    "\n",
    "# Sea Level Rise Severe Cases \n",
    "# ssh = Sandy_h + 0.253 #2020s\n",
    "# ssh = Sandy_h + 0.762 #2050s\n",
    "#ssh = Sandy_h + 1.905 #2100\n",
    "\n",
    "floodheight_div = np.full(18, ssh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-30T16:59:48.008409Z",
     "start_time": "2020-11-30T16:59:46.229317Z"
    }
   },
   "outputs": [],
   "source": [
    "# essentials\n",
    "ftm = 0.3048 # feet to meter\n",
    "nt  = 10     # split number of time history of storm surge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-30T16:59:48.359749Z",
     "start_time": "2020-11-30T16:59:48.022373Z"
    }
   },
   "outputs": [],
   "source": [
    "## Functions\n",
    "\n",
    "# curve fitting\n",
    "from scipy.optimize import curve_fit\n",
    "def func_fit(surfaceV,a,b):\n",
    "    return a*np.sqrt(surfaceV)+b*surfaceV\n",
    "\n",
    "def FloodParam(ssh,nt,loc,beta1_1_c,beta1_2_c,beta2_1_c,beta2_2_c,beta1_1_w,\n",
    "               beta1_2_w,beta2_1_w,beta2_2_w,alpha_c,beta_c,alpha_w,beta_w):\n",
    "    mu_c = beta1_1_c*ssh + beta1_2_c\n",
    "    sigma_c = abs(np.sqrt(beta2_1_c*ssh + beta2_2_c))\n",
    "    mu2_c = np.log(mu_c/np.sqrt(1+sigma_c**2/mu_c**2))\n",
    "    vari = np.log(1+sigma_c**2/mu_c**2)\n",
    "    sigma2_c = np.sqrt(vari)\n",
    "    \n",
    "    mu_w = beta1_1_w*ssh + beta1_2_w\n",
    "    sigma_w = abs(np.sqrt(beta2_1_w*ssh + beta2_2_w))\n",
    "    mu2_w = np.log(mu_w/np.sqrt(1+sigma_w**2/mu_w**2))\n",
    "    vari = np.log(1+sigma_w**2/mu_w**2)\n",
    "    sigma2_w = np.sqrt(vari)\n",
    "    \n",
    "    surge_duration_c=lognorm.rvs(s=sigma2_c,scale=np.exp(mu2_c))\n",
    "    surge_duration_w=lognorm.rvs(s=sigma2_w,scale=np.exp(mu2_w))\n",
    "\n",
    "    waitingtime_c = np.random.exponential(lambda_c)\n",
    "    waitingtime_w = np.random.exponential(lambda_w)\n",
    "    \n",
    "    riseratio_c = np.random.beta(a=alpha_c,b=beta_c)\n",
    "    riseratio_w = np.random.beta(a=alpha_w,b=beta_w)\n",
    "\n",
    "    tau_a_c = waitingtime_c - riseratio_c*surge_duration_c\n",
    "    tau_b_c = tau_a_c + surge_duration_c\n",
    "    tau_a_w = waitingtime_w - riseratio_w*surge_duration_w\n",
    "    tau_b_w = tau_a_w + surge_duration_w\n",
    "\n",
    "    chi_1_c = loc; chi_2_c = ssh; chi_1_w = loc; chi_2_w = ssh\n",
    "    omega1_c = 1/(riseratio_c*surge_duration_c)*np.log(ssh/loc)\n",
    "    omega2_c = -1/(1-riseratio_c)/surge_duration_c*np.log(ssh/loc)\n",
    "    omega1_w= 1/(riseratio_w*surge_duration_w)*np.log(ssh/loc)\n",
    "    omega2_w = -1/(1-riseratio_w)/surge_duration_w*np.log(ssh/loc)\n",
    "\n",
    "    time1_c = np.linspace(tau_a_c,waitingtime_c,nt)\n",
    "    time2_c = np.linspace(waitingtime_c,tau_b_c,nt)\n",
    "    time1_w = np.linspace(tau_a_w,waitingtime_w,nt)\n",
    "    time2_w = np.linspace(waitingtime_w,tau_b_w,nt)\n",
    "\n",
    "    cpi1_c = chi_1_c*np.exp(omega1_c*(time1_c-tau_a_c))\n",
    "    cpi2_c = chi_2_c*np.exp(omega2_c*(time2_c-waitingtime_c))\n",
    "    cpi1_w = chi_1_w*np.exp(omega1_w*(time1_w-tau_a_w))\n",
    "    cpi2_w = chi_2_w*np.exp(omega2_w*(time2_w-waitingtime_w))\n",
    "    \n",
    "    return time1_c,time2_c,time1_w,time2_w,cpi1_c,cpi2_c,cpi1_w,cpi2_w\n",
    "\n",
    "def SurfaceVolFunc(surfaceV,H):\n",
    "    \"\"\"\n",
    "    This function takes the surface volume values calculated using the Surface Volume Processor\n",
    "    and interpolates the values into a surface volume function.\n",
    "    \"\"\"\n",
    "    \n",
    "    a1 = 0;  a2 = 2;  b1 = 1;  b2 = 3;     c1 = 2;  c2 = 4;  d1 = 3;  d2 = 5\n",
    "    e1 = 4;  e2 = 6;  f1 = 5;  f2 = 7;     g1 = 6;  g2 = 8;  h1 = 7;  h2 = 9\n",
    "    i1 = 8;  i2 = 10; j1 = 9;  j2 = 11;    k1 = 10; k2 = 12; l1 = 11; l2 = 13\n",
    "    m1 = 12; m2 = 14; n1 = 13; n2 = 15;    o1 = 14; o2 = 16; p1 = 15; p2 = 17\n",
    "    q1 = 16; q2 = 18; r1 = 17; r2 = 19;    s1 = 18; s2 = 20; t1 = 19; t2 = 21\n",
    "    \n",
    "    popt1, pcov1 = curve_fit(func_fit, surfaceV[a1:a2], H[a1:a2])\n",
    "    popt2, pcov2 = curve_fit(func_fit, surfaceV[b1:b2], H[b1:b2])\n",
    "    popt3, pcov3 = curve_fit(func_fit, surfaceV[c1:c2], H[c1:c2])\n",
    "    popt4, pcov4 = curve_fit(func_fit, surfaceV[d1:d2], H[d1:d2])\n",
    "    popt5, pcov5 = curve_fit(func_fit, surfaceV[e1:e2], H[e1:e2])\n",
    "    popt6, pcov6 = curve_fit(func_fit, surfaceV[f1:f2], H[f1:f2])\n",
    "    popt7, pcov7 = curve_fit(func_fit, surfaceV[g1:g2], H[g1:g2])\n",
    "    popt8, pcov8 = curve_fit(func_fit, surfaceV[h1:h2], H[h1:h2])\n",
    "    popt9, pcov9 = curve_fit(func_fit, surfaceV[i1:i2], H[i1:i2])\n",
    "    popt10, pcov10 = curve_fit(func_fit, surfaceV[j1:j2], H[j1:j2])\n",
    "    popt11, pcov11 = curve_fit(func_fit, surfaceV[k1:k2], H[k1:k2])\n",
    "    popt12, pcov12 = curve_fit(func_fit, surfaceV[l1:l2], H[l1:l2])\n",
    "    popt13, pcov13 = curve_fit(func_fit, surfaceV[m1:m2], H[m1:m2])\n",
    "    popt14, pcov14 = curve_fit(func_fit, surfaceV[n1:n2], H[n1:n2])\n",
    "    popt15, pcov15 = curve_fit(func_fit, surfaceV[o1:o2], H[o1:o2])\n",
    "    popt16, pcov16 = curve_fit(func_fit, surfaceV[p1:p2], H[p1:p2])\n",
    "    popt17, pcov17 = curve_fit(func_fit, surfaceV[q1:q2], H[q1:q2])\n",
    "    popt18, pcov18 = curve_fit(func_fit, surfaceV[r1:r2], H[r1:r2])\n",
    "    popt19, pcov19 = curve_fit(func_fit, surfaceV[s1:s2], H[s1:s2])\n",
    "    popt20, pcov20 = curve_fit(func_fit, surfaceV[t1:t2], H[t1:t2])\n",
    "            \n",
    "    return popt1,popt2,popt3,popt4,popt5,popt6,popt7,popt8,popt9,popt10,popt11,popt12,popt13,popt14,popt15,popt16,popt17,popt18,popt19,popt20\n",
    "\n",
    "def FloodHeight(surfaceV,slope,roughness,SVf1,SVf2,SVf3,SVf4,SVf5,SVf6,SVf7,SVf8,SVf9,SVf10,SVf11,SVf12,SVf13,SVf14,SVf15,SVf16,SVf17,SVf18,SVf19,SVf20,time1_w,time2_w,cpi1_w,cpi2_w,nt,elev,fid,l,ssh,i):\n",
    "    \"\"\"\n",
    "    This function calculates flood heights for a given div without considering flood water redistribution\n",
    "    \"\"\"\n",
    "    sec1_w = time1_w * 60**2; sec2_w = time2_w * 60**2  #time in seconds\n",
    "    h1_w = cpi1_w;            h2_w = cpi2_w\n",
    "    V_new = 0;                v_new = np.zeros(nt*2)\n",
    "    \n",
    "    a1 = 0;  a2 = 2;  b1 = 1;  b2 = 3;     c1 = 2;  c2 = 4;  d1 = 3;  d2 = 5\n",
    "    e1 = 4;  e2 = 6;  f1 = 5;  f2 = 7;     g1 = 6;  g2 = 8;  h1 = 7;  h2 = 9\n",
    "    i1 = 8;  i2 = 10; j1 = 9;  j2 = 11;    k1 = 10; k2 = 12; l1 = 11; l2 = 13\n",
    "    m1 = 12; m2 = 14; n1 = 13; n2 = 15;    o1 = 14; o2 = 16; p1 = 15; p2 = 17\n",
    "    q1 = 16; q2 = 18; r1 = 17; r2 = 19;    s1 = 18; s2 = 20; t1 = 19; t2 = 21\n",
    "\n",
    "    for k in fid:\n",
    "        ds1 = np.tile(sec1_w[1] - sec1_w[0],(nt,1)).transpose() #delta between times in seconds\n",
    "        ds2 = np.tile(sec2_w[1] - sec2_w[0],(nt,1)).transpose()\n",
    "        ds  = np.concatenate((ds1,ds2),axis=1)\n",
    "\n",
    "        h1_new = h1_w - elev[k] #difference betwwen surge heights and elevations\n",
    "        h2_new = h2_w - elev[k]\n",
    "        h_new  = np.concatenate((h1_new,h2_new),axis=0)\n",
    "        v_new[h_new > 0] = ((l * h_new[h_new > 0]) / (l + (2*h_new[h_new > 0])))**(2./3.)*slope**.5 / roughness #Manning's Equation\n",
    "        v_new[h_new <= 0] = 0\n",
    "\n",
    "        V_new = V_new + np.sum(h_new*l*v_new*ds) #volume to plug into surface volume functions to determin flood heights\n",
    "    \n",
    "    #determine flood heights based on surface volume function\n",
    "    fld_h = np.zeros(np.shape(V_new))\n",
    "    fld_h [V_new == 0]= 0\n",
    "    fld_h [(V_new <= surfaceV[a2])  &  (V_new > 0)] = func_fit(V_new[(V_new <= surfaceV[a2])  &  (V_new > 0)],*SVf1)\n",
    "    fld_h [(V_new <= surfaceV[b2])  &  (V_new > surfaceV[a2])] = func_fit(V_new[(V_new <= surfaceV[b2])  &  (V_new > surfaceV[a2])],*SVf2)\n",
    "    fld_h [(V_new <= surfaceV[c2])  &  (V_new > surfaceV[b2])] = func_fit(V_new[(V_new <= surfaceV[c2])  &  (V_new > surfaceV[b2])],*SVf3)\n",
    "    fld_h [(V_new <= surfaceV[d2])  &  (V_new > surfaceV[c2])] = func_fit(V_new[(V_new <= surfaceV[d2])  &  (V_new > surfaceV[c2])],*SVf4)\n",
    "    fld_h [(V_new <= surfaceV[e2])  &  (V_new > surfaceV[d2])] = func_fit(V_new[(V_new <= surfaceV[e2])  &  (V_new > surfaceV[d2])],*SVf5)\n",
    "    fld_h [(V_new <= surfaceV[f2])  &  (V_new > surfaceV[e2])] = func_fit(V_new[(V_new <= surfaceV[f2])  &  (V_new > surfaceV[e2])],*SVf6)\n",
    "    fld_h [(V_new <= surfaceV[g2])  &  (V_new > surfaceV[f2])] = func_fit(V_new[(V_new <= surfaceV[g2])  &  (V_new > surfaceV[f2])],*SVf7)\n",
    "    fld_h [(V_new <= surfaceV[h2])  &  (V_new > surfaceV[f2])] = func_fit(V_new[(V_new <= surfaceV[h2])  &  (V_new > surfaceV[f2])],*SVf8)\n",
    "    fld_h [(V_new <= surfaceV[i2])  &  (V_new > surfaceV[h2])] = func_fit(V_new[(V_new <= surfaceV[i2])  &  (V_new > surfaceV[h2])],*SVf9)\n",
    "    fld_h [(V_new <= surfaceV[j2])  &  (V_new > surfaceV[i2])] = func_fit(V_new[(V_new <= surfaceV[j2])  &  (V_new > surfaceV[i2])],*SVf10)\n",
    "    fld_h [(V_new <= surfaceV[k2])  &  (V_new > surfaceV[j2])] = func_fit(V_new[(V_new <= surfaceV[k2])  &  (V_new > surfaceV[j2])],*SVf11)\n",
    "    fld_h [(V_new <= surfaceV[l2])  &  (V_new > surfaceV[k2])] = func_fit(V_new[(V_new <= surfaceV[l2])  &  (V_new > surfaceV[k2])],*SVf12)\n",
    "    fld_h [(V_new <= surfaceV[m2])  &  (V_new > surfaceV[l2])] = func_fit(V_new[(V_new <= surfaceV[m2])  &  (V_new > surfaceV[l2])],*SVf13)\n",
    "    fld_h [(V_new <= surfaceV[n2])  &  (V_new > surfaceV[m2])] = func_fit(V_new[(V_new <= surfaceV[n2])  &  (V_new > surfaceV[m2])],*SVf14)\n",
    "    fld_h [(V_new <= surfaceV[o2])  &  (V_new > surfaceV[n2])] = func_fit(V_new[(V_new <= surfaceV[o2])  &  (V_new > surfaceV[n2])],*SVf15)\n",
    "    fld_h [(V_new <= surfaceV[p2])  &  (V_new > surfaceV[o2])] = func_fit(V_new[(V_new <= surfaceV[p2])  &  (V_new > surfaceV[o2])],*SVf16)\n",
    "    fld_h [(V_new <= surfaceV[q2])  &  (V_new > surfaceV[p2])] = func_fit(V_new[(V_new <= surfaceV[q2])  &  (V_new > surfaceV[p2])],*SVf17)\n",
    "    fld_h [(V_new <= surfaceV[r2])  &  (V_new > surfaceV[q2])] = func_fit(V_new[(V_new <= surfaceV[r2])  &  (V_new > surfaceV[q2])],*SVf18)\n",
    "    fld_h [(V_new <= surfaceV[s2])  &  (V_new > surfaceV[r2])] = func_fit(V_new[(V_new <= surfaceV[s2])  &  (V_new > surfaceV[r2])],*SVf19)\n",
    "    fld_h [V_new >= surfaceV[s2]] = func_fit(V_new[V_new > surfaceV[s2]],*SVf20)\n",
    "    fld_h [fld_h < 0]= 0\n",
    "    \n",
    "    return fld_h,V_new\n",
    "\n",
    "def FloodHeightWall(surfaceV,ParWall,slope,roughness,SVf1,SVf2,SVf3,SVf4,SVf5,SVf6,SVf7,SVf8,SVf9,SVf10,SVf11,SVf12,SVf13,SVf14,SVf15,SVf16,SVf17,SVf18,SVf19,SVf20,time1_w,time2_w,cpi1_w,cpi2_w,nt,elev,fid,l,ssh,i):\n",
    "    \"\"\"\n",
    "    This function calculates flood heights with an altered Manning's Equation to account for the height of the wall.\n",
    "    \"\"\"\n",
    "    \n",
    "    sec1_w = time1_w * 60**2; sec2_w = time2_w * 60**2\n",
    "    h1_w = cpi1_w;            h2_w = cpi2_w\n",
    "    V_new = 0; v_new = np.zeros(nt*2)\n",
    "    \n",
    "    a1 = 0;  a2 = 2;  b1 = 1;  b2 = 3;     c1 = 2;  c2 = 4;  d1 = 3;  d2 = 5\n",
    "    e1 = 4;  e2 = 6;  f1 = 5;  f2 = 7;     g1 = 6;  g2 = 8;  h1 = 7;  h2 = 9\n",
    "    i1 = 8;  i2 = 10; j1 = 9;  j2 = 11;    k1 = 10; k2 = 12; l1 = 11; l2 = 13\n",
    "    m1 = 12; m2 = 14; n1 = 13; n2 = 15;    o1 = 14; o2 = 16; p1 = 15; p2 = 17\n",
    "    q1 = 16; q2 = 18; r1 = 17; r2 = 19;    s1 = 18; s2 = 20; t1 = 19; t2 = 21\n",
    "\n",
    "    for k in fid:\n",
    "        if k in ParWall:\n",
    "            ds1 = np.tile(sec1_w[1] - sec1_w[0],(nt,1)).transpose()\n",
    "            ds2 = np.tile(sec2_w[1] - sec2_w[0],(nt,1)).transpose()\n",
    "            ds  = np.concatenate((ds1,ds2),axis=1)\n",
    "\n",
    "            h1_new = h1_w - elev[k]\n",
    "            h2_new = h2_w - elev[k]\n",
    "            h_new  = np.concatenate((h1_new,h2_new),axis=0)\n",
    "            cwr = 0.611 + 0.075*h_new[h_new > 0]/elev[k]            \n",
    "            v_new[h_new > 0] = ((l * h_new[h_new > 0]) / (l + (2*h_new[h_new > 0])))**(2./3.)*slope**.5 / roughness * cwr\n",
    "            v_new[h_new <= 0] = 0\n",
    "            V_new = V_new + np.sum(h_new*l*v_new*ds)\n",
    "            \n",
    "        else:\n",
    "            ds1 = np.tile(sec1_w[1] - sec1_w[0],(nt,1)).transpose()\n",
    "            ds2 = np.tile(sec2_w[1] - sec2_w[0],(nt,1)).transpose()\n",
    "            ds  = np.concatenate((ds1,ds2),axis=1)\n",
    "\n",
    "            h1_new = h1_w - elev[k]\n",
    "            h2_new = h2_w - elev[k]\n",
    "            h_new  = np.concatenate((h1_new,h2_new),axis=0)\n",
    "            v_new[h_new > 0] = ((l * h_new[h_new > 0]) / (l + (2*h_new[h_new > 0])))**(2./3.)*slope**.5 / roughness\n",
    "            v_new[h_new <= 0] = 0\n",
    "            V_new = V_new + np.sum(h_new*l*v_new*ds)\n",
    "        \n",
    "    fld_h = np.zeros(np.shape(V_new))\n",
    "    fld_h [V_new == 0]= 0\n",
    "    fld_h [(V_new <= surfaceV[a2])  &  (V_new > 0)] = func_fit(V_new[(V_new <= surfaceV[a2])  &  (V_new > 0)],*SVf1)\n",
    "    fld_h [(V_new <= surfaceV[b2])  &  (V_new > surfaceV[a2])] = func_fit(V_new[(V_new <= surfaceV[b2])  &  (V_new > surfaceV[a2])],*SVf2)\n",
    "    fld_h [(V_new <= surfaceV[c2])  &  (V_new > surfaceV[b2])] = func_fit(V_new[(V_new <= surfaceV[c2])  &  (V_new > surfaceV[b2])],*SVf3)\n",
    "    fld_h [(V_new <= surfaceV[d2])  &  (V_new > surfaceV[c2])] = func_fit(V_new[(V_new <= surfaceV[d2])  &  (V_new > surfaceV[c2])],*SVf4)\n",
    "    fld_h [(V_new <= surfaceV[e2])  &  (V_new > surfaceV[d2])] = func_fit(V_new[(V_new <= surfaceV[e2])  &  (V_new > surfaceV[d2])],*SVf5)\n",
    "    fld_h [(V_new <= surfaceV[f2])  &  (V_new > surfaceV[e2])] = func_fit(V_new[(V_new <= surfaceV[f2])  &  (V_new > surfaceV[e2])],*SVf6)\n",
    "    fld_h [(V_new <= surfaceV[g2])  &  (V_new > surfaceV[f2])] = func_fit(V_new[(V_new <= surfaceV[g2])  &  (V_new > surfaceV[f2])],*SVf7)\n",
    "    fld_h [(V_new <= surfaceV[h2])  &  (V_new > surfaceV[f2])] = func_fit(V_new[(V_new <= surfaceV[h2])  &  (V_new > surfaceV[f2])],*SVf8)\n",
    "    fld_h [(V_new <= surfaceV[i2])  &  (V_new > surfaceV[h2])] = func_fit(V_new[(V_new <= surfaceV[i2])  &  (V_new > surfaceV[h2])],*SVf9)\n",
    "    fld_h [(V_new <= surfaceV[j2])  &  (V_new > surfaceV[i2])] = func_fit(V_new[(V_new <= surfaceV[j2])  &  (V_new > surfaceV[i2])],*SVf10)\n",
    "    fld_h [(V_new <= surfaceV[k2])  &  (V_new > surfaceV[j2])] = func_fit(V_new[(V_new <= surfaceV[k2])  &  (V_new > surfaceV[j2])],*SVf11)\n",
    "    fld_h [(V_new <= surfaceV[l2])  &  (V_new > surfaceV[k2])] = func_fit(V_new[(V_new <= surfaceV[l2])  &  (V_new > surfaceV[k2])],*SVf12)\n",
    "    fld_h [(V_new <= surfaceV[m2])  &  (V_new > surfaceV[l2])] = func_fit(V_new[(V_new <= surfaceV[m2])  &  (V_new > surfaceV[l2])],*SVf13)\n",
    "    fld_h [(V_new <= surfaceV[n2])  &  (V_new > surfaceV[m2])] = func_fit(V_new[(V_new <= surfaceV[n2])  &  (V_new > surfaceV[m2])],*SVf14)\n",
    "    fld_h [(V_new <= surfaceV[o2])  &  (V_new > surfaceV[n2])] = func_fit(V_new[(V_new <= surfaceV[o2])  &  (V_new > surfaceV[n2])],*SVf15)\n",
    "    fld_h [(V_new <= surfaceV[p2])  &  (V_new > surfaceV[o2])] = func_fit(V_new[(V_new <= surfaceV[p2])  &  (V_new > surfaceV[o2])],*SVf16)\n",
    "    fld_h [(V_new <= surfaceV[q2])  &  (V_new > surfaceV[p2])] = func_fit(V_new[(V_new <= surfaceV[q2])  &  (V_new > surfaceV[p2])],*SVf17)\n",
    "    fld_h [(V_new <= surfaceV[r2])  &  (V_new > surfaceV[q2])] = func_fit(V_new[(V_new <= surfaceV[r2])  &  (V_new > surfaceV[q2])],*SVf18)\n",
    "    fld_h [(V_new <= surfaceV[s2])  &  (V_new > surfaceV[r2])] = func_fit(V_new[(V_new <= surfaceV[s2])  &  (V_new > surfaceV[r2])],*SVf19)\n",
    "    fld_h [V_new >= surfaceV[s2]] = func_fit(V_new[V_new > surfaceV[s2]],*SVf20)\n",
    "    fld_h [fld_h < 0]= 0\n",
    "             \n",
    "    return fld_h,V_new\n",
    "\n",
    "def FloodTravel(ssh,V_new,SVf1,SVf2,SVf3,SVf4,SVf5,SVf6,SVf7,SVf8,SVf9,SVf10,SVf11,SVf12,SVf13,SVf14,SVf15,SVf16,SVf17,SVf18,SVf19,SVf20):\n",
    "    \"\"\"\n",
    "    This function takes the volumes redistributed in FloodTravelSectGroup and calculates the new flood heights for each div based on those volumes.\n",
    "    \"\"\"\n",
    "    \n",
    "    a1 = 0;  a2 = 2;  b1 = 1;  b2 = 3;     c1 = 2;  c2 = 4;  d1 = 3;  d2 = 5\n",
    "    e1 = 4;  e2 = 6;  f1 = 5;  f2 = 7;     g1 = 6;  g2 = 8;  h1 = 7;  h2 = 9\n",
    "    i1 = 8;  i2 = 10; j1 = 9;  j2 = 11;    k1 = 10; k2 = 12; l1 = 11; l2 = 13\n",
    "    m1 = 12; m2 = 14; n1 = 13; n2 = 15;    o1 = 14; o2 = 16; p1 = 15; p2 = 17\n",
    "    q1 = 16; q2 = 18; r1 = 17; r2 = 19;    s1 = 18; s2 = 20; t1 = 19; t2 = 21\n",
    "    \n",
    "    fld_h = np.zeros(np.shape(V_new))\n",
    "    fld_h [V_new == 0]= 0\n",
    "    fld_h [(V_new <= surfaceV[a2])  &  (V_new > 0)] = func_fit(V_new[(V_new <= surfaceV[a2])  &  (V_new > 0)],*SVf1)\n",
    "    fld_h [(V_new <= surfaceV[b2])  &  (V_new > surfaceV[a2])] = func_fit(V_new[(V_new <= surfaceV[b2])  &  (V_new > surfaceV[a2])],*SVf2)\n",
    "    fld_h [(V_new <= surfaceV[c2])  &  (V_new > surfaceV[b2])] = func_fit(V_new[(V_new <= surfaceV[c2])  &  (V_new > surfaceV[b2])],*SVf3)\n",
    "    fld_h [(V_new <= surfaceV[d2])  &  (V_new > surfaceV[c2])] = func_fit(V_new[(V_new <= surfaceV[d2])  &  (V_new > surfaceV[c2])],*SVf4)\n",
    "    fld_h [(V_new <= surfaceV[e2])  &  (V_new > surfaceV[d2])] = func_fit(V_new[(V_new <= surfaceV[e2])  &  (V_new > surfaceV[d2])],*SVf5)\n",
    "    fld_h [(V_new <= surfaceV[f2])  &  (V_new > surfaceV[e2])] = func_fit(V_new[(V_new <= surfaceV[f2])  &  (V_new > surfaceV[e2])],*SVf6)\n",
    "    fld_h [(V_new <= surfaceV[g2])  &  (V_new > surfaceV[f2])] = func_fit(V_new[(V_new <= surfaceV[g2])  &  (V_new > surfaceV[f2])],*SVf7)\n",
    "    fld_h [(V_new <= surfaceV[h2])  &  (V_new > surfaceV[f2])] = func_fit(V_new[(V_new <= surfaceV[h2])  &  (V_new > surfaceV[f2])],*SVf8)\n",
    "    fld_h [(V_new <= surfaceV[i2])  &  (V_new > surfaceV[h2])] = func_fit(V_new[(V_new <= surfaceV[i2])  &  (V_new > surfaceV[h2])],*SVf9)\n",
    "    fld_h [(V_new <= surfaceV[j2])  &  (V_new > surfaceV[i2])] = func_fit(V_new[(V_new <= surfaceV[j2])  &  (V_new > surfaceV[i2])],*SVf10)\n",
    "    fld_h [(V_new <= surfaceV[k2])  &  (V_new > surfaceV[j2])] = func_fit(V_new[(V_new <= surfaceV[k2])  &  (V_new > surfaceV[j2])],*SVf11)\n",
    "    fld_h [(V_new <= surfaceV[l2])  &  (V_new > surfaceV[k2])] = func_fit(V_new[(V_new <= surfaceV[l2])  &  (V_new > surfaceV[k2])],*SVf12)\n",
    "    fld_h [(V_new <= surfaceV[m2])  &  (V_new > surfaceV[l2])] = func_fit(V_new[(V_new <= surfaceV[m2])  &  (V_new > surfaceV[l2])],*SVf13)\n",
    "    fld_h [(V_new <= surfaceV[n2])  &  (V_new > surfaceV[m2])] = func_fit(V_new[(V_new <= surfaceV[n2])  &  (V_new > surfaceV[m2])],*SVf14)\n",
    "    fld_h [(V_new <= surfaceV[o2])  &  (V_new > surfaceV[n2])] = func_fit(V_new[(V_new <= surfaceV[o2])  &  (V_new > surfaceV[n2])],*SVf15)\n",
    "    fld_h [(V_new <= surfaceV[p2])  &  (V_new > surfaceV[o2])] = func_fit(V_new[(V_new <= surfaceV[p2])  &  (V_new > surfaceV[o2])],*SVf16)\n",
    "    fld_h [(V_new <= surfaceV[q2])  &  (V_new > surfaceV[p2])] = func_fit(V_new[(V_new <= surfaceV[q2])  &  (V_new > surfaceV[p2])],*SVf17)\n",
    "    fld_h [(V_new <= surfaceV[r2])  &  (V_new > surfaceV[q2])] = func_fit(V_new[(V_new <= surfaceV[r2])  &  (V_new > surfaceV[q2])],*SVf18)\n",
    "    fld_h [(V_new <= surfaceV[s2])  &  (V_new > surfaceV[r2])] = func_fit(V_new[(V_new <= surfaceV[s2])  &  (V_new > surfaceV[r2])],*SVf19)\n",
    "    fld_h [V_new >= surfaceV[s2]] = func_fit(V_new[V_new > surfaceV[s2]],*SVf20)\n",
    "    fld_h [fld_h < 0]= 0\n",
    "    if fld_h > ssh:\n",
    "        fld_h = ssh\n",
    "    \n",
    "    return fld_h\n",
    "\n",
    "def FloodTravelSectGroup(ssh,sect0,sect1,sect2,sect3,sect_3,sect_2,sect_1,files,V,\n",
    "                         SVfg1,SVfg2,SVfg3,SVfg4,SVfg5,SVfg6,SVfg7,SVfg8,SVfg9,SVfg10,SVfg11,SVfg12,SVfg13,SVfg14,SVfg15,SVfg16,SVfg17,SVfg18,SVfg19,SVfg20):\n",
    "    \"\"\"\n",
    "    This function propagates the flooding to different divisions within a group and calls FloodTravel to resditribute the water\n",
    "    \"\"\"\n",
    "    V_sect = np.zeros(ndiv18)\n",
    "    fld_h_sect = np.zeros(ndiv18)\n",
    "    V_sect_avr = np.zeros(ndiv18)\n",
    " \n",
    "    # ID3~14\n",
    "    for i in range(12):\n",
    "        m = int(sect3[i][3])\n",
    "        for s in sect3[i]:\n",
    "            s = int(s)\n",
    "            V_sect[m] = V_sect[m] + V[s]\n",
    "        V_sect_avr[m] = V_sect[m]/7\n",
    "\n",
    "    #ID0\n",
    "    i = 10; \n",
    "    for s in sect0:\n",
    "        s = int(s)\n",
    "        V_sect[i] = V_sect[i] + V[s]\n",
    "    V_sect_avr[i] = V_sect[i]/4\n",
    "\n",
    "    #ID1\n",
    "    i = 11; \n",
    "    for s in sect1:\n",
    "        s = int(s)\n",
    "        V_sect[i] = V_sect[i] + V[s]\n",
    "    V_sect_avr[i] = V_sect[i]/5\n",
    "\n",
    "    #ID2\n",
    "    i = 5;\n",
    "    for s in sect2:\n",
    "        s = int(s)\n",
    "        V_sect[i] = V_sect[i] + V[s]\n",
    "    V_sect_avr[i] = V_sect[i]/6\n",
    "\n",
    "    #ID15\n",
    "    i = 7;\n",
    "    for s in sect_3:\n",
    "        s = int(s)\n",
    "        V_sect[i] = V_sect[i] + V[s]\n",
    "    V_sect_avr[i] = V_sect[i]/6\n",
    "\n",
    "    #ID16\n",
    "    i = 2;\n",
    "    for s in sect_2:\n",
    "        s = int(s)\n",
    "        V_sect[i] = V_sect[i] + V[s]\n",
    "    V_sect_avr[i] = V_sect[i]/5\n",
    "\n",
    "    #ID17\n",
    "    i = 6;\n",
    "    for s in sect_1:\n",
    "        s = int(s)\n",
    "        V_sect[i] = V_sect[i] + V[s]\n",
    "    V_sect_avr[i] = V_sect[i]/4\n",
    "\n",
    "    for i in range(ndiv18):    \n",
    "        # redistribute the water\n",
    "        fld_h_sect[i] = FloodTravel(ssh,V_sect[i],SVfg1[i,:],SVfg2[i,:],SVfg3[i,:],SVfg4[i,:],SVfg5[i,:],SVfg6[i,:],SVfg7[i,:],SVfg8[i,:],SVfg9[i,:],SVfg10[i,:],SVfg11[i,:],SVfg12[i,:],SVfg13[i,:],SVfg14[i,:],SVfg15[i,:],SVfg16[i,:],SVfg17[i,:],SVfg18[i,:],SVfg19[i,:],SVfg20[i,:])\n",
    "\n",
    "    fld_h_sect[fld_h_sect<0] = 0\n",
    "    \n",
    "    return fld_h_sect,V_sect_avr\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-30T16:59:54.513406Z",
     "start_time": "2020-11-30T16:59:48.362498Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Mimi\\anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n"
     ]
    }
   ],
   "source": [
    "# Storm Setup Values\n",
    "loc = 0.68;        ftm = 0.3048;       nt = 100\n",
    "xi_c = 0.08;       delta_c = 0.14;     lambda_c = 4.11\n",
    "xi_w = 0.45;       delta_w = 0.13;     lambda_w = 1.52\n",
    "beta1_1_c = 27.84; beta1_2_c = -18.58; beta2_1_c = 10.46; beta2_2_c = -6.98\n",
    "beta1_1_w = 19.23; beta1_2_w = -11.90; beta2_1_w = 4.96;  beta2_2_w = -1.89\n",
    "alpha_c = 3;       beta_c = 3.1;       alpha_w = 2.56;    beta_w = 2.07\n",
    "g = 9.81 # m/s2    \n",
    "nt = 10\n",
    "\n",
    "# Topography Data\n",
    "topo   = pd.read_csv(r\"Data\\LMN_div_low19.csv\")\n",
    "elev   = topo[\"MEAN2\"]\n",
    "fid    = topo[\"FID\"]\n",
    "div18  = topo[\"DIV18\"]\n",
    "ndiv18 = np.unique(div18).size\n",
    "\n",
    "roughness = pd.read_csv(r\"Data\\LMN_Roughness.csv\")[\"Roughness\"]\n",
    "slope     = pd.read_csv(r\"Data\\LMN_Slope.csv\")[\"Slope\"]\n",
    "\n",
    "# Import Surface Volume\n",
    "files = glob.glob(r'Data\\NewSurfaceVolumeCombined\\LMN_div18_new_*.csv') \n",
    "groupcsvfiles = glob.glob(r'Data\\NewSurfaceVolumeGrouped\\LMN_div18_new_grouped_*.csv') \n",
    "\n",
    "H = np.append(np.linspace(0,3,13),np.linspace(3.5,7,8))\n",
    "\n",
    "l        = 100  # segment unit length\n",
    "sections = [10,11,5,12,1,13,16,4,17,15,3,14,0,9,8,7,2,6] #div number, in sequential geographic order from east: north to south, then west: north to south\n",
    "\n",
    "# Setup - Division Connectivities\n",
    "# these sections give the groupings of different divisions when propagating flood water\n",
    "sect3 = np.zeros([len(sections)-6,7])\n",
    "k = 3\n",
    "for i in sections[3:-3]:\n",
    "    sect3[k-3] = [sections[k-3],sections[k-2],sections[k-1],i,sections[k+1],sections[k+2],sections[k+3]]\n",
    "    k = k + 1\n",
    "\n",
    "sect0  = [sections[0],sections[1],sections[2],sections[3]] \n",
    "sect1  = [sections[0],sections[1],sections[2],sections[3],sections[4]]\n",
    "sect2  = [sections[0],sections[1],sections[2],sections[3],sections[4],sections[5]]\n",
    "sect_3 = [sections[-6],sections[-5],sections[-4],sections[-3],sections[-2],sections[-1]]\n",
    "sect_2 = [sections[-5],sections[-4],sections[-3],sections[-2],sections[-1]]\n",
    "sect_1 = [sections[-4],sections[-3],sections[-2],sections[-1]]\n",
    "\n",
    "#Surface Volume functions for grouped divisions\n",
    "SVfg1 = np.zeros([ndiv18,2]);  SVfg2 = np.zeros([ndiv18,2]);  SVfg3 = np.zeros([ndiv18,2]);  SVfg4 = np.zeros([ndiv18,2])\n",
    "SVfg5 = np.zeros([ndiv18,2]);  SVfg6 = np.zeros([ndiv18,2]);  SVfg7 = np.zeros([ndiv18,2]);  SVfg8 = np.zeros([ndiv18,2])\n",
    "SVfg9 = np.zeros([ndiv18,2]);  SVfg10 = np.zeros([ndiv18,2]); SVfg11 = np.zeros([ndiv18,2]); SVfg12 = np.zeros([ndiv18,2])\n",
    "SVfg13 = np.zeros([ndiv18,2]); SVfg14 = np.zeros([ndiv18,2]); SVfg15 = np.zeros([ndiv18,2]); SVfg16 = np.zeros([ndiv18,2])\n",
    "SVfg17 = np.zeros([ndiv18,2]); SVfg18 = np.zeros([ndiv18,2]); SVfg19 = np.zeros([ndiv18,2]); SVfg20= np.zeros([ndiv18,2])\n",
    "\n",
    "i = 0\n",
    "\n",
    "for f in groupcsvfiles:\n",
    "    surfaceVg = pd.read_csv(f)[\"volume\"]\n",
    "    SVfg1[i,:],SVfg2[i,:],SVfg3[i,:],SVfg4[i,:],SVfg5[i,:],SVfg6[i,:],SVfg7[i,:],SVfg8[i,:],SVfg9[i,:],SVfg10[i,:],SVfg11[i,:],SVfg12[i,:],SVfg13[i,:],SVfg14[i,:],SVfg15[i,:],SVfg16[i,:],SVfg17[i,:],SVfg18[i,:],SVfg19[i,:],SVfg20[i,:] = SurfaceVolFunc(surfaceVg,H)\n",
    "    i = i + 1 \n",
    "\n",
    "#Surface Volume functions for individual divisions\n",
    "SVf1 = np.zeros([ndiv18,2]);  SVf2 = np.zeros([ndiv18,2]);  SVf3 = np.zeros([ndiv18,2]);  SVf4 = np.zeros([ndiv18,2])\n",
    "SVf5 = np.zeros([ndiv18,2]);  SVf6 = np.zeros([ndiv18,2]);  SVf7 = np.zeros([ndiv18,2]);  SVf8 = np.zeros([ndiv18,2])\n",
    "SVf9 = np.zeros([ndiv18,2]);  SVf10 = np.zeros([ndiv18,2]); SVf11 = np.zeros([ndiv18,2]); SVf12 = np.zeros([ndiv18,2])\n",
    "SVf13 = np.zeros([ndiv18,2]); SVf14 = np.zeros([ndiv18,2]); SVf15 = np.zeros([ndiv18,2]); SVf16 = np.zeros([ndiv18,2])\n",
    "SVf17 = np.zeros([ndiv18,2]); SVf18 = np.zeros([ndiv18,2]); SVf19 = np.zeros([ndiv18,2]); SVf20= np.zeros([ndiv18,2])\n",
    "\n",
    "SV_all = []\n",
    "\n",
    "i = 0\n",
    "for f in files:\n",
    "    surfaceV_height = pd.read_csv(f)\n",
    "    surfaceV = surfaceV_height[\"volume\"]\n",
    "    SV_all = np.append(SV_all,surfaceV)\n",
    "    SVf1[i,:],SVf2[i,:],SVf3[i,:],SVf4[i,:],SVf5[i,:],SVf6[i,:],SVf7[i,:],SVf8[i,:],SVf9[i,:],SVf10[i,:],SVf11[i,:],SVf12[i,:],SVf13[i,:],SVf14[i,:],SVf15[i,:],SVf16[i,:],SVf17[i,:],SVf18[i,:],SVf19[i,:],SVf20[i,:] = SurfaceVolFunc(surfaceV,H)\n",
    "    i = i+1\n",
    "    \n",
    "SV_all = SV_all.reshape(18,21)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flood Simulation - No Protective Measures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flood simulation takes in storm peaks, surge heights, and time, as well as information about the terrain, and the surface volume functions from the Import Data cell to obtain the floodheights of individual cells and redistribute the flooding to adjacent cells. The resulting flood heights of the storm for cold and warm storms are outputted in a csv.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-30T16:59:54.899978Z",
     "start_time": "2020-11-30T16:59:54.515008Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed Time:  0.3480689525604248 sec\n",
      "Complete!\n"
     ]
    }
   ],
   "source": [
    "s_i = 0 # hurricane data \n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "# Storm Surge imports\n",
    "peak_c = ssh\n",
    "peak_w = ssh\n",
    "[time1_c,time2_c,time1_w,time2_w,cpi1_c,cpi2_c,cpi1_w,cpi2_w]=FloodParam(ssh,nt,loc,beta1_1_c,beta1_2_c,beta2_1_c,beta2_2_c,beta1_1_w,beta1_2_w,beta2_1_w,beta2_2_w,alpha_c,beta_c,alpha_w,beta_w)\n",
    "\n",
    "fld_h_c = np.zeros(ndiv18);  fld_h_w = np.zeros(ndiv18)\n",
    "V_c     = np.zeros(ndiv18);  V_w     = np.zeros(ndiv18)\n",
    "\n",
    "#Function calls\n",
    "\n",
    "#Flood height calculation\n",
    "i = 0\n",
    "for i in range(18):\n",
    "    fld_h_w[i],V_w[i] = FloodHeight(SV_all[i],slope[i],roughness[i],SVf1[i,:],SVf2[i,:],SVf3[i,:],SVf4[i,:],SVf5[i,:],SVf6[i,:],SVf7[i,:],SVf8[i,:],SVf9[i,:],SVf10[i,:],SVf11[i,:],SVf12[i,:],SVf13[i,:],SVf14[i,:],SVf15[i,:],SVf16[i,:],SVf17[i,:],SVf18[i,:],SVf19[i,:],SVf20[i,:],time1_w,time2_w,cpi1_w,cpi2_w,nt,elev[div18==i],fid[div18==i],l,peak_w,i)\n",
    "    fld_h_c[i],V_c[i] = FloodHeight(SV_all[i],slope[i],roughness[i],SVf1[i,:],SVf2[i,:],SVf3[i,:],SVf4[i,:],SVf5[i,:],SVf6[i,:],SVf7[i,:],SVf8[i,:],SVf9[i,:],SVf10[i,:],SVf11[i,:],SVf12[i,:],SVf13[i,:],SVf14[i,:],SVf15[i,:],SVf16[i,:],SVf17[i,:],SVf18[i,:],SVf19[i,:],SVf20[i,:],time1_c,time2_c,cpi1_c,cpi2_c,nt,elev[div18==i],fid[div18==i],l,peak_c,i)\n",
    "\n",
    "# redistribution - group\n",
    "fld_h_w_sect_g,V_w_sect_avr = FloodTravelSectGroup(peak_w,sect0,sect1,sect2,sect3,sect_3,sect_2,sect_1,files,V_w,SVfg1,SVfg2,SVfg3,SVfg4,SVfg5,SVfg6,SVfg7,SVfg8,SVfg9,SVfg10,SVfg11,SVfg12,SVfg13,SVfg14,SVfg15,SVfg16,SVfg17,SVfg18,SVfg19,SVfg20)\n",
    "fld_h_c_sect_g,V_c_sect_avr = FloodTravelSectGroup(peak_c,sect0,sect1,sect2,sect3,sect_3,sect_2,sect_1,files,V_c,SVfg1,SVfg2,SVfg3,SVfg4,SVfg5,SVfg6,SVfg7,SVfg8,SVfg9,SVfg10,SVfg11,SVfg12,SVfg13,SVfg14,SVfg15,SVfg16,SVfg17,SVfg18,SVfg19,SVfg20)\n",
    "\n",
    "with open('Sample_Output_Data/single_storm_flood_heights_w_sandy.csv','w') as out:\n",
    "    write = csv.writer(out)\n",
    "    write.writerow(fld_h_w_sect_g)\n",
    "with open('Sample_Output_Data/single_storm_flood_heights_c_sandy.csv','w') as out:\n",
    "    write = csv.writer(out)\n",
    "    write.writerow(fld_h_c_sect_g)\n",
    "\n",
    "elapsed_time = time.time() - start_time\n",
    "print(\"Elapsed Time: \",elapsed_time, \"sec\")            \n",
    "print(\"Complete!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flood Simulation - BigU Lower East Side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This cell takes in storm parameters and surface volumes as well as parameters for a sea wall called the BigU and calculates the flood heights given with the effects of the BigU. The resulting output are csv files for warm and cold storm flood heights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-11-30T16:59:55.107425Z",
     "start_time": "2020-11-30T16:59:54.901974Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed Time:  0.182511568069458 sec\n",
      "Complete!\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-7-0aa8d48e0337>:10: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  elev[BigU_les] = elev_wall[BigU_les] + wallh_bigU\n"
     ]
    }
   ],
   "source": [
    "s_i = 0 # hurricane data \n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "topo = pd.read_csv(r\"Data\\LMN_div_low19.csv\")\n",
    "elev = topo[\"MEAN2\"]\n",
    "elev_wall = topo[\"MEAN3\"] #MEAN3 is the elevation at the coast, where the wall will be placed\n",
    "BigU_les = pd.read_csv(r\"Data\\BigU_LES_all.csv\")[\"ID\"] #this csv specifies the locations of the wall\n",
    "wallh_bigU = 15 * ftm #the \"15\" specifies the height of the wall at 15 ft\n",
    "elev[BigU_les] = elev_wall[BigU_les] + wallh_bigU\n",
    "\n",
    "# Storm Surge imports\n",
    "peak_c = ssh\n",
    "peak_w = ssh\n",
    "[time1_c,time2_c,time1_w,time2_w,cpi1_c,cpi2_c,cpi1_w,cpi2_w]=FloodParam(ssh,nt,loc,beta1_1_c,beta1_2_c,beta2_1_c,beta2_2_c,beta1_1_w,beta1_2_w,beta2_1_w,beta2_2_w,alpha_c,beta_c,alpha_w,beta_w)\n",
    "\n",
    "fld_h_c = np.zeros(ndiv18);  fld_h_w = np.zeros(ndiv18)\n",
    "V_c     = np.zeros(ndiv18);  V_w     = np.zeros(ndiv18)\n",
    "\n",
    "#Function Calls\n",
    "\n",
    "#Flood height calculation with walls\n",
    "i = 0\n",
    "for i in range(18):\n",
    "    fld_h_w[i],V_w[i] = FloodHeightWall(SV_all[i],BigU_les,slope[i],roughness[i],SVf1[i,:],SVf2[i,:],SVf3[i,:],SVf4[i,:],SVf5[i,:],SVf6[i,:],SVf7[i,:],SVf8[i,:],SVf9[i,:],SVf10[i,:],SVf11[i,:],SVf12[i,:],SVf13[i,:],SVf14[i,:],SVf15[i,:],SVf16[i,:],SVf17[i,:],SVf18[i,:],SVf19[i,:],SVf20[i,:],time1_w,time2_w,cpi1_w,cpi2_w,nt,elev[div18==i],fid[div18==i],l,peak_w,i)\n",
    "    fld_h_c[i],V_c[i] = FloodHeightWall(SV_all[i],BigU_les,slope[i],roughness[i],SVf1[i,:],SVf2[i,:],SVf3[i,:],SVf4[i,:],SVf5[i,:],SVf6[i,:],SVf7[i,:],SVf8[i,:],SVf9[i,:],SVf10[i,:],SVf11[i,:],SVf12[i,:],SVf13[i,:],SVf14[i,:],SVf15[i,:],SVf16[i,:],SVf17[i,:],SVf18[i,:],SVf19[i,:],SVf20[i,:],time1_c,time2_c,cpi1_c,cpi2_c,nt,elev[div18==i],fid[div18==i],l,peak_c,i)\n",
    "\n",
    "# redistribution - group\n",
    "fld_h_w_sect_g,V_w_sect_avr = FloodTravelSectGroup(peak_w,sect0,sect1,sect2,sect3,sect_3,sect_2,sect_1,files,V_w,SVfg1,SVfg2,SVfg3,SVfg4,SVfg5,SVfg6,SVfg7,SVfg8,SVfg9,SVfg10,SVfg11,SVfg12,SVfg13,SVfg14,SVfg15,SVfg16,SVfg17,SVfg18,SVfg19,SVfg20)\n",
    "fld_h_c_sect_g,V_c_sect_avr = FloodTravelSectGroup(peak_c,sect0,sect1,sect2,sect3,sect_3,sect_2,sect_1,files,V_c,SVfg1,SVfg2,SVfg3,SVfg4,SVfg5,SVfg6,SVfg7,SVfg8,SVfg9,SVfg10,SVfg11,SVfg12,SVfg13,SVfg14,SVfg15,SVfg16,SVfg17,SVfg18,SVfg19,SVfg20)\n",
    "\n",
    "with open('Sample_Output_Data/single_storm_flood_heights_bigU_w_sandy.csv','w') as out:\n",
    "    write = csv.writer(out)\n",
    "    write.writerow(fld_h_w_sect_g)\n",
    "with open('Sample_Output_Data/single_storm_flood_heights_bigU_c_sandy.csv','w') as out:\n",
    "    write = csv.writer(out)\n",
    "    write.writerow(fld_h_c_sect_g)\n",
    "\n",
    "elapsed_time = time.time() - start_time\n",
    "print(\"Elapsed Time: \",elapsed_time, \"sec\")            \n",
    "print(\"Complete!\")\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "289.302px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "390.179px",
    "left": "177.581px",
    "right": "20px",
    "top": "135.979px",
    "width": "351.819px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}