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

# # Tracking the Smoke Caused by the fires
#
# In this example we show how to use [HRRR Smoke Experimental dataset](https://data.planetos.com/datasets/noaa_hrrr_wrf_smoke) to analyse smoke in the US and we will also download historical fire data from [Cal Fire](https://www.fire.ca.gov/incidents) web page to visualize burned area since 2013.
#
# [The High-Resolution Rapid Refresh Smoke (HRRRSmoke)](https://data.planetos.com/datasets/noaa_hrrr_wrf_smoke) is a three-dimensional model that allows simulation of mesoscale flows and smoke dispersion over complex terrain, in the boundary layer and aloft at high spatial resolution over the CONUS domain. The smoke model comprises a suite of fire and environmental products for forecasters during the fire weather season. Products derived from the HRRRSmoke model include the Fire Radiative Power (FRP), Near-Surface Smoke (PM2.5), and Vertically Integrated Smoke, to complement the 10-meter winds, 1-hour precipitation, 2-meter temperature and surface visibility experimental forecast products. Keep in mind, that this dataset is EXPERIMENTAL. Therefore, they should not be used to make decisions regarding safety of life or property.
#
# HRRR Smoke has many different weather parameters, two of them are directly smoke related - Column-integrated mass density and Mass density (concentration) @ Specified height level above ground. First is vertically-integrated smoke and second smoke on the lowest model level (8 m). 
#
# As we would like to have data about the Continental United States we will download data by using Package API. Then we will create a widget where you can choose timestamp by using a slider. After that, we will also save the same data as a GIF to make sharing the results with friends and colleagues more fun. And finally, we will compare smoke data with CAMS particulate matter data to find out if there's a colleration between them as HRRR Smoke is still experimental.

# %matplotlib notebook
# %matplotlib inline
import numpy as np
import dh_py_access.lib.datahub as datahub
import xarray as xr
import matplotlib.pyplot as plt
import ipywidgets as widgets
from mpl_toolkits.basemap import Basemap,shiftgrid
import dh_py_access.package_api as package_api
import matplotlib.colors as colors
import warnings
import datetime
import shutil
import imageio
import seaborn as sns
import pandas as pd
import os
import matplotlib as mpl
import wget
warnings.filterwarnings("ignore")

# <font color='red'>Please put your datahub API key into a file called APIKEY and place it to the notebook folder or assign your API key directly to the variable API_key!</font>

server = 'api.planetos.com'
API_key = open('APIKEY').readlines()[0].strip() #'<YOUR API KEY HERE>'
version = 'v1'

# At first, we need to define the dataset name and a variable we want to use.

dh = datahub.datahub(server,version,API_key)
dataset = 'noaa_hrrr_wrf_smoke'
variable_name1 = 'Mass_density_concentration_height_above_ground'

# Then we define spatial range. We decided to analyze US, where unfortunately catastrofic wildfires are taking place at the moment and influeces air quality. 

# +
# reftime = datetime.datetime.strftime(datetime.datetime.today(), '%Y-%m-%d') + 'T00:00:00'
area_name = 'usa'
today_hr = datetime.datetime.strftime(datetime.datetime.today(),'%Y%m%dT%H')

latitude_north = 49; longitude_west = -127
latitude_south = 26; longitude_east = -70.5
# -



# ### Download the data with package API
#
# 1. Create package objects
# 2. Send commands for the package creation
# 3. Download the package files

package_hrrr = package_api.package_api(dh,dataset,variable_name1,longitude_west,longitude_east,latitude_south,latitude_north,area_name=area_name+today_hr)

package_hrrr.make_package()

package_hrrr.download_package()

# ### Work with the downloaded files
#
# We start with opening the files with xarray. After that, we will create a map plot with a time slider, then make a GIF using the images, then we will do the same thing for closer area - California; and finally, we will download csv file about fires in California to visualize yearly incidents data as a bar chart.

dd1 = xr.open_dataset(package_hrrr.local_file_name)
dd1['longitude'] = ((dd1.lon+180) % 360) - 180
dd1[variable_name1].data[dd1[variable_name1].data < 0] = 0
dd1[variable_name1].data[dd1[variable_name1].data == np.nan] = 0



