{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import folium\n",
    "import geopandas as gpd\n",
    "import rasterio\n",
    "from rasterio import plot\n",
    "import numpy as np\n",
    "from rasterio.mask import mask\n",
    "import os\n",
    "from datetime import date, datetime, timedelta\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## NIFC has replaced GeoMAC as a the official source of perimeters in the US.\n",
    "## There are 2 files, active and archived. The active file below has the most recent perimeters for the current year.\n",
    "fires = gpd.read_file(r'https://opendata.arcgis.com/datasets/5da472c6d27b4b67970acc7b5044c862_0.geojson')\n",
    "# Alternative when NIFC down https://services9.arcgis.com/RHVPKKiFTONKtxq3/ArcGIS/rest/services/USA_Wildfires_v1/FeatureServer/1 (manual download)\n",
    "\n",
    "# Import county data and filter by state\n",
    "counties = gpd.read_file(r'D:\\data\\us_counties.geojson')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fire and state where it is located\n",
    "firename = 'Glass'\n",
    "state = 'CA'\n",
    "fireyear = '2020'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Update the basepath based on the most recent data update maintained in the list below\n",
    "current = '2019'\n",
    "prev = '2018'\n",
    "\n",
    "# Where the score tifs are located\n",
    "baset = r'D:/Geotiff/'\n",
    "\n",
    "# Use string, not os.path.join, for gdal command later in the script\n",
    "if os.path.exists(baset + current):\n",
    "    basetif = baset + current\n",
    "else:\n",
    "    basetif = baset + prev\n",
    "\n",
    "# Where to store the outputs\n",
    "base = r'D:\\\\data'\n",
    "\n",
    "# Score tif\n",
    "score1file = os.path.join(basetif, state, state + '_score1.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def invalidGeom(infile):\n",
    "    #if not infile.is_valid.bool():\n",
    "    f = infile['geometry'].buffer(0)\n",
    "    infile['geometry'] = f\n",
    "    return infile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Format the fire names in the NIFC file, clip to get the county/state\n",
    "county = counties.loc[counties.state == state]\n",
    "fires.IncidentName=fires.IncidentName.str.title().str.strip()\n",
    "\n",
    "invalidGeom(fires)\n",
    "\n",
    "firestate = gpd.clip(fires,county)\n",
    "firestate['state'],firestate['county'] = state,county.fl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the fire by name\n",
    "fire = firestate.loc[firestate.IncidentName == firename]\n",
    "                     \n",
    "# Check that the correct fire is selected\n",
    "if len(fire) < 1:\n",
    "    print('Search for a manually created perimeter file')\n",
    "    try:\n",
    "        fire = gpd.read_file(os.path.join(r'Z:\\FirePerimeters', fireyear, state,firename, firename.replace(' ','_')+'.gpkg'))\n",
    "        print('Selected a manually created file')\n",
    "    except:\n",
    "        print('Try complex name')\n",
    "    if len(fire) < 1:\n",
    "        try:\n",
    "        # Try for a manually created perimeter file (Remember to create 'IncidentName' and 'DateCurrent' fields)\n",
    "            fire = fires.loc[(~fires.ComplexName.isnull()) & (fires.ComplexName.str.contains(firename))]\n",
    "            fire.IncidentName = fire.ComplexName\n",
    "            fire['GISAcres'] = fire['GISAcres'].sum()\n",
    "            fire = fire.dissolve(by='ComplexName')\n",
    "            print('Selected by Complex Name')\n",
    "        except:\n",
    "            print('Try partial string')\n",
    "        if len(fire) < 1:\n",
    "            try:\n",
    "                fire = firestate.loc[firestate.IncidentName.str.contains(firename.title().strip())]\n",
    "            except:\n",
    "                print('No matches')\n",
    "\n",
    "elif len(fire) > 1:\n",
    "    try:\n",
    "        fire = firestate.loc[firestate.IncidentName.str.contains(firename.title().strip())]\n",
    "        print('Multiple fires selected')\n",
    "    except:\n",
    "        print('No matches')\n",
    "\n",
    "    \n",
    "\n",
    "# Sometimes the NIFC data has invalid geometries, use 0 width buffer to fix\n",
    "if len(fire) == 1:\n",
    "    invalidGeom(fire)\n",
    "\n",
    "print(str(len(fire)) + ' fire(s) selected')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If multiple fires, select the correct fire\n",
    "#fire = fire.iloc[[1]] #Double bracket to maintain datatype when using iloc with geopandas\n",
    "\n",
    "fire.iloc[0].DateCurrent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "def boundsfire(df):\n",
    "    # Get the bounds of the fire for the folium map\n",
    "    minmax = ['minx','miny','maxx','maxy']\n",
    "    for f in minmax:\n",
    "        df[f] = df.bounds[f]\n",
    "    return df[f] \n",
    "\n",
    "# Create the date field for the marketing text\n",
    "def dateNIFC(datefield):\n",
    "    datefield = pd.to_datetime(datefield)\n",
    "    dt_string = datefield.strftime(\"%Y-%m-%d %H:%M:%S\")\n",
    "    dt = datefield.strftime(\"%m-%d-%Y\")\n",
    "    return dt\n",
    "\n",
    "date = dateNIFC(fire.iloc[0].DateCurrent)\n",
    "\n",
    "boundsfire(fire)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a gpkg for the fire to clip the score tif data\n",
    "# Remove any spaces from fire name for output files\n",
    "x = str(fire.iloc[0].IncidentName).replace(' ','_').replace('-','_')\n",
    "\n",
    "# fire data output\n",
    "basepath = os.path.join(base,fireyear,state,fire.iloc[0].IncidentName)\n",
    "if not os.path.exists(basepath):\n",
    "    os.makedirs(basepath)\n",
    "\n",
    "fireshp = os.path.join(basepath,x+'.gpkg')\n",
    "fireshpshp = os.path.join(basepath,x+'.shp')\n",
    "\n",
    "fire.to_file(fireshp, driver=\"GPKG\",OVERWRITE='YES')\n",
    "fire.to_file(fireshpshp,OVERWRITE='YES')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# paths for gdal command\n",
    "outfirepath = basepath +'\\\\'+ x +'.tif'\n",
    "infireshppath = basepath +'\\\\'+ x +'.gpkg'\n",
    "\n",
    "outfire = str('\"' +outfirepath+'\"')\n",
    "inscore1 = str('\"' +score1file+'\"')\n",
    "infireshp = str('\"' +infireshppath+'\"')\n",
    "\n",
    "# Run gdal cmd to cut the raster to the fire polygon\n",
    "# Ensure connection through vpn\n",
    "cmd = \"gdalwarp -cutline \"+ infireshp +\" -crop_to_cutline -dstalpha -dstnodata 255  -overwrite \" + inscore1 + \" \" +  outfire\n",
    "os.system(cmd)\n",
    "print(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reclassify(src):\n",
    "    # Reclassfy the score for rendering in the map\n",
    "    global img\n",
    "    img = src.read(1)\n",
    "    img[img==0] = 0 #Negligible\n",
    "    img[img==1] = 1 #Low\n",
    "    img[np.where((img >=2) & (img <=3))] = 2 #Moderate\n",
    "    img[np.where((img>=4) & (img<=12))] = 3 #High\n",
    "    img[np.where((img>=13) & (img<=30))] = 4 #Extreme\n",
    "    return img\n",
    "\n",
    "def pixelCount(r, vals):\n",
    "    # Create a list of the raster values to summarize % of score categories\n",
    "    val_list = []\n",
    "    size = float(r.size)\n",
    "    for v in vals:\n",
    "        count = np.sum(r == v)\n",
    "        val_list.append(count)#/size for pct - can't use with alpha band/nodata\n",
    "    return val_list\n",
    "\n",
    "def pixelCountTest(r, vals):\n",
    "    # Create a list of the raster values to summarize % of score categories\n",
    "    val_list = []\n",
    "    size = float(r.size)\n",
    "    for v in vals:\n",
    "        count = np.sum(r == v)\n",
    "        pct = float(count/size)\n",
    "        val_list.append(pct)#/size for pct - can't use with alpha band/nodata\n",
    "    return val_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basemaps\n",
    "def cartoBasemap(miny, maxx):\n",
    "    global m\n",
    "    m = folium.Map([fire.miny, fire.maxx], zoom_start=5, tiles='cartodbpositron')\n",
    "    return m\n",
    "\n",
    "def esriBasemap(miny, maxx):\n",
    "    global m\n",
    "    tile = \"https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}\"\n",
    "    attribu = 'Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community'\n",
    "    m = folium.Map([fire.miny, fire.maxx], zoom_start=9, tiles=tile,attr=attribu)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "esriBasemap(fire.iloc[0].miny, fire.iloc[0].maxx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# score1 display colors \n",
    "# r,g,b,a - change a to 0 for translucent, 1 for opaque\n",
    "score1_color = {\n",
    " 0: (233,233,238,1),\n",
    " 1: (243,242,6,1),\n",
    " 2: (253,187,15,1),\n",
    " 3: (244,119,33,1),\n",
    " 4: (208,22,38,1),\n",
    " 255: (0, 0, 0, 0)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add fire raster to map and calculate the % in each score category\n",
    "# Create folium map\n",
    "with rasterio.open(outfirepath) as src:\n",
    "    reclassify(src)\n",
    "    \n",
    "    ### Temporary workaround - lambda function assigns colormap incorrectly if the raster does not start at 0\n",
    "    ### Add 1 '0' value cell in west corner\n",
    "    img[0,0] = 0\n",
    "        \n",
    "    ## Folium map - optional for viewing fire\n",
    "    folium.raster_layers.ImageOverlay(img, colormap=(lambda x: score1_color[x]),\n",
    "                                 bounds =[[src.bounds.top, src.bounds.left], [src.bounds.bottom, src.bounds.right]]).add_to(m)\n",
    "    # calculate the pixel count for each score category\n",
    "    results = []\n",
    "    val_list = [0,1, 2, 3, 4] # score reclass values \n",
    "    results.append(pixelCount(img, val_list))\n",
    "    print(results)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Folium map - optional\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the % of WFHC for the fire\n",
    "def roundCalc(rs,x,total):\n",
    "    if rs[0][x] == 0:\n",
    "        y = 0\n",
    "    else:\n",
    "        y = (rs[0][x]/total)*100\n",
    "    x = round(y)\n",
    "    return x\n",
    "\n",
    "# Round to 1 decimal place\n",
    "def rCalc(rs,x,total):\n",
    "    if rs[0][x] == 0:\n",
    "        y = 0\n",
    "    else:\n",
    "        y = (rs[0][x]/total)*100\n",
    "    x = round(y,1)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Stats for blurb\n",
    "\n",
    "# Tif file names\n",
    "append = ['_score1.tif','_score2.tif','_score3.tif','_score4.tif']\n",
    "firename_string = x+\"_\"+str(date)\n",
    "\n",
    "# Create fire tifs\n",
    "for v in append:\n",
    "    tif, outtif = os.path.join(basetif,state,state+v), os.path.join(basepath,firename_string+v)\n",
    "    intif, outfire = str('\"' +tif+'\"'),str('\"' +outtif+'\"')\n",
    "    cmd = \"gdalwarp -cutline \"+ infireshp +\" -crop_to_cutline -dstnodata 255 -overwrite \" + intif + \" \" +  outfire \n",
    "    # remove -dstalpha alpha band for symbology in arcmap, add for visualizing with folium\n",
    "    os.system(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Maintain order of the values/variables in the lists below\n",
    "## Be sure to check for significant non-burnable/negligible areas - remove from analysis\n",
    "\n",
    "# calculate the pixel counts for the categorical rasters \n",
    "# Order of the score below correspond to raster values 0,1,2,3,4,5,6. 'null' for values not used in each score category\n",
    "score1 = ['score11','score12','score13','score14','score15','null']\n",
    "score2 = ['score21','score22','null','score23','score24','score25','score26']\n",
    "score3 = ['null','score31','score32','score33','null','score34','null']\n",
    "score4 = ['score41','score42','null','score43','null','score44','null']\n",
    "\n",
    "variables = [score1,score2,score3,score4]\n",
    "variable_list = ['score1','score2','score3','score4'] # Values are placeholders for calculated stats\n",
    "\n",
    "# Open each score raster and calculate the score distribution\n",
    "n = 0\n",
    "print_list = pd.DataFrame(columns=['Score','Value','pct'])\n",
    "for v in append:\n",
    "    varname_list = variables[n]\n",
    "    outtif = os.path.join(basepath,firename_string+v)\n",
    "    if os.path.exists(outtif):\n",
    "        with rasterio.open(outtif) as src:\n",
    "            if varname_list == score1:# Reclassify score1 only\n",
    "                reclassify(src)\n",
    "            else:\n",
    "                img = src.read(1)\n",
    "\n",
    "            # calculate the pixel count for each score category\n",
    "            r = []\n",
    "            val_list = [0, 1, 2, 3, 4, 5, 6] \n",
    "            r.append(pixelCount(img, val_list))\n",
    "\n",
    "            # Get total number of pixels to divide by\n",
    "            t = r[0][0]+r[0][1]+r[0][2]+r[0][3]+r[0][4]+r[0][5]+r[0][6]\n",
    "\n",
    "            # Calculate the % of each value and reassign each score category to variable_list\n",
    "            v0,v1,v2,v3,v4,v5,v6 = rCalc(r,0,t),rCalc(r,1,t),rCalc(r,2,t),rCalc(r,3,t),rCalc(r,4,t),rCalc(r,5,t),rCalc(r,6,t)\n",
    "            vname = variable_list[n] \n",
    "            variable_list[n] = [v0,v1,v2,v3,v4,v5,v6]   \n",
    "\n",
    "            # Print the % in each score category formatted for marketing\n",
    "            y = 0\n",
    "            for var in varname_list:\n",
    "                if varname_list[y] == 'null':\n",
    "                    a=1 # skip 'null'\n",
    "                else:\n",
    "                    print(str(varname_list[y])+\": \"+str(variable_list[n][y])+\"%\")\n",
    "                    print_list = print_list.append({'Score': vname,'Value': (str(varname_list[y])),'pct': (str(variable_list[n][y]))}, ignore_index=True)\n",
    "\n",
    "                # Get high & extreme pct\n",
    "                if varname_list[y] == 'score14':\n",
    "                    sc1 = variable_list[n][y]\n",
    "                if varname_list[y] == 'score15':\n",
    "                    sc2 = variable_list[n][y]\n",
    "                y = y + 1\n",
    "    else:\n",
    "        print(state + v + ' does not exist')\n",
    "            \n",
    "    n = n + 1\n",
    "    print('')\n",
    "\n",
    "# Combine High & Extreme for marketing snippet\n",
    "sc1sc2 = round(sc1 + sc2,2)\n",
    "\n",
    "acres = 'As of ' + date + ', the '+ str(firename)+ ' Fire has burned approximately ' + str(round(fire.iloc[0].GISAcres,0))+ ' acres in ' + state +' according to NIFC perimeter data.'\n",
    "risk = 'X has determined that approximately ' +str(sc1sc2)+'% of the area is in X risk categories.'\n",
    "print(marketing,risk)\n",
    "\n",
    "print_list = print_list.append({'Score': acres,'Value': risk}, ignore_index=True)\n",
    "print_list.to_csv(os.path.join(basepath, firename_string+'_scores.csv'), index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reset fire name for command line arguments in MXD creation script\n",
    "print(firename)\n",
    "if firename != fire.iloc[0].IncidentName:\n",
    "    firename = fire.iloc[0].IncidentName"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the MXD and pdf map with arcpy scripts\n",
    "score_1 = \"c:/pythonenv/arcgis10.6/python D:\\scripts\\score1map.py \" + fireyear +\" \"+ state +\" \"+ '\"'+firename+ '\"'\n",
    "score_2 = \"c:/pythonenv/arcgis10.6/python D:\\scripts\\score2map.py \" + fireyear +\" \"+ state +\" \"+ '\"'+firename+ '\"'\n",
    "\n",
    "os.system(score_1)\n",
    "os.system(score_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}