Location_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from scipy import stats\n",
    "from scipy import interpolate\n",
    "import matplotlib.cm as cm\n",
    "import numpy.polynomial.polynomial as poly\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('../WORK/surf_bT_ret_171031.txt', header=None, delimiter=r\"\\s+\", names = [' ', 'UTC', 'MY', 'L_s', 'LTST', 'Surf_lat', 'Surf_lon', 'Surf_rad', 'Surf_elev', 'T_surf', 'T_surf_err', 'Dust_column', 'Dust_column_err', 'H2Oice_column', 'H2Oice_column_err', 'p_surf', 'p_surf_err', 'P_qual', 'T_qual'])\n",
    "#read data into dataframe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above reads the text file into a dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [
     11
    ]
   },
   "outputs": [],
   "source": [
    "########\n",
    "# pco2 #\n",
    "########\n",
    "# --------------------------------------------\n",
    "# Equilibrium vapor pressure over solid CO2\n",
    "# Brown and Ziegler (1980)\n",
    "# --------------------------------------------\n",
    "# Input:\n",
    "#    T = temperature of solid [K]\n",
    "# Output:\n",
    "#    vapor pressure [Pa]\n",
    "def pco2(T):\n",
    "    A0 = 2.13807649e1\n",
    "    A1 = -2.57064700e3\n",
    "    A2 = -7.78129489e4\n",
    "    A3 = 4.32506256e6\n",
    "    A4 = -1.20671368e8\n",
    "    A5 = 1.34966306e9\n",
    "    # Pressure in torr\n",
    "    ptorr = np.exp(A0 + (A1/(T)) + (A2/(T**2)) + (A3/(T**3)) + (A4/(T**4)) + (A5/(T**5)))\n",
    "    # Pressure in Pa\n",
    "    p = ptorr*133.3223684211\n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": [
     12
    ]
   },
   "outputs": [],
   "source": [
    "########\n",
    "# tco2 #\n",
    "########\n",
    "#-------------------------------------------------------------------------\n",
    "#Calculate (approximately) the condensation temperature of CO2 at the\n",
    "#specified pressure(s) (Pa), using the empirical law from Brown and \n",
    "#Ziegler (1980). This function calls the pco2() routine.\n",
    "#-------------------------------------------------------------------------\n",
    "# Input:\n",
    "#    p = vapor pressure [Pa]\n",
    "# Output:\n",
    "#    temperature of solid [K]\n",
    "def tco2(p):\n",
    "    Trange = np.arange(30,350)\n",
    "    T = interpolate.pchip_interpolate(pco2(Trange),Trange,p)\n",
    "    return T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both Formulas above are to calculate the equilibrium vapor pressure and temperature for carbon dioxide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Function: ls2sol\n",
    "# Purpose: convert Mars \"Lsubs\" to sol (day of year)\n",
    "# Input: \n",
    "#    ls = areocentric longitude of the Sun [array]\n",
    "#    n = degree of polynomial fit [scalar]\n",
    "def ls2sol(ls, n):\n",
    "    # data\n",
    "    lsdata = (0,30,60,90,120,150,180,210,240,270,300,330,360) # L_s array\n",
    "    soldata = (0,61.2,126.6,193.3,257.8,317.5,371.9,421.6,468.5,514.6,562.0,612.9,668.6) # Sol array\n",
    "    \n",
    "    # polynomial fit\n",
    "    p = np.polyfit(lsdata, soldata, n)\n",
    "    \n",
    "    # sol for Ls input\n",
    "    sol = np.polyval(p, ls)\n",
    "    \n",
    "    # return result\n",
    "    return sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formula to convert solar longitude to Sols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Sol'] = ls2sol(df['L_s'],8)\n",
    "dfday = df[df['LTST'] >= .5]\n",
    "dfnight = df[df['LTST'] <= .5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cond = (df['MY'],df['L_s'],pco2(df['T_surf']), df['p_surf'], tco2(df['p_surf']), df['T_surf'], df['LTST'], df['Surf_lat'], df['Surf_lon'], df['Dust_column'] , df['Sol'] )\n",
    "cond_right= np.transpose(cond)\n",
    "CO2 = pd.DataFrame(list(cond_right))\n",
    "CO2.columns = ['MY','L_s','Eq_Vap_P','p_surf','Frost_T','T_surf','LTST', 'Surf_lat', 'Surf_lon', 'Dust_column', 'Sol']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "newCO2 = CO2[CO2['Eq_Vap_P'] <= CO2['p_surf']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "CO2day = CO2[CO2['LTST'] >= 0.5]\n",
    "CO2night = CO2[CO2['LTST'] <= 0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def MYretrieval(dataframe):\n",
    "    dataframes = []\n",
    "    for i in range(8):\n",
    "        dataframes.append(dataframe[(dataframe['MY'] == i + 28)])\n",
    "    return dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def CreateSubFrame(MarsDF, RanL_s, RanLat, RanLong):\n",
    "    '''Function to divide previously created data frames in smaller data frames by providing ranges\n",
    "    Input:\n",
    "    MarsYearDF: Mars Data frames \n",
    "    RanL_s: range of Solar Longitude example [0,360]\n",
    "    RanAlbedo: Albedo Range example [0.5,1.0]\n",
    "    RanSZA Solar Zenith angle range example [0,45]\n",
    "    RanLat and RanLong are set as default\n",
    "    Output: Several Subframes'''\n",
    "    L_s_MYDF = MarsDF['L_s']\n",
    "    Lat_MYDF = MarsDF['Surf_lat']\n",
    "    Long_MYDF = MarsDF['Surf_lon']\n",
    "    MY_SUB = MarsDF[(L_s_MYDF <= RanL_s[-1])  & (L_s_MYDF >= RanL_s[0]) &\n",
    "                        (Lat_MYDF <= RanLat[-1] ) & (Lat_MYDF >= RanLat[0]) & \n",
    "                        (Long_MYDF <= RanLong[-1]) & (Long_MYDF >= RanLong[0])]\n",
    "    return MY_SUB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def column(MYDict, Parameter):\n",
    "    '''Function to run statistics on data on the subframes from the dictionaries created by the SubDict function\n",
    "    Input:\n",
    "    MYDict: Dictionary Created by the SubDict Function\n",
    "    Output:\n",
    "    SubAvg: Average of the Albedo for each subframe in the dictionary\n",
    "    SubStd: Standard Deviation of the Albedo for each subframe in the dictionary\n",
    "    Sub_L_s: Average of the  Solar Longitude range for each subframe in the dictionary'''\n",
    "    SubframeNum = len(MYDict.keys())\n",
    "    SubTemp = []\n",
    "    SubPres =[]\n",
    "    SubMY = []\n",
    "    SubLat =[]\n",
    "    SubLong = []\n",
    "    SubDust =[]\n",
    "    SubSol = []\n",
    "    SubStr = 'DataFrame{}'\n",
    "    for i in range(SubframeNum):\n",
    "        Subframe = MYDict[SubStr.format(i)]\n",
    "        SubTemp.append((Subframe['T_surf'].values))\n",
    "        SubPres.append((Subframe['p_surf'].values))\n",
    "        SubMY.append((Subframe['MY'].values))\n",
    "        SubLat.append((Subframe['Surf_lat'].values))\n",
    "        SubLong.append((Subframe['Surf_lon'].values))\n",
    "        SubDust.append((Subframe['Dust_column'].values))\n",
    "        SubSol.append((Subframe['Sol'].values))\n",
    "    SubTemp = np.array(SubTemp)\n",
    "    SubPres = np.array(SubPres)\n",
    "    SubMY = np.array(SubMY)\n",
    "    SubLat = np.array(SubLat)\n",
    "    SubLong = np.array(SubLong)\n",
    "    SubDust = np.array(SubDust)\n",
    "    SubSol = np.array(SubSol)\n",
    "    if Parameter =='Temperature':\n",
    "        return SubTemp\n",
    "    elif Parameter =='Pressure': \n",
    "        return SubPres\n",
    "    elif Parameter =='Martian Year':\n",
    "        return SubMY \n",
    "    elif Parameter =='Latitude':\n",
    "        return SubLat \n",
    "    elif Parameter =='Longitude':\n",
    "        return SubLong\n",
    "    elif Parameter =='Dust Column':\n",
    "        return SubDust\n",
    "    elif Parameter =='Sol':\n",
    "        return SubSol "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def IRFlux(T):\n",
    "    '''Calculate the luminosity of the thing.\n",
    "    Input:\n",
    "    T = Temperature\n",
    "    Output:\n",
    "    boltz = luminosity'''\n",
    "    sigma = 5.67e-8 # W/m**2/K**4\n",
    "    boltz = sigma * T**4\n",
    "    return boltz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def emi(temp,frosttemp):\n",
    "    sigma = 5.67e-8\n",
    "    a = IRFlux(temp)\n",
    "    b = IRFlux(frosttemp)\n",
    "    emi = a/b\n",
    "    return emi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def square(lat,lon,size=59.157935, type = 'km'): # 1 degree latitude (Co-created with Tyler Horvath)\n",
    "    '''Function to draw equal length trapezoid that are latitude dependent\n",
    "    Inputs:\n",
    "    Latitude = array or float, latitude of the locations\n",
    "    Longitude = array or float, longitude of the locations\n",
    "    Size = size of the box default is equal to one degree by one degree could be kilometers or degrees\n",
    "    type = string, if input is km the size argument is inputed as kilometers else, deg is inputed as degrees\n",
    "    Outputs:\n",
    "    lower_lat = lower latitude \n",
    "    upper_lat = upper latitude \n",
    "    right_lon = right longitude corner\n",
    "    left_lon = left longitude corner'''\n",
    "    radius = 3389.5 # mars radius in km\n",
    "    d2r = np.pi/180\n",
    "    if type == 'km': # treat size input as km if >= 59....\n",
    "        lat_deg = size*180/(np.pi*radius)\n",
    "        upper_lat = lat + .5 * lat_deg\n",
    "        lower_lat = lat - .5 * lat_deg\n",
    "        lon_deg = size / (np.pi * np.cos(upper_lat*d2r) / 180)\n",
    "        left_lon = lon - .5 * lon_deg * size\n",
    "        right_lon = lon + .5 * lon_deg * size\n",
    "        return(lower_lat, upper_lat, right_lon, left_lon)\n",
    "    elif type == 'deg':\n",
    "        upper_lat = lat + .5 * size\n",
    "        lower_lat = lat - .5 * size\n",
    "        leng = radius * np.pi * np.cos(upper_lat*d2r)*2\n",
    "        lon_deg = 59.157935*size/leng*360\n",
    "        left_lon = lon - .5 * lon_deg \n",
    "        right_lon = lon + .5 * lon_deg\n",
    "        return(lower_lat, upper_lat, right_lon, left_lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def xerr(a,b,covmat):\n",
    "    '''Function to calculate the x-intercept of a fitted line in the form y = ax+b with its error\n",
    "    Inputs:\n",
    "    a = float, the slope of the fitted line\n",
    "    b = float, the x intercept of the line\n",
    "    covmat = array, the covariance matrix of the fitted line\n",
    "    Outputs:\n",
    "    x = float, the x-intercept of the line\n",
    "    stdx = float , error of the line '''\n",
    "    x = -b/a\n",
    "    stdx = x*np.sqrt((covmat[1,1]/(b**2))+(covmat[0,0]/(a**2))-2*((covmat[0,1]/(a*b))))\n",
    "    return x, stdx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 361,
   "metadata": {
    "code_folding": []
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-361-f04cbf79e9e9>, line 15)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-361-f04cbf79e9e9>\"\u001b[0;36m, line \u001b[0;32m15\u001b[0m\n\u001b[0;31m    for i in range(len())\u001b[0m\n\u001b[0m                         ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "def LocFunc(lat,lon,size=59.157935,type = 'km'):\n",
    "    Location = square(lat, lon, size=2, type = 'deg')\n",
    "    Day = MYretrieval(dfday)\n",
    "    Night = MYretrieval(dfnight)\n",
    "    CO2D = MYretrieval(CO2day)\n",
    "    CO2N = MYretrieval(CO2night)\n",
    "    for i in range(len(Day)):\n",
    "        Day[i].insert(10 , 'Frost_T', CO2D[i]['Frost_T'])\n",
    "        Night[i].insert(10 , 'Frost_T', CO2D[i]['Frost_T'])\n",
    "    LocationDay = []\n",
    "    LocationNight =[]\n",
    "    for i in range(len(Day)):\n",
    "        LocationDay.append(CreateSubFrame(Day[i], RanL_s=np.linspace(0,360), RanLat=np.linspace(Location[0],Location[1]), RanLong=np.linspace(Location[3],Location[2])))\n",
    "        LocationNight.append(CreateSubFrame(Night[i], RanL_s=np.linspace(0,360), RanLat=np.linspace(Location[0],Location[1]), RanLong=np.linspace(Location[3],Location[2])))\n",
    "    for i in range(len())\n",
    "    MaxTempDay = []\n",
    "    for i in range(len(LocationDay)):\n",
    "        MaxTempDay.append(np.max(LocationDay[i]['T_surf'])[i])\n",
    "    #FirstFrostDay = []\n",
    "    #for i in range(len(LocationDay)):\n",
    "    #    if lat >= 58.0: \n",
    "    #        FirstFrostDay.append(np.min(LocationDay[i]['L_s'][(LocationDay[i]['L_s'] >= 180.0) & (LocationDay[i]['T_surf'] <= LocationDay[i]['Frost_T'])]))\n",
    "    #    elif (lat < 58.0 and lat > -58.0):\n",
    "    #        FirstFrostDay.append(np.nan)\n",
    "    #    elif lat <= -58.0: \n",
    "    #        FirstFrostDay.append(np.min(LocationDay[i]['L_s'][(LocationDay[i]['L_s'] <= 180.0) & (LocationDay[i]['T_surf'] <= LocationDay[i]['Frost_T'])]))\n",
    "    #LastFrostDay = []\n",
    "    #for i in range(len(LocationDay)):\n",
    "    #    if lat >= 58.0: \n",
    "    #        FirstFrostDay.append(np.max(LocationDay[i]['L_s'][(LocationDay[i]['L_s'] >= 180.0) & (LocationDay[i]['T_surf'] <= LocationDay[i]['Frost_T'])]))\n",
    "    #    elif (lat < 58.0 and lat > -58.0):\n",
    "    #        FirstFrostDay.append(np.nan)\n",
    "    #    elif lat <= -58.0: \n",
    "    #        FirstFrostDay.append(np.max(LocationDay[i]['L_s'][(LocationDay[i]['L_s'] <= 180.0) & (LocationDay[i]['T_surf'] <= LocationDay[i]['Frost_T'])]))\n",
    "    ##FirstFrostNight = []\n",
    "    ##FirstFrostNight = []\n",
    "    #    #if lat >= 0\n",
    "    #    #    FirstFrostNight.append(np.min(LocationNight[i]['L_s'][(LocationNight[i]['L_s'] >= 180.0) & (LocationNight[i]['T_surf'] <= LocationNight[i]['Frost_T'])]))\n",
    "    #    #elif lat <= -58.0: \n",
    "    #    #    FirstFrostNight.append(np.min(LocationNight[i]['L_s'][(LocationNight[i]['L_s'] <= 180.0) & (LocationNight[i]['T_surf'] <= LocationNight[i]['Frost_T'])]))\n",
    "    return LocationDay, LocationNight, MaxTempDay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 359,
   "metadata": {},
   "outputs": [],
   "source": [
    "BuzzelDay, BuzzelNight, FirstFrostDay = LocFunc(-58.0, -53.2, size=2, type = 'deg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 360,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[243.53200000000001,\n",
       " 273.31599999999997,\n",
       " 273.435,\n",
       " 275.09100000000001,\n",
       " 275.35700000000003,\n",
       " 277.88099999999997,\n",
       " 160.893,\n",
       " nan]"
      ]
     },
     "execution_count": 360,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FirstFrostDay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = BuzzelDay[1]['T_surf']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pandas.core.series.Series"
      ]
     },
     "execution_count": 345,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 309,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[nan, nan, nan, nan, nan, nan, nan, nan]"
      ]
     },
     "execution_count": 309,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FirstFrostNight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.min(BuzzelDay[1]['L_s'][(BuzzelDay[1]['L_s'] <= 180.0) & (BuzzelDay[1]['T_surf'] <= BuzzelCO2D[1]['Frost_T'])])\n",
    "b = np.min(BuzzelDay[1]['L_s'][(BuzzelDay[1]['L_s'] >= 180.0) & (BuzzelDay[1]['T_surf'] <= BuzzelCO2D[1]['Frost_T'])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2654700000000001"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "187.99021999999999"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}