# Here we are making a Basemap of the US that we will use for showing the data.

m = Basemap(projection='merc', lat_0 = 55, lon_0 = -4,
         resolution = 'h', area_thresh = 0.05,
         llcrnrlon=longitude_west, llcrnrlat=latitude_south,
         urcrnrlon=longitude_east, urcrnrlat=latitude_north)
lons,lats = np.meshgrid(dd1.longitude.data,dd1.lat.data)
lonmap,latmap = m(lons,lats)

# Now it is time to plot all the data. A great way to do it is to make an interactive widget, where you can choose time stamp by using a slider. 
#
# As the minimum and maximum values are very different, we are using logarithmic colorbar to visualize it better.
#
# On the map we can see that the areas near fires have more smoke, but it travels pretty far. Depending on when the notebook is run, we can see very different results.
#
# But first we define minimum, maximum and also colormap.

# +
vmax = np.nanmax(dd1[variable_name1].data)
vmin = 2

cmap = mpl.cm.twilight.colors[:-100]
tmap = mpl.colors.LinearSegmentedColormap.from_list('twilight_edited', cmap)

# -



def loadimg(k):
    fig=plt.figure(figsize=(10,7))
    ax = fig.add_subplot(111)
    pcm = m.pcolormesh(lonmap,latmap,dd1[variable_name1].data[k][0],
                norm = colors.LogNorm(vmin=vmin, vmax=vmax),cmap = tmap)
    ilat,ilon = np.unravel_index(np.nanargmax(dd1[variable_name1].data[k][0]),dd1[variable_name1].data[k][0].shape)
    cbar = plt.colorbar(pcm,fraction=0.024, pad=0.040,ticks=[10**0, 10**1, 10**2,10**3])
    cbar.ax.set_yticklabels([0,10,100,1000]) 
    ttl = plt.title('Near Surface Smoke ' + str(dd1[variable_name1].time[k].data)[:-10],fontsize=20,fontweight = 'bold')
    ttl.set_position([.5, 1.05])
    cbar.set_label(dd1[variable_name1].units)

    m.drawcountries()
    m.drawstates()
    m.drawcoastlines()
    print("Maximum: ","%.2f" % np.nanmax(dd1[variable_name1].data[k][0]))

    plt.show()
widgets.interact(loadimg, k=widgets.IntSlider(min=0,max=len(dd1[variable_name1].data)-1,step=1,value=0, layout=widgets.Layout(width='100%')))

# Let's include an image from the last time-step as well, because GitHub Preview doesn't show the time slider images. 

loadimg(9)


# With the function below we will save images you saw above to the local filesystem as a GIF, so it is easily shareable with others. 

def make_ani(m,lonmap,latmap,aniname,smaller_area=False):
    if smaller_area==True:
        fraction = 0.035
        fontsize = 13
    else:
        fraction = 0.024
        fontsize = 20
        
    folder = './anim/'
    for k in range(len(dd1[variable_name1])):
        filename = folder + 'ani_' + str(k).rjust(3,'0') + '.png'
        if not os.path.exists(filename):
            fig=plt.figure(figsize=(10,7))
            ax = fig.add_subplot(111)
            pcm = m.pcolormesh(lonmap,latmap,dd1[variable_name1].data[k][0],
            norm = colors.LogNorm(vmin=vmin, vmax=vmax),cmap = tmap)

            m.drawcoastlines()
            m.drawcountries()
            m.drawstates()
            cbar = plt.colorbar(pcm,fraction=fraction, pad=0.040,ticks=[10**0, 10**1, 10**2,10**3])
            cbar.ax.set_yticklabels([0,10,100,1000]) 
            ttl = plt.title('Near Surface Smoke ' + str(dd1[variable_name1].time[k].data)[:-10],fontsize=fontsize,fontweight = 'bold')
            ttl.set_position([.5, 1.05])
            cbar.set_label(dd1[variable_name1].units)
            ax.set_xlim()
            if not os.path.exists(folder):
                os.mkdir(folder)
            plt.savefig(filename,bbox_inches = 'tight',dpi=150)
            plt.close()

    files = sorted(os.listdir(folder))
    images = []
    for file in files:
        if not file.startswith('.'):
            filename = folder + file
            images.append(imageio.imread(filename))
    kargs = { 'duration': 0.3,'quantizer':2,'fps':5.0}
    imageio.mimsave(aniname, images, **kargs)
    print ('GIF is saved as {0} under current working directory'.format(aniname))
    shutil.rmtree(folder)
