# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="W2zYi5O_Mcgy"
# # Extract daily weather variables based on unique ID 
#

# + [markdown] id="fZCPD9SxMzMe"
# ## Setup

# + [markdown] id="Z2A9PoqQUye1"
# ### Connect to Drive

# + id="eMCHZj6PpL3w" colab={"base_uri": "https://localhost:8080/", "height": 71} executionInfo={"status": "ok", "timestamp": 1601659777602, "user_tz": 420, "elapsed": 585, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="827ea2cd-e5ce-4307-9d3c-8d256c4501b0"
from google.colab import drive # import drive from google colab
ROOT = "/content/drive"     # default location for the drive
print(ROOT)                 # print content of ROOT (Optional)

drive.mount(ROOT)           # we mount the google drive at /content/drive

# + id="y237qvlb5kwV"
import os
rootPath = "/content/drive/My Drive/California FireTrends (2012-2020)"
os.chdir(rootPath)

# + [markdown] id="Ag6HtRJpNNWC"
# ### Connect to GEE API

# + id="9P_WfKSZNbJ4"
# initialize and connect to GEE 
from google.colab import auth
auth.authenticate_user()
# !earthengine authenticate
import ee 
ee.Initialize()

# + id="OaUQbtDpNbPB" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601569780648, "user_tz": 420, "elapsed": 30338, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="f8bf74c0-f2b0-4e00-af5c-c25edbf84a8d"
# Installs geemap package
import subprocess

try:
    import geemap
except ImportError:
    print('geemap package not installed. Installing ...')
    subprocess.check_call(["python", '-m', 'pip', 'install', 'geemap'])

# Checks whether this notebook is running on Google Colab
try:
    import google.colab
    import geemap.eefolium as emap
except:
    import geemap as emap

# + [markdown] id="NBuFTbvcU3Bk"
# ### Connect to Cloud Bucket

# + id="Omn5x960NbTD"
# Connect to google cloud 
# ! gcloud auth login

# + [markdown] id="6dZbRPY8U-9L"
# ### Load Libraries 

# + id="o66RFau97lMx"
# %%time 

# Important library for many geopython libraries
# !apt install gdal-bin python-gdal python3-gdal 
# Install rtree - Geopandas requirment
# !apt install python3-rtree 
# Install Geopandas
# !pip install git+git://github.com/geopandas/geopandas.git
# Install descartes - Geopandas requirment
# !pip install descartes 
# Install Folium for Geographic data visualization
# !pip install folium
# Install plotlyExpress
# !pip install plotly_express

# + id="Q1DSJekd-tgI"
import pandas as pd
import numpy as np
import geopandas as gpd
from shapely.geometry import Point
import matplotlib
import matplotlib.pyplot as plt 
import folium
import rtree


# + id="EGByx1pSVHPV"
def listFiles_ByExt(rootPath, ext):
  '''
  retrieve file path + names based on extension
  '''
  file_list = []
  root = rootPath
  for path, subdirs, files in os.walk(root):
      for names in files: 
          if names.endswith(ext):
              file_list.append(os.path.join(path,names))
  return file_list

def createFolder(rootPath, folderName): 
  '''
  Create new folder in root path 
  '''
  folderPath = os.path.join(rootPath, folderName) 
  if not os.path.exists(folderPath):
      os.makedirs(folderPath)
  return folderPath 


def listFiles_ALL(rootPath):
  '''
  retrieve file path + names 
  '''
  file_list = []
  root = rootPath
  for path, subdirs, files in os.walk(root):
    for names in files: 
      file_list.append(os.path.join(path,names))
  return file_list


# + [markdown] id="H9iaY0vqQuFT"
# ## Prepare Fire Perimeters

# + [markdown] id="MC3zkmMvM0i1"
# ### Combine all interpolated surfaces into one file

# + id="iaUfwQqPXCPc"
rootPath = os.path.join(r'Products', 'By_Fire')
inter_files = listFiles_ByExt(rootPath,'.geojson')

# + id="wvkoi7PCM6oj"
outpath = os.path.join(r'Products', 'ALL_2012_2020')

gdf = pd.concat([
    gpd.read_file(shp)
    for shp in inter_files
    ]).pipe(gpd.GeoDataFrame)

# gdf.to_file(os.path.join(outpath, f'ALL_2012_2020.shp'))
# gdf.to_file(os.path.join(outpath, f'ALL_2012_2020.geojson'), driver='GeoJSON')

# + [markdown] id="J_Ng0tgSM6-Q"
# ### Clip to CA Extent
#
# Currently the study area is a bit outside state lines 

# + id="IqaNTpV4aC0J"
gdf.crs

# + id="Bph5UHfbM-S_"
# check; set CRS 
state_path = '/content/drive/My Drive/COVID_FireTrends/Data/US_states/'
states = gpd.read_file(state_path)
CA_states = states[(states.STATE_NAME == 'California')]
CA_states = CA_states.to_crs('EPSG:3310')
CA_states_buff = CA_states
CA_states_buff['geometry'] = CA_states.geometry.buffer(5000)

# + id="i2N424YyklAs"
CA_states_buff.to_file('Data/CA_Extent/California_5kmbuff.shp')

# + id="o0hIfFZyaNhU"
## clip to california extent 
Interpolated_CA = gpd.overlay(gdf, CA_states_buff, how='intersection')

# + id="J6hS2SFRagaJ" colab={"base_uri": "https://localhost:8080/", "height": 71} executionInfo={"status": "ok", "timestamp": 1601524096613, "user_tz": 420, "elapsed": 28816, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="dd33f912-5c3d-498f-8ca3-ef8179f721b0"
# Interpolated_CA.to_file(os.path.join(outpath, f'ALL_2012_2020.shp'))
# Interpolated_CA.to_file(os.path.join(outpath, f'ALL_2012_2020.geojson'), driver='GeoJSON')

# + [markdown] id="-k3b8LyTM-r4"
# ### Create unique ID column (numeric)
#
# {FIREID}{YR}{JD}

# + id="SEzTFWtzNCV4"
Interpolated_CA.sort_values(by=['Year', 'Fire'])
# Interpolated_CA['FireID'] = Interpolated_CA.groupby(['Fire', 'Year']).cumcount()
Interpolated_CA['FireID'] = 1
mask = Interpolated_CA.groupby(['Fire','Year'])['FireID'].transform(lambda x : len(x)>0)
Interpolated_CA.loc[mask,'FireID'] = Interpolated_CA.loc[mask,['Fire','Year']].astype(str).sum(1).factorize()[0] + 1

# + id="ZNNBv5o0eYag"
Interpolated_CA['Unique_ID'] =  Interpolated_CA['YRJD'].astype('str').str[2:] + Interpolated_CA['FireID'].astype('str').str.zfill(3)
Interpolated_CA['Unique_ID'] = Interpolated_CA['Unique_ID'].astype('int64')

# + id="-U4R3_zji76Y" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601589438701, "user_tz": 420, "elapsed": 705, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="b008f1db-4287-4e0f-d689-9cecda669504"
Interpolated_CA['FireID'].nunique()
Interpolated_CA.groupby(['Fire', 'Year']).ngroups

