{ "cells": [ { "cell_type": "code", "execution_count": 3, "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": 4, "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", "# Lots of issues with invalid geometries - remove them\n", "#fires = fires.loc[fires.is_valid]\n", "\n", "# Import county data and filter by state\n", "counties = gpd.read_file(r'Z:\\FirePerimeters\\fl_counties.geojson')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Fire and state where it is located\n", "firename = 'Glass'\n", "state = 'CA'\n", "fireyear = '2020'" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Update the basepath based on the most recent FireLine update maintained in the list below\n", "current = '2019'\n", "prev = '2018'\n", "\n", "# Where the score tifs are located\n", "baset = r'Z:\\ArcFireLine\\FireRing\\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'Z:\\\\FirePerimeters'\n", "\n", "# Score tif\n", "wfhsfile = os.path.join(basetif, state, state + '_wfhs.tif')" ] }, { "cell_type": "code", "execution_count": 7, "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": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Python38\\lib\\site-packages\\geopandas\\geoseries.py:358: UserWarning: GeoSeries.notna() previously returned False for both missing (None) and empty geometries. Now, it only returns False for missing values. Since the calling GeoSeries contains empty geometries, the result has changed compared to previous versions of GeoPandas.\n", "Given a GeoSeries 's', you can use '~s.is_empty & s.notna()' to get back the old behaviour.\n", "\n", "To further ignore this warning, you can do: \n", "import warnings; warnings.filterwarnings('ignore', 'GeoSeries.notna', UserWarning)\n", " return self.notna()\n" ] } ], "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": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 fire(s) selected\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":4: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame.\n", "Try using .loc[row_indexer,col_indexer] = value instead\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", " infile['geometry'] = f\n" ] } ], "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 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": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'2020-09-30T03:59:52'" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "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": 11, "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": 12, "metadata": {}, "outputs": [], "source": [ "# Create a gpkg for the fire to clip the score tif data\n", "# Remove any spaces from fire name\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": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "gdalwarp -cutline \"Z:\\\\FirePerimeters\\2020\\CA\\Glass\\Glass.gpkg\" -crop_to_cutline -dstalpha -dstnodata 255 -overwrite \"Z:\\ArcFireLine\\FireRing\\Geotiff/2019\\CA\\CA_wfhs.tif\" \"Z:\\\\FirePerimeters\\2020\\CA\\Glass\\Glass.tif\"\n" ] } ], "source": [ "# paths for gdal command\n", "outfirepath = basepath +'\\\\'+ x +'.tif'\n", "infireshppath = basepath +'\\\\'+ x +'.gpkg'\n", "\n", "outfire = str('\"' +outfirepath+'\"')\n", "inwfhs = str('\"' +wfhsfile+'\"')\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 \" + inwfhs + \" \" + outfire\n", "os.system(cmd)\n", "print(cmd)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "def reclassify(src):\n", " # Reclassfy the WFHS 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": 15, "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 © Esri — 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": 16, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "esriBasemap(fire.iloc[0].miny, fire.iloc[0].maxx)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# WFHS display colors \n", "# r,g,b,a - change a to 0 for translucent, 1 for opaque\n", "wfhs_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": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1424, 3347, 24891, 183990, 30738]]\n" ] } ], "source": [ "# Add fire raster to map and calculate the % in each score category\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: wfhs_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] # WFHS reclass values \n", " results.append(pixelCount(img, val_list))\n", " print(results)\n", " " ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## Folium map - optional\n", "m" ] }, { "cell_type": "code", "execution_count": 20, "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": 21, "metadata": {}, "outputs": [], "source": [ "## Stats for marketing\n", "\n", "# Tif file names\n", "append = ['_wfhs.tif','_fuel.tif','_slope.tif','_access.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": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Negligible: 0.6%\n", "Low: 1.4%\n", "Moderate: 10.2%\n", "High: 75.3%\n", "Extreme: 12.6%\n", "\n", "F0: 0%\n", "F1: 40.2%\n", "F3: 32.8%\n", "NF: 2.3%\n", "F5: 24.5%\n", "WA: 0.3%\n", "\n", "S1: 3.3%\n", "S2: 28.1%\n", "S3: 49.9%\n", "S5: 18.8%\n", "\n", "A0: 97.6%\n", "A1: 1.6%\n", "A3: 0.4%\n", "A5: 0.4%\n", "\n", "As of 09-30-2020, the Glass Fire has burned approximately 50962.0 acres in CA according to NIFC perimeter data. FireLine has determined that approximately 87.9% of the area is in High or Extreme risk categories.\n" ] } ], "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", "# Order of the score below correspond to raster values 0,1,2,3,4,5,6. 'null' if score/value does not exist\n", "wfhs = ['Negligible','Low','Moderate','High','Extreme','null']\n", "fuel = ['F0','F1','null','F3','NF','F5','WA']\n", "slope = ['null','S1','S2','S3','null','S5','null']\n", "access = ['A0','A1','null','A3','null','A5','null']\n", "\n", "variables = [wfhs,fuel,slope,access]\n", "variable_list = ['wfhs','fuel','slope','access'] # 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 == wfhs:# Reclassify WFHS 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] == 'High':\n", " High = variable_list[n][y]\n", " if varname_list[y] == 'Extreme':\n", " Ext = 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", "highext = round(High + Ext,2)\n", "\n", "marketing = '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 = 'FireLine has determined that approximately ' +str(highext)+'% of the area is in High or Extreme risk categories.'\n", "print(marketing,risk)\n", "\n", "print_list = print_list.append({'Score': marketing,'Value': risk}, ignore_index=True)\n", "print_list.to_csv(os.path.join(basepath, firename_string+'_scores.csv'), index=False)\n", "\n", "# If Negligible risk category accounts for a large portion of the perimeter, check for/remove non-burnable areas" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Glass\n" ] } ], "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\n" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create the MXD and pdf map\n", "fuel = \"c:/python27/arcgis10.5/python Z:\\FirePerimeters\\CreateMXD_Fuels.py \" + fireyear +\" \"+ state +\" \"+ '\"'+firename+ '\"'\n", "wfhs = \"c:/python27/arcgis10.5/python Z:\\FirePerimeters\\CreateMXD_WFHS.py \" + fireyear +\" \"+ state +\" \"+ '\"'+firename+ '\"'\n", "\n", "os.system(wfhs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.3" } }, "nbformat": 4, "nbformat_minor": 4 }