make_ani(m,lonmap,latmap,'hrrr_smoke.gif')

# As we are interested in California fires right now, it would make sense to make animation of only California area as well. So people can be prepared when smoke hits their area. The model has pretty good spatial resolution as well - 3 km, which makes tracking the smoke easier. 

# +
latitude_north_cal = 43; longitude_west_cal = -126.
latitude_south_cal = 30.5; longitude_east_cal = -113

m2 = Basemap(projection='merc', lat_0 = 55, lon_0 = -4,
         resolution = 'h', area_thresh = 0.05,
         llcrnrlon=longitude_west_cal, llcrnrlat=latitude_south_cal,
         urcrnrlon=longitude_east_cal, urcrnrlat=latitude_north_cal)
lons2,lats2 = np.meshgrid(dd1.longitude.data,dd1.lat.data)
lonmap_cal,latmap_cal = m2(lons2,lats2)
# -

make_ani(m2,lonmap_cal,latmap_cal,'hrrr_smoke_california.gif',smaller_area=True)

# Finally, we will remove the package we downloaded. 

os.remove(package_hrrr.local_file_name)

# # Data about Burned Area from Cal Fire
# Now we will download csv file from Cal Fire web page and illustrate how many acres each year was burnt since 2013. 

if not os.path.exists('acres_burned.csv'):
    wget.download('https://www.fire.ca.gov/imapdata/mapdataall.csv',out='acres_burned.csv')

datain = pd.read_csv('acres_burned.csv')

# Here we convert `incident_dateonly_created` column to datetime, so it's easier to group data by year.

datain['incident_dateonly_created'] = pd.to_datetime(datain['incident_dateonly_created'])

# Below you can see the data from `acres_burned.csv` file. It has information about each incident. This time we only compute total acres burned each year. 

datain

# Computing yearly sums. In some reason there's many years without much data, so we will filter it out. Also, reseting index, as we don't want dates to be as an index and making year column. 

burned_acres_yearly = datain.resample('1AS', on='incident_dateonly_created')['incident_acres_burned'].sum()
burned_acres_yearly = burned_acres_yearly[burned_acres_yearly.index > datetime.datetime(2012,1,1)]
burned_acres_yearly = burned_acres_yearly.reset_index()
burned_acres_yearly['year'] = pd.DatetimeIndex(burned_acres_yearly.incident_dateonly_created).year

# We can see the computed data below.

burned_acres_yearly

# Finally we will make a bar chart of the data. We are using [seaborn](https://seaborn.pydata.org/) this time for plotting the data and to visualize it better, we added colormap to bar chart as well. 
# Image will be saved into the working directory. 

# +
fig,ax = plt.subplots(figsize=(10,6))
pal = sns.color_palette("YlOrRd_r", len(burned_acres_yearly))
rank = burned_acres_yearly['incident_acres_burned'].argsort().argsort()

sns.barplot(x='year',y='incident_acres_burned',data=burned_acres_yearly,ci=95,ax=ax,palette=np.array(pal[::-1])[rank])
ax.set_xlabel('Year',fontsize=15)

ax.set_ylabel('Burned Area [acres]',fontsize=15)
ax.grid(color='#C3C8CE',alpha=1)
ax.set_axisbelow(True)
ax.spines['bottom'].set_color('#C3C8CE')
ax.spines['top'].set_color('#C3C8CE')
ax.spines['left'].set_color('#C3C8CE')
ax.spines['right'].set_color('#C3C8CE')
ttl = ax.set_title('Burned Area in California',fontsize=20,fontweight = 'bold')
ttl.set_position([.5, 1.05])
ax.tick_params(labelsize=15,length=0)
plt.savefig('acres_burned_cali.png',dpi=300)

StarcoderdataJupyterScriptsDedupFiltered