# + id="2l6M6Tz1efoO" colab={"base_uri": "https://localhost:8080/", "height": 71} executionInfo={"status": "ok", "timestamp": 1601589508236, "user_tz": 420, "elapsed": 28665, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="40e6fa7f-7039-4140-e0ce-e4e6e09596fb"
Interpolated_CA.to_file(os.path.join(outpath, f'ALL_2012_2020.shp'))
Interpolated_CA.to_file(os.path.join(outpath, f'ALL_2012_2020.geojson'), driver='GeoJSON')

# + [markdown] id="Cd0W0QenQmQJ"
# ### Export to Bucket and to Asset

# + id="0VV3oGqVlMwW" colab={"base_uri": "https://localhost:8080/", "height": 156} executionInfo={"status": "ok", "timestamp": 1601589544144, "user_tz": 420, "elapsed": 3502, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AO<KEY>", "userId": "04192059296941555235"}} outputId="4c071bce-2369-40e2-f120-e6a455f632ee"
# ! gsutil -m cp -R 'Data/CA_Extent' gs://daily_fire_surfaces

# + id="kekwfc8OlW84"
# ! earthengine upload table --asset_id=users/escaduto/CA_FireTrends/CA_Extent gs://daily_fire_surfaces/CA_Extent/California_5kmbuff.shp

# + id="bL9LUT36B4J4"
# ! gsutil -m cp -R 'Products/ALL_2012_2020' gs://daily_fire_surfaces

# + id="e_flNFxQCI0k"
# ! earthengine upload table --asset_id=users/escaduto/CA_FireTrends/DerivedPerimeters_2012_2020 gs://daily_fire_surfaces/ALL_2012_2020/ALL_2012_2020.shp

# + id="S3HaWn9hU_od"
# import feature collection asset 
fire_perimeters = ee.FeatureCollection('users/escaduto/CA_FireTrends/DerivedPerimeters_2012_2020')
CA_Extent = ee.FeatureCollection('users/escaduto/CA_FireTrends/CA_Extent')

# + id="LH2rEmbFh2PT"
outpath = os.path.join(r'Products', 'ALL_2012_2020')
Interpolated_CA = gpd.read_file(os.path.join(outpath, f'ALL_2012_2020.geojson'))

# + id="EgsxL2XqiBMu" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601589599529, "user_tz": 420, "elapsed": 531, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="5704c400-609c-4937-f30d-067c32d29f26"
Interpolated_CA['FireID'].nunique()


# + [markdown] id="xGilwZkINGAz"
# ## Functions to GET Layers from GEE

# + [markdown] id="75YkHhPGvywJ"
# #### Cloud Masks

# + id="CjM6oY1xAy6y"
def maskL8sr(image):  
    # Bits 3 and 5 are cloud shadow and cloud, respectively.
    cloudShadowBitMask = (1 << 3);
    cloudsBitMask = (1 << 5);
    # Get the pixel QA band.
    qa = image.select('pixel_qa');
    # Both flags should be set to zero, indicating clear conditions.
    mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) and (qa.bitwiseAnd(cloudsBitMask).eq(0));
    return image.updateMask(mask);


# + [markdown] id="nvp7gSrPvtMw"
# Landsat 7 Mask

# + id="mPpMSfrzZaUI"
# cloud m
def cloudMaskL457(image):
  cloudShadowBitMask = (1 << 3);
  cloudsBitMask = (1 << 5);
  cloudsconfidenceBitMask = (1 << 7);
  qa = image.select('pixel_qa');
  mask = qa.bitwiseAnd(cloudsconfidenceBitMask).eq(0) and (qa.bitwiseAnd(cloudsBitMask).eq(0)) or (qa.bitwiseAnd(cloudShadowBitMask).eq(0));
  return image.updateMask(mask);



# + [markdown] id="96dfvXAsNIFn"
# ### A. Vegetation Indices (30m) 
#
# 1) MONTHLY NDVI, NDWI, NDMI
#
# 2) NDMI DELTA
#
# 3) MONTHLY EVI
#
# 4) ANNUAL NDVI, NDWI, NDMI
#
# 5) ANNUAL EVI
#
#
#

# + [markdown] id="7GutnSqkF8ga"
# Landsat 8 Mask

# + [markdown] id="rcuMJ1ClBQeU"
# #### [1] MONTHLY NDVI, NDWI, NDMI

# + id="9WbHgSEfOx2C"
from datetime import datetime, timedelta

def getStartEndDates(yr, mth):
  date_time_obj = datetime.strptime(str(yr) + "-" + str(mth), '%Y-%m') # create datetime object
  d = date_time_obj - timedelta(days=60) # subtract 2 months
  newyr = str(d.year)
  newmth = str(d.month)
  startDate = newyr + "-" + newmth + "-01"
  endDate = str(yr) + "-" + str(mth) + "-01"
  return startDate, endDate


# + id="vBI6DU_4BA3E"
def getMonthlyLandsat(yr, mth, LandsatCollection, Cloudmask):
  startDate, endDate = getStartEndDates(yr, mth)

  LandsatImage = ee.Image(ee.ImageCollection(LandsatCollection)
                  .filterDate(startDate, endDate)    
                  .sort('CLOUD_COVER')
                  .filterBounds(CA_Extent)
                  .map(Cloudmask)
                  .median());

  LandsatClip = LandsatImage.clip(CA_Extent)
  return LandsatClip 


def computeMonthlyVegIndices(yr, mth, LandsatCollection, Cloudmask, EVILayerName, swir_band, nir_band, red_band):
  Landsat = getMonthlyLandsat(yr, mth, LandsatCollection, Cloudmask)
  
  red = Landsat.select(red_band)
  nir = Landsat.select(nir_band)
  swir = Landsat.select(swir_band)

  NDVI = nir.subtract(red).divide(nir.add(red)).rename('NDVI'); # nir-red 
  NDMI = nir.subtract(swir).divide(nir.add(swir)).rename('NDMI'); # nir-swir
  NDMI_Delta = computeNDMIDelta(NDMI, mth, swir_band, nir_band)
  EVI = getMonthlylEVI(yr, mth, EVILayerName)

  Landsat = Landsat.addBands(NDMI).addBands(EVI).addBands(NDVI).addBands(NDMI_Delta)
  Landsat = Landsat.select(['NDVI', 'NDMI', 'EVI', 'NDMIDelta'])
  return (Landsat)


# + [markdown] id="lT1ZATbXCXHg"
# #### [2] NDMI DELTA

# + id="VqCaFB35CbM-"
def computeNDMIDelta(NDMI, mth, swir_band, nir_band):
  baseLine = ee.Image(ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
                  .filter(ee.Filter.calendarRange(2012 - 5, 2012,'year'))
                  .filter(ee.Filter.calendarRange(int(mth),int(mth),'month'))
                  .map(cloudMaskL457)
                  .filterBounds(CA_Extent)
                  .median());                                                                                         
                                                                                            
  baseLineClip = baseLine.clip(CA_Extent)

  swir = baseLineClip.select(swir_band);
  nir = baseLineClip.select(nir_band);
  
  baseLineNDMI = nir.subtract(swir).divide(nir.add(swir)).rename('NDMI_Base');  
  NDMIDelta = NDMI.select('NDMI').subtract(baseLineNDMI.select('NDMI_Base')).rename('NDMIDelta')      
  return NDMIDelta


# + [markdown] id="9jI7CUR3CAPY"
# #### [3] MONTHLY EVI

# + id="FPZcCEkQCBhT"
def getMonthlylEVI(yr, mth, EVILayerName):   #'LANDSAT/LC08/C01/T1_32DAY_EVI'
  startDate, endDate = getStartEndDates(yr, mth)
  EVI_monthly_collection = ee.Image(ee.ImageCollection(EVILayerName)
                          .filterDate(startDate, endDate)
                          .filterBounds(CA_Extent)
                          .select('EVI')
                          .median());
                                  
  EVI_Monthly = EVI_monthly_collection.clip(CA_Extent).rename('EVI');  
  return EVI_Monthly    


# + [markdown] id="7pmAuzaSGBYf"
# #### [4] ANNUAL NDVI, NDWI, NDMI

# + id="GNUUAx9dGGXL"
def getAnnualLandsat(yr, LandsatCollection, Cloudmask):
  start_yr = int(yr) - 1 
  Landsat_annual_Image = ee.Image(ee.ImageCollection(LandsatCollection)
                  .sort('CLOUD_COVER')
                  .filter(ee.Filter.date(f'{start_yr}-10-01', f'{yr}-04-30'))
                  .filterBounds(CA_Extent)
                  .map(Cloudmask)
                  .max());

  Landsat_Annual_Clip = Landsat_annual_Image.clip(CA_Extent)
  return Landsat_Annual_Clip


# + id="SgPleTXjxma4"
def computeAnnualVegIndices(yr, LandsatCollection, Cloudmask, EVILayerName, swir_band, nir_band, red_band):
  Landsat = getAnnualLandsat(yr, LandsatCollection, Cloudmask)

  red = Landsat.select(red_band)
  nir = Landsat.select(nir_band)
  swir = Landsat.select(swir_band)

  NDVI = nir.subtract(red).divide(nir.add(red)).rename('Annual_NDVI');
  NDMI = nir.subtract(swir).divide(nir.add(swir)).rename('Annual_NDMI');
  EVI = getAnnualEVI(yr, EVILayerName)

  Landsat = Landsat.addBands(NDVI).addBands(NDMI).addBands(EVI)
  Landsat = Landsat.select(['Annual_NDVI', 'Annual_NDMI', 'Annual_EVI'])
  return (Landsat)


# + [markdown] id="fiiozfa6GHrs"
# #### [5] ANNUAL EVI

# + id="rMx1jCkgGJOv"
def getAnnualEVI(yr, EVILayerName): #'LANDSAT/LC08/C01/T1_32DAY_EVI'
  start_yr = int(yr) - 1 
  EVI_annual_collection = ee.Image(ee.ImageCollection(EVILayerName)
                          .filter(ee.Filter.date(f'{start_yr}-10-01', f'{yr}-04-30'))
                          .filterBounds(CA_Extent)
                          .select('EVI')
                          .max());
                                  
  EVI_Annual = EVI_annual_collection.clip(CA_Extent).rename('Annual_EVI');   
  return EVI_Annual


# + [markdown] id="bh-qjAvDNfXi"
# ### B. Daily Weather (resample to 30m) 
#
# 1) Gridmet 
#
# 2) Gridmet FFWI
#
# 3) Daymet 
#
# 4) Daymet VPD & RH
#

# + [markdown] id="gVwWuugF_SGw"
# #### [1] GridMET

# + id="e3VZHx9MNr_6"
def getGRIDMET(yr, jd):
  gridmet_layers = ['pr', 'rmax', 'rmin', 'sph', 'srad', 'th', 
                    'tmmn','tmmx', 'vs', 'erc', 'bi', 'eto',
                    'fm100', 'fm1000', 'etr', 'vpd']

  gridmet_collection = ee.Image(ee.ImageCollection('IDAHO_EPSCOR/GRIDMET')
                            .filter(ee.Filter.calendarRange(int(yr), int(yr),'year'))
                            .filter(ee.Filter.calendarRange(int(jd),int(jd),'DAY_OF_YEAR'))
                            .select(gridmet_layers)
                            .mosaic())

  gridmet = gridmet_collection.clip(CA_Extent) 

  gridmet = getFFWI(gridmet)
  return gridmet


# + [markdown] id="NniVRG4MES0a"
# #### [2] FFWI

# + id="bP6NQ-vJERww"
def getFFWI(gridmet):
  FFWI = gridmet.expression(
          "(b('rmax') < 10) ? (1 - (2 *( (0.03229 + (0.281073 * b('rmax')) - (0.000578 * b('rmax') * b('tmmx'))) / (30 * 1.0))) + ((1.5 * ( (0.03229 + (0.281073 * b('rmax')) - (0.000578 * b('rmax') * b('tmmx'))) /(30*1.0))**2)) - (0.5 * ((0.03229 + (0.281073 * b('rmax')) - (0.000578 * b('rmax') * b('tmmx'))) / (30 * 1.0))**3)) * (((1+b('vs')**2)**0.5) / 0.3002)" +
          ": (b('rmax') > 10 && b('rmax') <= 50) ? (1 - (2 *( (2.22749 + (0.160107* b('rmax')) - (0.01478* b('tmmx'))) / (30 * 1.0))) + ((1.5 * ( (2.22749 + (0.160107* b('rmax')) - (0.01478* b('tmmx'))) /(30*1.0))**2)) - (0.5 * ((2.22749 + (0.160107* b('rmax')) - (0.01478* b('tmmx'))) / (30 * 1.0))**3)) * (((1+b('vs')**2)**0.5) / 0.3002)" +
          ": (b('rmax') > 50) ? (1 - (2 *((21.0606 + (0.005565* (b('rmax')**2)) - (0.00035 * b('rmax') * b('tmmx')) - (0.483199 * b('rmax'))) / (30 * 1.0))) + ((1.5 * ( (21.0606 + (0.005565* (b('rmax')**2)) - (0.00035 * b('rmax') * b('tmmx')) - (0.483199 * b('rmax'))) /(30*1.0))**2)) - (0.5 * ((21.0606 + (0.005565* (b('rmax')**2)) - (0.00035 * b('rmax') * b('tmmx')) - (0.483199 * b('rmax'))) / (30 * 1.0))**3)) * (((1+b('vs')**2)**0.5) / 0.3002)" +
          ": 0").rename('FFWI')

  gridmet_full = gridmet.addBands(FFWI)
  return gridmet_full


# + [markdown] id="Z3lLhrui_V68"
# #### [3] DayMET

# + id="CxwYWYGB_X5I"
def getDAYMET(yr, jd):
  daymet_layers = ['tmax', 'srad', 'vp', 'tmin', 'prcp']

  daymet_collection = ee.Image(ee.ImageCollection('NASA/ORNL/DAYMET_V3')
                            .filter(ee.Filter.calendarRange(int(yr), int(yr),'year'))
                            .filter(ee.Filter.calendarRange(int(jd),int(jd),'DAY_OF_YEAR'))
                            .select(daymet_layers)
                            .mosaic())

  daymet = daymet_collection.clip(CA_Extent)

  daymet = getVPD_RH(daymet)

  return daymet 


# + [markdown] id="M_q5SNKR_xtz"
# #### [4] VPD & RH

# + id="_j3t0rM5_0wI"
def getVPD_RH(daymet):
  VP_amb = daymet.expression("610.78**((17.269 * tmin) / (237.3 + tmin))", {
    'tmin': daymet.select('tmin')
  })
    
  VP_sat = daymet.expression("610.78**((17.269 * (0.394 * tmin + 0.606 * tmax)) / (237.3 + (0.394 * tmin + 0.606 * tmax)))", {
    'tmin': daymet.select('tmin'),
    'tmax': daymet.select('tmax')
  })

  VP_Deficit = VP_sat.subtract(VP_amb).divide(1000).rename('DayMET_VPD')
  Rel_Humidity = VP_amb.divide(VP_sat).multiply(100).rename('DayMET_RH')

  daymet_full = daymet.addBands(VP_Deficit).addBands(Rel_Humidity);
  return daymet_full


# + [markdown] id="NJJ4hKGWN3fx"
# ### C. Topography (resample to 30m)
#
# 1) DEM, Slope, Aspect
#
# 2) Multiscale TPI 
#
# 3) Topo Diversity 

# + [markdown] id="8w7LVLleeQYR"
# #### [1] DEM, Slope, Aspect

# + id="Lqri41ALeTIO"
DEM_Image = ee.Image('USGS/NED').select('elevation');

DEM = DEM_Image.clip(CA_Extent).rename('DEM')
Slope = ee.Terrain.slope(DEM).rename('Slope')
Aspect = ee.Terrain.aspect(DEM).rename('Aspect')

topo = DEM.addBands(Slope).addBands(Aspect);

# + [markdown] id="qSxnUOiwdo8q"
# #### [2] TPI 

# + id="jQG8qoBtOS1g"
# 270 m (Valley (-) Ridge (+))
TPI_Image = ee.Image('CSP/ERGo/1_0/US/mTPI').select('elevation');

TPI = TPI_Image.clip(CA_Extent) 

# ridge > + 1 STDEV
# 2 upper slope > 0.5 STDV =< 1 STDV
# 3 middle slope> -0.5 STDV, < 0.5 STDV, slope > 5 deg
# 4 flats slope >= -0.5 STDV, =< 0.5 STDV , slope <= 5 deg
# 5 lower slopes >= -1.0 STDEV, < 0.5 STDV
# 6 valleys < -1.0 STDV

# + id="IeRCUSiV3GEk"
def getStats(operation, ras_Image, band_name):
  out = ras_Image.reduceRegion(**{
      'reducer': operation, 
      'bestEffort': True,
      'geometry': CA_Extent.geometry()
      });
  return out.getNumber(f'{band_name}').getInfo()


TPI_ST = getStats(ee.Reducer.stdDev(), TPI, 'elevation')
TPI_mean = getStats(ee.Reducer.mean(), TPI, 'elevation')
TPI_max = getStats(ee.Reducer.max(), TPI, 'elevation')
TPI_min = getStats(ee.Reducer.min(), TPI, 'elevation')

# + id="bhRSkYis1WD4"
thresholds = ee.Image([TPI_min, TPI_mean - TPI_ST, 
                       TPI_mean - TPI_ST/2, TPI_mean + TPI_ST/2, TPI_mean + TPI_ST, TPI_max]);

TPI_Zones = TPI.expression(f'b(0) <= {TPI_mean - (2*(TPI_ST))} ? 1 : \
                             b(0) < {TPI_mean - (1/2*(TPI_ST))} ? 2 : \
                             b(0) < {TPI_mean + (1/2*(TPI_ST))} ? 3 : \
                             b(0) < {TPI_mean + (1*(TPI_ST))} ? 4 : \
                             b(0) < {TPI_max} ? 5 : 0').rename('TPI')

ridge = TPI_Zones.select(['TPI']).eq(1)
valley = TPI_Zones.select(['TPI']).eq(5)

# + [markdown] id="suef83RieTW5"
# ####[3] Topo Diversity

# + id="N9g0rfzDegms"
topo_diversity_Image = ee.Image('CSP/ERGo/1_0/Global/ALOS_topoDiversity').select('constant');

topo_diversity= topo_diversity_Image.clip(CA_Extent).rename('TOPO_DIV')
topo = topo.addBands(topo_diversity)

# + [markdown] id="fBSbuDQVOTA9"
# ### D. LandFire (30m)
#
# USE: ee.Reducer.countDistinct()
#
# 1) Vegetation Type 
#
# 2) Vegetation Height 

# + [markdown] id="_aMhewxPhgdq"
# #### [1] Vegetation Type

# + id="VnRzIEcoOixL"
EVT_Image = ee.Image(ee.ImageCollection('LANDFIRE/Vegetation/EVT/v1_4_0').select('EVT').mosaic())

EVT = EVT_Image.clip(CA_Extent).rename('EVT')

# + id="tHWLZGD7qjZM" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601589651882, "user_tz": 420, "elapsed": 1166, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="e34fadfd-883f-49ff-e558-bfc0e4873bbd"
print(EVT.bandNames().getInfo())

# + [markdown] id="xZnubPeOhjxm"
# #### [2] Vegetation Height

# + id="w6WMXYDUOsLO"
EVH_Image = ee.Image(ee.ImageCollection('LANDFIRE/Vegetation/EVH/v1_4_0').select('EVH').mosaic())

EVH = EVH_Image.clip(CA_Extent).rename('EVH')

# + id="a8QGDIlJ17iJ"
Veg_TPI = TPI_Zones.addBands(EVH).addBands(EVT);

# + [markdown] id="CRS-vEt8QcHx"
# ### E. Cost Distance to human settlement (resample to 30m)

# + id="jv5f03MGQe6t"
roads = ee.FeatureCollection('TIGER/2016/Roads');

# // Create a source image (palm oil mill=1, others=0)
sources = ee.Image().toByte().paint(roads, 1);

# // Mask the sources image with itself.
sources = sources.updateMask(sources);
sources = sources.clip(CA_Extent)

# // The cost data is generated from slope 
cover = Slope.select(0);

# // Variable cost
cost = cover.lte(10) and (cover.gt(0)).multiply(1).add(
       cover.lte(30) and (cover.gt(10)).multiply(2)).add(
       cover.lte(50) and (cover.gt(30)).multiply(3)).add(
       cover.lte(60) and (cover.gt(50)).multiply(4)).add(
       cover.lte(100) and (cover.gt(60)).multiply(5))

# // Compute the cumulative cost to traverse the land cover.
cumulativeCost = cost.cumulativeCost(**{
  'source': sources,
  'maxDistance': 10 * 1000 
}).rename('travelCost')


# + [markdown] id="rvRIcMy5PEvZ"
# ## For each instance of fire, GET layers, Mosaic, Reduce, Save
#
# Make sure all layers are aligned, with correct CRS projection and resolution

# + id="LKMhZ8JeiNUV"
# import feature collection asset 
fire_perimeters = ee.FeatureCollection('users/escaduto/CA_FireTrends/DerivedPerimeters_2012_2020')
CA_Extent = ee.FeatureCollection('users/escaduto/CA_FireTrends/CA_Extent')

# + id="AYUjurEhiNUa"
outpath = os.path.join(r'Products', 'ALL_2012_2020')
Interpolated_CA = gpd.read_file(os.path.join(outpath, f'ALL_2012_2020.geojson'))

# + id="Uxtig6HIiNUd" colab={"base_uri": "https://localhost:8080/", "height": 102} executionInfo={"status": "ok", "timestamp": 1601589672094, "user_tz": 420, "elapsed": 4117, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="00264785-06da-4251-c035-0d5b3793e6da"
Interpolated_CA.columns


# + id="fvNxmwhIOmmw"
def exportGEE(features, var, yrjd):
  exporttask = ee.batch.Export.table.toDrive(**{
  'collection': features,
  'folder': "GoogleEE_VariableExtraction_10012020",
  'fileNamePrefix': f'{yrjd}_{var}',
  # 'selectors':(["Aspect","DEM","Slope"]),
  });

  exporttask.start()
  print(f'exporting {yrjd}_{var}.csv')



# + id="edc3YEm5O8bv"
def GEEReducer(IMG, collection, reducer, scale):
  reduced = IMG.reduceRegions(**{
    'collection': collection,
    'reducer': reducer,
    'scale': scale,
  })
  return reduced


# + id="6SjHfXhOiPan"
# 'LANDSAT/LE07/C01/T1_SR'
# 'LANDSAT/LC08/C01/T1_SR'

def zonalStatsByDay(yr_jd_list, LandsatCollection, Cloudmask, EVILayerName, swir_band, nir_band, red_band):
  for yrjd in yr_jd_list:
    yr_jd_dt = datetime.strptime(yrjd, '%Y%j') 
    year = yr_jd_dt.year
    dy = yr_jd_dt.day
    mnth = yr_jd_dt.month
    mnth = str(mnth).zfill(2)
    jd = yrjd[-3:]
    print(jd, mnth, year)

    day_perimeter = fire_perimeters.filter(ee.Filter.eq('YRJD', yrjd))
    
    # Weather 
    daymet = getDAYMET(year, jd)
    gridmet = getGRIDMET(year, jd)

    # Vegetation 
    monthly_indices = computeMonthlyVegIndices(year, mnth, LandsatCollection, Cloudmask, EVILayerName, swir_band, nir_band, red_band)
    annual_indices = computeAnnualVegIndices(year, LandsatCollection, Cloudmask, EVILayerName, swir_band, nir_band, red_band)

    allLayers = monthly_indices.addBands(annual_indices).addBands(gridmet).addBands(topo).addBands(cumulativeCost)
    all_Features = GEEReducer(allLayers, day_perimeter, ee.Reducer.mean(), 10)

    exportGEE(all_Features, 'dailyfire_zonalstats', yrjd)



# + [markdown] id="LJk0iouPyanj"
# #### Daymet 2012-2019

# + id="hnpl1d4yyFSQ"
def daymet_zonalStatsByDay(yr_jd_list):
  for yrjd in yr_jd_list:
    yr_jd_dt = datetime.strptime(yrjd, '%Y%j') 
    year = yr_jd_dt.year
    dy = yr_jd_dt.day
    mnth = yr_jd_dt.month
    mnth = str(mnth).zfill(2)
    jd = yrjd[-3:]
    print(jd, mnth, year)

    day_perimeter = fire_perimeters.filter(ee.Filter.eq('YRJD', yrjd))
    
    # Weather 
    daymet = getDAYMET(year, jd)
    daymet_Features = GEEReducer(daymet, day_perimeter, ee.Reducer.mean(), 10)

    exportGEE(daymet_Features, 'DAYMET', yrjd)


# + id="PjsjIjq4yqP0"
daymet_fires = Interpolated_CA[Interpolated_CA['Year'] <= 2019]
yr_jd_list  = daymet_fires['YRJD'].unique().tolist()
yr_jd_list.sort(reverse=True)

# + id="sM9zCYxuymdu"
daymet_zonalStatsByDay(yr_jd_list[158:])

# + [markdown] id="9ubEUklIzI_l"
# #### Count Distinct (TPI, EVH, EVT) [**ARCPY for now]

# + id="p_F_C72AZtGb"
#[TODO] Fuel Type, Height, TPI, Fuel Model 
# LandFire + TPI (Percent by unique value)
#Veg_TPI_Features = GEEReducer(TPI_Zones, day_perimeter, ee.Reducer.countDistinct(), 30)
#ee.Reducer.group(**{'reducer': ee.Reducer.count(),'groupField':0,'groupName':'TPI'})

# exportGEE(Veg_TPI_Features, 'Veg_TPI', yrjd)

# + [markdown] id="23rRYylFVfSu"
# #### Landsat 8 2013-2020

# + id="ufEtbPJmQL60"
L8_fires = Interpolated_CA[Interpolated_CA['Year'] >= 2013]
yr_jd_list  = L8_fires['YRJD'].unique().tolist()
yr_jd_list.sort(reverse=True)

# + id="dyfx3an7VmVd"
LandsatCollection = 'LANDSAT/LC08/C01/T1_SR'
Cloudmask = maskL8sr
EVILayerName='LANDSAT/LC08/C01/T1_32DAY_EVI'
swir_band = 'B6'
nir_band = 'B5'
red_band = 'B4'

zonalStatsByDay(yr_jd_list, LandsatCollection, Cloudmask, EVILayerName, swir_band, nir_band, red_band)

# + [markdown] id="u69laz5RWB1B"
# #### Landsat 7 2012-2013

# + id="vJqcrg92V2u_"
L7_fires = Interpolated_CA[Interpolated_CA['Year'] < 2013]
yr_jd_list = L7_fires['YRJD'].unique().tolist()
yr_jd_list.sort(reverse=True)

# + id="MEODZmm3V73C"
LandsatCollection = 'LANDSAT/LE07/C01/T1_SR'
Cloudmask = cloudMaskL457
EVILayerName= 'LANDSAT/LE07/C01/T1_32DAY_EVI'
swir_band = 'B5'
nir_band = 'B4'
red_band = 'B3'

zonalStatsByDay(yr_jd_list, LandsatCollection, Cloudmask, EVILayerName, swir_band, nir_band, red_band)

# + id="CSH0MspRwOj5" colab={"base_uri": "https://localhost:8080/", "height": 35} executionInfo={"status": "ok", "timestamp": 1600550551736, "user_tz": 420, "elapsed": 549, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="5e136157-91c6-4e33-828c-1f5f948c7be7"
os.getcwd()

# + [markdown] id="iBH_YhyhWr1y"
# ## Visualize Layers

# + id="lRtrJPozi7d-" colab={"base_uri": "https://localhost:8080/", "height": 757} executionInfo={"status": "ok", "timestamp": 1600651899067, "user_tz": 420, "elapsed": 1912, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="86fdb9f3-378d-49db-dc42-19e7b117b882"
visParams = {
  'bands': ['B4', 'B3', 'B2'],
  'min': 0,
  'max': 3000,
  'gamma': 1.4,
};

ndviParams ={
  'min': -1,
  'max': 1,
  'palette': ['ffff00', '330033','7fffd4']
};

Map = emap.Map(center=[38.50178453635526,-122.74843617724784], zoom=10)
# Map.addLayer(TPI_Zones, {'palette': ['ffff00', '330033','7fffd4', 
#                                                          'b99879','cc8e7f', '268b07']}, 'TPI_class')
Map.addLayer(landsat, visParams , 'landsat')
# Map.addLayer(valley, {'min':0, 'max':1, 'palette': ['ffff00', '330033']} , 'valley')
# Map.addLayer(ridge, {'min':0, 'max':1, 'palette': ['ffff00', '330033']} , 'ridge')
Map.addLayer(cumulativeCost, {'min': 0, 'max': 5e3, 'palette': ['34a853', 'fbc034', 'ea4335', '4285f4']} , 'cumulativeCost')
Map.addLayerControl()
Map

# + [markdown] id="iSj8q-bwz19Z"
# ## Post-Processing: CSV
#
# - By year 
# - All 
#

# + [markdown] id="APrklNLwDz0u"
# ### Merge daily files

# + id="W5vj37QqUs5B"
rootPath = "/content/drive/My Drive/GoogleEE_VariableExtraction_10012020"

# + id="6wXCumI1z4uU"
import glob, os

# gridmet, indices, topo, cost distance
vars_csv = glob.glob(os.path.join(rootPath, '*dailyfire_zonalstats.csv'))
combined_vars = pd.concat([pd.read_csv(f) for f in vars_csv])

# + id="4-gSJdOFWvJc"
outpth_all = 'Variables/CSVs'
combined_vars.to_csv(os.path.join(outpth_all, "combinedvars_2012_2020.csv"))

# + id="doRW9DOkTZ1j"
# daymet 
daymet_csv = glob.glob(os.path.join(rootPath, '*_DAYMET.csv'))
combined_daymet = pd.concat([pd.read_csv(f) for f in daymet_csv])
combined_daymet.to_csv(os.path.join(outpth_all, "daymet_2012_2019.csv"))

# + [markdown] id="PGa24zQ_8I99"
# Percent Count do with arcpy

# + id="bNIdc6RJInRS"
import os
# TPI, EVH, EVT, fuel model (convert to decimal ratio)
outpath = os.path.join(r'Products', 'ALL_2012_2020')
Interpolated_CA = gpd.read_file(os.path.join(outpath, f'ALL_2012_2020.geojson'))

# + id="DnGZ3JzS8NnT"
Interpolated_CA_4326 = Interpolated_CA.to_crs('EPSG:4326')

# + id="bDm9jcwoMSqW"
Interpolated_CA_4326['Group_ID'] = Interpolated_CA_4326['Unique_ID'].astype('str').str.zfill(8)

# + id="TG87V1IP91DS" colab={"base_uri": "https://localhost:8080/", "height": 71} executionInfo={"status": "ok", "timestamp": 1601659878208, "user_tz": 420, "elapsed": 34645, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="f57aa152-cd01-488a-99e7-7b82fd1e0dea"
Interpolated_CA_4326.to_file(os.path.join(outpath, f'ALL_2012_2020_WGS84.shp'))
Interpolated_CA_4326.to_file(os.path.join(outpath, f'ALL_2012_2020_WGS84.geojson'), driver='GeoJSON')

# + [markdown] id="cE_2wlzFD9DI"
# ### Cleanup & Compute Ratios (TPI, FBFM, EVT, EVH)

# + [markdown] id="nHfzKD9V-chU"
# #### CBD

# + id="9t-mMvG9-e_-"
CBD_df = pd.read_csv('Variables/ExtractVar_CSV/CBD/CBD_Extracted.csv', index_col=0).reset_index()
CBD_df.columns = ['OID_', 'Group_ID', 'ZONE_CODE', 'COUNT', 'AREA', 'CBD']
CBD_df = CBD_df[['Group_ID', 'CBD']]
CBD_df['Unique_ID'] = CBD_df['Group_ID'].astype('int64')
CBD_df = CBD_df.drop(columns=['Group_ID'])

# + id="DfSvRM2IT5wI" colab={"base_uri": "https://localhost:8080/", "height": 68} executionInfo={"status": "ok", "timestamp": 1601662732366, "user_tz": 420, "elapsed": 392, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="75cbaa6b-c552-4149-c1ea-dc49b6d1a07d"
CBD_df.dtypes

# + [markdown] id="umeuNE1HBo7B"
# #### TPI

# + id="kbH72d8J-fFZ"
TPI_df = pd.read_csv('Variables/ExtractVar_CSV/TPI/TPI_Extracted.csv', index_col=0).reset_index()
TPI_df.columns = ['OID_', 'Group_ID', 'Valley', 'LowerSlope', 'MidSlope', 'UpperSlope', 'Ridge']
TPI_df = TPI_df[['Group_ID', 'Valley', 'LowerSlope', 'MidSlope', 'UpperSlope', 'Ridge']]
TPI_df['Unique_ID'] = TPI_df['Group_ID'].astype('int64')
TPI_df = TPI_df.drop(columns=['Group_ID'])

# + id="8MIjB6UGDsyE"
# Percent Ratio 
TPI_df = TPI_df.set_index('Unique_ID')
TPI_df = TPI_df.div(TPI_df.sum(axis=1), axis=0)
TPI_df = TPI_df.reset_index()

# + id="7i1o4jilT7HY" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601662736135, "user_tz": 420, "elapsed": 295, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="db548c1a-c787-48f9-bda9-1816b53614db"
len(TPI_df)

# + [markdown] id="Hm9Cd1SyCo_O"
# #### FBFM

# + id="Uvm-IlwYCAYJ"
FBFM_df = pd.read_csv('Variables/ExtractVar_CSV/FBFM/FBFM_Extracted.csv', index_col=0).reset_index()
FBFM_df = FBFM_df.drop(columns=['OID_', 'URBAN',	'SNOW_ICE',	'AGRICULTUR',	'WATER',	'BARREN']) 
FBFM_df['Unique_ID'] = FBFM_df['GROUP_ID'].astype('int64')
FBFM_df = FBFM_df.drop(columns=['GROUP_ID'])

# + id="VDkqj5QuCT_E"
# Percent Ratio 
FBFM_df = FBFM_df.set_index('Unique_ID')
FBFM_df = FBFM_df.div(FBFM_df.sum(axis=1), axis=0)
FBFM_df = FBFM_df.reset_index()

# + id="NyhPiSPSeMpB" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601662740921, "user_tz": 420, "elapsed": 325, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="0fd69324-f50e-4351-dde3-5bf9af7d88bb"
len(FBFM_df)

# + [markdown] id="seaPGUSsCz4M"
# #### EVH

# + id="S-nEBMNXC3Z0"
EVH_df = pd.read_csv('Variables/ExtractVar_CSV//EVH/EVH_Extracted.csv', index_col=0).reset_index()
EVH_df = EVH_df.drop(columns=['OID_']) 
EVH_df['Unique_ID'] = EVH_df['GROUP_ID'].astype('int64')
EVH_df = EVH_df.drop(columns=['GROUP_ID'])

# + id="r5bcJRIKDq_e"
# Percent Ratio 
EVH_df['shortHeight'] = EVH_df['VALUE_101'] + EVH_df['VALUE_104']
EVH_df['medShortHeight'] = EVH_df['VALUE_102'] + EVH_df['VALUE_105']
EVH_df['medHeight'] = EVH_df['VALUE_103'] + EVH_df['VALUE_104'] + EVH_df['VALUE_107']
EVH_df['TreeHeight_5m'] = EVH_df['VALUE_108']
EVH_df['TreeHeight_10m'] = EVH_df['VALUE_109']
EVH_df['TreeHeight_25m'] = EVH_df['VALUE_110']
EVH_df['TreeHeight_50m'] = EVH_df['VALUE_111']
EVH_df['TreeHeight_50m+'] = EVH_df['VALUE_112']

EVH_df = EVH_df[['Unique_ID', 'shortHeight', 'medShortHeight', 'medHeight', 'TreeHeight_5m', 
        'TreeHeight_10m', 'TreeHeight_25m', 'TreeHeight_50m', 'TreeHeight_50m+']]

# + id="WyA8hfdKgilz"
# Percent Ratio 
EVH_df = EVH_df.set_index('Unique_ID')
EVH_df = EVH_df.div(EVH_df.sum(axis=1), axis=0)
EVH_df = EVH_df.reset_index()

# + id="5tYfHWWIgpRO" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601662749130, "user_tz": 420, "elapsed": 497, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="058e1ec3-d02b-4f54-eb32-a78eb2b6dfd3"
len(EVH_df)

# + [markdown] id="AMndlmuRCz1X"
# #### EVT

# + id="FyKWsujRC2xS"
EVT_df = pd.read_csv('Variables/ExtractVar_CSV/EVT/EVT_Extracted.csv', index_col=0).reset_index()
EVT_meta = pd.read_csv('Variables/ExtractVar_CSV/EVT/US_130EVT_02092015.csv', index_col=0).reset_index()
EVT_df['Unique_ID'] = EVT_df['GROUP_ID'].astype('int64')
EVT_df = EVT_df.drop(columns=['OID_', 'GROUP_ID'])

# + id="tIgnkt3whCe8"
EVT_meta = EVT_meta[['VALUE', 'EVT_PHYS', 'EVT_CLASS']]
EVT_meta['EVT_PHYS'] = np.where(EVT_meta['EVT_PHYS'] == 'Agricultural', EVT_meta['EVT_CLASS'], EVT_meta['EVT_PHYS'])
EVT_meta = EVT_meta[['VALUE', 'EVT_PHYS']]

# + id="TFALaZ4_ijXC"
EVT_meta = EVT_meta.replace(['Conifer', 'Conifer-Hardwood', 'Open tree canopy'], 'Conifer')
EVT_meta = EVT_meta.replace(['Hardwood', 'Hardwood-Conifer'], 'Hardwood')
EVT_meta = EVT_meta.replace(['Barren', 'Non-vegetated'], 'Barren')
EVT_meta = EVT_meta.replace(['Barren', 'Non-vegetated'], 'Barren')
EVT_meta = EVT_meta.replace(['Grassland', 'Herbaceous - grassland', 'Exotic Herbaceous', 'Herbaceous / Nonvascular-dominated'], 'Grassland')
EVT_meta = EVT_meta.replace(['Shrubland', 'Exotic Tree-Shrub'], 'Shrubland')
EVT_meta = EVT_meta.replace(['Riparian'], 'Riparian')
EVT_meta = EVT_meta.replace(['Sparsely Vegetated'], 'Sparsely Vegetated')
EVT_meta = EVT_meta.replace(['Snow-Ice', 'Nodata', 'Open Water', 'Quarries-Strip Mines-Gravel Pits', 'Other'], 'Other')
EVT_meta = EVT_meta.replace(['Developed', 'Developed-High Intensity', 'Developed-Low Intensity', 'Developed-Medium Intensity', 'Developed-Roads'], 'Developed')

evt_value_info = EVT_meta.groupby(['EVT_PHYS'])['VALUE'].apply(list).reset_index()


# + id="wcArYooCQ2oF"
def prepend(lst, str): 
    # Using format() 
    str += '{0}'
    lst = [str.format(i) for i in lst] 
    return(lst) 


# + id="4PUCSmPoRNv4"
column_list = EVT_df.columns.tolist()

def updateEVTColumns(column_list, evt_value_info, EVT_df):
  for index, row in evt_value_info.iterrows():
    ID_list = prepend(row.VALUE, 'VALUE_')
    ID_list = [x for x in ID_list if x in column_list]
    EVT_df[row.EVT_PHYS] = EVT_df[ID_list].sum(axis=1)
  return EVT_df


# + id="AlpXdWXjXxb2"
list_col = evt_value_info['EVT_PHYS'].tolist()
list_col.append('Unique_ID')

# + id="qri-D6FtWuIs"
EVT_df = updateEVTColumns(column_list, evt_value_info, EVT_df)
EVT_df = EVT_df[list_col]

# + id="jj-zDe49Xab2"
# Percent Ratio 
EVT_df = EVT_df.set_index('Unique_ID')
EVT_df = EVT_df.div(EVT_df.sum(axis=1), axis=0)
EVT_df = EVT_df.reset_index()

# + id="AghzrLJzWDh4" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601664534232, "user_tz": 420, "elapsed": 285, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="3ee00389-8e91-4bbf-b12a-ea41ad6e0b7d"
len(EVT_df)

# + [markdown] id="NXUWwm-ODpTa"
# ### Merge ALL

# + [markdown] id="ZfUvYaCaGWjt"
# ####2012-2020 (w/out Daymet)

# + id="EyWru9ncEMYK"
vars_df = pd.read_csv('Variables/CSVs/combinedvars_2012_2020.csv', index_col=0).reset_index()

# + id="56uzyb39ZDPz"
import pandas as pd
from functools import reduce

df_list = [vars_df, CBD_df, FBFM_df, TPI_df, EVH_df, EVT_df]

gridmet_2012_2020 = reduce(lambda  left,right: pd.merge(left,right,on=['Unique_ID'],
                                            how='left'), df_list)

# + id="yVksFKl7JBm3" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601664566010, "user_tz": 420, "elapsed": 350, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="0dc99db0-b938-44cc-c3d9-a196989696da"
gridmet_2012_2020['FireID'].nunique()
# 408

# + id="jSdRFOpDpesS"
gridmet_2012_2020 = gridmet_2012_2020.fillna(0)
gridmet_2012_2020 = gridmet_2012_2020.drop(columns=['index', 'system:index', 'Area (ac_1', 'DRAWSEQ'])

# + id="jcUC4JxndGxX"
gridmet_2012_2020.to_csv('Variables/Final_Variables/gridmet_2012_2020.csv')

# + [markdown] id="7e-rY4i8GaPY"
# ####2012-2019 (w/Daymet)

# + id="SIqkDZNFGbff"
daymet_df = pd.read_csv('Variables/CSVs/daymet_2012_2019.csv', index_col=0).reset_index()
#len 361

# + id="3-ChvhHdHkUY"
daymet_df = daymet_df[['DayMET_RH', 'DayMET_VPD', 
            'Unique_ID', 'FireID', 'prcp', 'srad',
            'tmax', 'tmin', 'vp']]

# + id="NNa5CyI3HgCa"
df_list = [daymet_df, CBD_df, FBFM_df, TPI_df, EVH_df, EVT_df]

daymet_2012_2019 = reduce(lambda  left,right: pd.merge(left,right,on=['Unique_ID'],
                                            how='inner'), df_list)
#len 359

# + id="ByEeZZK_IukI" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601575330407, "user_tz": 420, "elapsed": 358, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="f8e97797-5e1b-4130-b0c1-83e5062455ea"
len(daymet_2012_2019['FireID'].unique())
# len 351

# + id="ucI16jsSdRYM"
daymet_2012_2019.to_csv('Variables/Final_Variables/daymet_2012_2019.csv')

# + [markdown] id="NDLG7plwGmHp"
# #### w/RAWS

# + id="7Xkf85ir0nyn"
raws_csv_list = listFiles_ByExt(f"Variables/WeatherStation_FireSpread/data/Final_CSV", '.csv')
combined_csv = pd.concat([pd.read_csv(f) for f in raws_csv_list ])

# + id="5d0kkc1RqmSq"
combined_csv

# + id="iGoGE44LB4kY" colab={"base_uri": "https://localhost:8080/", "height": 119} executionInfo={"status": "ok", "timestamp": 1601664877750, "user_tz": 420, "elapsed": 625, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="1793443c-352b-4e88-8ec4-27e12aa57203"
combined_csv = combined_csv[combined_csv['Date'].notna()]
combined_csv['Date'] = combined_csv['Date'].astype('int64')

# + id="3KAxfP0cujEP"
import math
def calculateVPD(Temp_val, RH_val):
  '''
  returns VPD based on temp and rel humidity values
  '''
  VP_sat = 0.61078 * math.exp((Temp_val * 17.2694) / (Temp_val + 238.3))
  VP_air = (VP_sat * RH_val) / 100 
  VPD = VP_sat - VP_air 
  return VPD


# + id="Ntc68FTqwicf"
def convertTemp_FtoC(Temp_val):
  '''
  converts F to C temperature
  '''
  C_TEMP = ((Temp_val - 32) * 5 ) / 9 
  return C_TEMP


# + id="if_PyJLjMQTl"
night_list = []
day_list = ['6 am', '7 am', '8 am', '9 am', '10 am','11 am', '12 pm','1 pm', '2 pm','3 pm', '4 pm', '5 pm', '6 pm',
            '7 pm', '8 pm', '9 pm', '10 pm','11 pm', '12 am','1 am', '2 am','3 am', '4 am', '5 am']

def getDailyRAWSByFire(rim_df):
  '''
  converts hourly RAWS to daily (night vs day). keeps only day.
  '''
  rim_df = rim_df.dropna(subset=['Date'])
  criteria = [rim_df['Hour'].isin(day_list), rim_df['Hour'].isin(night_list)]
  values = ['Day', 'Night']
  rim_df['Day_Night']  = np.select(criteria, values, 0)
  day_raws_df = rim_df.groupby(by=['Day_Night', 'Date', 'Fire']).agg({'Total_Solar_Rad': 'mean', 
                                                                    'Wind_Avg_mph': 'mean',
                                                                    'Wind_Dir_Deg': 'mean',
                                                                    'Wind_Max_mph': ['mean', 'max'],
                                                                    'Air_Temp_Avg': ['mean', 'max'],
                                                                    'Fuel_Temp_Avg': 'mean',
                                                                    'Fuel_Moist_Per': 'mean',
                                                                    'Rel_Hum_Per': 'mean',
                                                                    'Dew_Point_Deg': 'mean',
                                                                    'Wet_Bulb': 'mean',
                                                                    'Total_Precip': 'mean'}).reset_index()
  day_raws_df.columns = ['Day_Night', 'Date', 'Fire', 'Total_Solar_Rad', 'Wind_Avg_mph',
       'Wind_Dir_Deg', 'Wind_Max_Avg', 'Wind_Max_mph', 'Air_Temp_Avg', 'Air_Temp_Max',
       'Fuel_Temp_Avg', 'Fuel_Moist_Per', 'Rel_Hum_Per', 'Dew_Point_Deg',
       'Wet_Bulb', 'Total_Precip']

  day_raws_df['Air_Temp_Celcius'] = convertTemp_FtoC(day_raws_df['Air_Temp_Avg'])
  day_raws_df['Air_Temp_Max_Celcius'] = convertTemp_FtoC(day_raws_df['Air_Temp_Max'])
  day_raws_df['VPD'] = day_raws_df.apply(lambda x: calculateVPD(x['Air_Temp_Celcius'],x['Rel_Hum_Per']),axis=1)
  day_raws_df = day_raws_df[day_raws_df['Day_Night'] == 'Day']
  day_raws_df['Date'] = day_raws_df['Date'].astype('str')
  day_raws_df['Date'] = pd.to_datetime(day_raws_df['Date'], format= '%Y%m%d')
  day_raws_df['Year'] =  day_raws_df['Date'].apply(lambda x: int(x.strftime('%Y')))
  day_raws_df['Date'] =  day_raws_df['Date'].apply(lambda x: x.strftime('%Y-%m-%d'))
  return day_raws_df


# + id="i9d6sxemMQ21"
dailyRAWS_2012_2020 = getDailyRAWSByFire(combined_csv)

# + id="9M9YDGtumHBE" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601665348234, "user_tz": 420, "elapsed": 397, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="850b2f4a-3f02-42c8-a164-c6d429b4b988"
len(dailyRAWS_2012_2020['Fire'].unique())
dailyRAWS_2012_2020.groupby(['Fire', 'Year']).ngroups

# + id="tOcas2MvOnNS"
df_list = [dailyRAWS_2012_2020, gridmet_2012_2020]

merged_RAWS_2012_2020 = reduce(lambda  left,right: pd.merge(left,right,on=['Fire', 'Year', 'Date'],
                                            how='inner'), df_list)

# + id="fkteNW2UNj8F" colab={"base_uri": "https://localhost:8080/", "height": 34} executionInfo={"status": "ok", "timestamp": 1601665427067, "user_tz": 420, "elapsed": 380, "user": {"displayName": "<NAME>", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gh_Xg3Y--EQFQ6RVDUS0SoJsVKUb6AOaaVkYu-7=s64", "userId": "04192059296941555235"}} outputId="ec3e9aee-ba0a-447c-dcc4-5df8e4ff5d63"
merged_RAWS_2012_2020['FireID'].nunique()

# + id="1m2nlm-ydZSl"
merged_RAWS_2012_2020.to_csv('Variables/Final_Variables/RAWS_fullday_2012_2020.csv')

# + id="jaEum7l7AFl8"


StarcoderdataJupyterScriptsDedupFiltered