Notebook

Visualize all Cascade Volcanoes¶

In [1]:

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
import matplotlib.pyplot as plt
import xarray as xr
from datetime import datetime
import matplotlib
import ipyleaflet
import sys
import os
import dask_gateway
import planetary_computer
sys.path.append('../sar_snowmelt_timing')
import s1_rtc_bs_utils
import contextily as ctx
import rioxarray as rxr
import pathlib
import glob
import re
from matplotlib.pyplot import cm
import scipy
import matplotlib.cm as cm
from matplotlib.colors import Normalize

In [2]:

names = ['MtBaker',
         'GlacierPeak',
         'MtRainier',
         'MtAdams',
         'MtStHelens',
         'MtHood',
         'MtJefferson',
         'ThreeSisters',
         #'CraterLake',
         'MtShasta',
         'LassenPeak',]

geojsons = ['../input/shapefiles/mt_baker.geojson',
            '../input/shapefiles/glacier_peak.geojson',
            '../input/shapefiles/mt_rainier.geojson',
            '../input/shapefiles/mt_adams.geojson',
            '../input/shapefiles/mt_st_helens.geojson',
            '../input/shapefiles/mt_hood.geojson',
            '../input/shapefiles/mt_jefferson.geojson',
            '../input/shapefiles/three_sisters.geojson',
            #'../input/shapefiles/crater_lake.geojson',
            '../input/shapefiles/mt_shasta.geojson',
            '../input/shapefiles/lassen_peak.geojson',]

In [3]:

for nm in names:
    geotiff_list = glob.glob(f'../output/cascades/{nm}/runoff_onset_*_2*.tif')
    
    year_list = []
    for geotiff in geotiff_list:
        year_list.append(re.search("([0-9]{4})", geotiff).group(0))
    year_list = [int(year) for year in year_list]

 #Create variable used for time axis
    time_var = xr.Variable('time', year_list)

# Load in and concatenate all individual GeoTIFFs
    runoffs_allyears = xr.concat([rxr.open_rasterio(i) for i in geotiff_list],dim=time_var).squeeze().sortby('time')
    
    plots = runoffs_allyears.plot(col='time')   
    
    for ax in plots.axes.flatten():
        ax.set_xlabel('')
        ax.set_ylabel('')
        ax.axis('off')

In [4]:

f,ax=plt.subplots(figsize=(10,10))

for nm in names:
    geotiff_list = glob.glob(f'../output/cascades/{nm}/runoff_onset_*_2*.tif')
    
    year_list = []
    for geotiff in geotiff_list:
        year_list.append(re.search("([0-9]{4})", geotiff).group(0))
    year_list = [int(year) for year in year_list]

 #Create variable used for time axis
    time_var = xr.Variable('time', year_list)

# Load in and concatenate all individual GeoTIFFs
    runoffs_allyears = xr.concat([rxr.open_rasterio(i) for i in geotiff_list],dim=time_var).squeeze().sortby('time')
    
    #runoffs_allyears.plot(col='time')
    
    means = runoffs_allyears.mean(dim=['x','y'])
    
    ax.plot(means,label=f'{nm}')
    
    print(f'For {nm}, 2015 snowmelt: {means[0]:.2f} Median snowmelt: {means.median():.2f}')
    print(f'On average 2015 snowmelt occured {means.median() - means[0] :.2f} days before the median snowmelt.')
    
ax.legend()
ax.set_title('Average DOY of Snowmelt Runoff Onset vs Year')
ax.set_label('Average DOY OF Snowmelt Runoff Onset')

For MtBaker, 2015 snowmelt: 108.69 Median snowmelt: 142.01
On average 2015 snowmelt occured 33.33 days before the median snowmelt.
For GlacierPeak, 2015 snowmelt: 110.45 Median snowmelt: 142.49
On average 2015 snowmelt occured 32.04 days before the median snowmelt.
For MtRainier, 2015 snowmelt: 114.68 Median snowmelt: 143.59
On average 2015 snowmelt occured 28.91 days before the median snowmelt.
For MtAdams, 2015 snowmelt: 104.03 Median snowmelt: 140.01
On average 2015 snowmelt occured 35.99 days before the median snowmelt.
For MtStHelens, 2015 snowmelt: 100.04 Median snowmelt: 103.36
On average 2015 snowmelt occured 3.32 days before the median snowmelt.
For MtHood, 2015 snowmelt: 104.39 Median snowmelt: 131.28
On average 2015 snowmelt occured 26.89 days before the median snowmelt.
For MtJefferson, 2015 snowmelt: 105.40 Median snowmelt: 133.43
On average 2015 snowmelt occured 28.03 days before the median snowmelt.
For ThreeSisters, 2015 snowmelt: 103.06 Median snowmelt: 136.32
On average 2015 snowmelt occured 33.26 days before the median snowmelt.
For MtShasta, 2015 snowmelt: 102.15 Median snowmelt: 136.31
On average 2015 snowmelt occured 34.16 days before the median snowmelt.
For LassenPeak, 2015 snowmelt: 100.96 Median snowmelt: 124.19
On average 2015 snowmelt occured 23.23 days before the median snowmelt.

In [5]:

f,ax = plt.subplots(5,2,figsize=(6,10))

cmap_min = 80
cmap_max = 230


cmap = cm.viridis
norm = Normalize(vmin=cmap_min, vmax=cmap_max)

for nm,axs in zip(names,ax.ravel()):
    filename = glob.glob(f'../output/cascades/{nm}/*allyears_median*.tif')
    plots = rxr.open_rasterio(filename[0])
    plots.plot(ax=axs,vmin=cmap_min,vmax=cmap_max,add_colorbar=False)
    axs.set_title(f'{nm}')
    axs.set_aspect('equal')
    axs.set_xlabel('')
    axs.set_ylabel('')
    axs.axis('off')
    ctx.add_basemap(ax=axs, crs=plots.rio.crs, source=ctx.providers.Stamen.Terrain, attribution=False)
f.suptitle(f'2015-2022 Median Snowmelt Runoff Onset Day of Year')
f.tight_layout()


f.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax,label='2015-2022 Median Runoff Onset [DOY]')

Out[5]:

<matplotlib.colorbar.Colorbar at 0x7f6eb5fe0a60>

In [6]:

cmap = cm.viridis
norm = Normalize(vmin=80, vmax=200)

f,ax = plt.subplots(5,2,figsize=(6,10))

for nm,axs in zip(names,ax.ravel()):
    filename_median = glob.glob(f'../output/cascades/{nm}/*allyears_median*.tif')
    median_runoff = rxr.open_rasterio(filename_median[0]).squeeze()
    
    filename_dem = glob.glob(f'../output/cascades/{nm}/*dem*.tif')
    dem = rxr.open_rasterio(filename_dem[0]).squeeze().where(median_runoff>-1)
    
    
    stats_df = pd.DataFrame(columns=['runoff_dates','elevation'])
    stats_df['runoff_dates'] = median_runoff.data.reshape(-1)
    stats_df['elevation'] = dem.data.reshape(-1)
    
    stats_df = stats_df.dropna()
    stats_df.groupby(pd.cut(stats_df["elevation"], 20)).count()['runoff_dates'].plot(ax=axs,kind='barh',width=1,color=cmap(norm(stats_df.groupby(pd.cut(stats_df["elevation"], 20)).median()['runoff_dates'])))
    axs.set_title(f'{nm}')
    axs.set_xticks([]) 
    axs.set_yticks([])
    axs.set_ylabel('Elevation Bins')
    axs.set_xlabel('Count of Pixels in Elevation Bin')
    
f.suptitle(f'Hypsometry Colored by 2015-2022 Median DOY Runoff Onset')
f.tight_layout()
f.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax,label='2015-2022 Median Runoff Onset [DOY]')

Out[6]:

<matplotlib.colorbar.Colorbar at 0x7f6ecdfbde40>

In [7]:

f,ax = plt.subplots(5,4,figsize=(6,10),gridspec_kw={'width_ratios': [3, 1, 3, 1]})

cmap_min = 80
cmap_max = 230


cmap = cm.viridis
norm = Normalize(vmin=cmap_min, vmax=cmap_max)

duplicated_names = []
for name in names:
    duplicated_names.append(name)
    duplicated_names.append(name)
    
for i,(nm,axs) in enumerate(zip(duplicated_names,ax.ravel())):
    if i % 2 == 0:
        filename = glob.glob(f'../output/cascades/{nm}/*allyears_median*.tif')
        plots = rxr.open_rasterio(filename[0])
        plots.plot(ax=axs,vmin=cmap_min,vmax=cmap_max,add_colorbar=False)
        axs.set_title(f'{nm}')
        axs.set_aspect('equal')
        axs.set_xlabel('')
        axs.set_ylabel('')
        axs.axis('off')
        ctx.add_basemap(ax=axs, crs=plots.rio.crs, source=ctx.providers.Stamen.Terrain, attribution=False)
    else: 
        filename_median = glob.glob(f'../output/cascades/{nm}/*allyears_median*.tif')
        median_runoff = rxr.open_rasterio(filename_median[0]).squeeze()
    
        filename_dem = glob.glob(f'../output/cascades/{nm}/*dem*.tif')
        dem = rxr.open_rasterio(filename_dem[0]).squeeze().where(median_runoff>-1)
    
    
        stats_df = pd.DataFrame(columns=['runoff_dates','elevation'])
        stats_df['runoff_dates'] = median_runoff.data.reshape(-1)
        stats_df['elevation'] = dem.data.reshape(-1)
    
        stats_df = stats_df.dropna()
        stats_df.groupby(pd.cut(stats_df["elevation"], 20)).count()['runoff_dates'].plot(ax=axs,kind='barh',width=1,color=cmap(norm(stats_df.groupby(pd.cut(stats_df["elevation"], 20)).median()['runoff_dates'])))
        axs.set_title('')
        axs.set_xticks([]) 
        axs.set_yticks([])
        axs.set_ylabel('')
    
f.suptitle(f'2015-2022 Median Snowmelt Runoff Onset Day of Year with Hypsometry')

f.tight_layout()
f.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax,label='2015-2022 Median Runoff Onset [DOY]')

Out[7]:

<matplotlib.colorbar.Colorbar at 0x7f6ecc184a60>

In [8]:

for nm in names:
    filename_median = glob.glob(f'../output/cascades/{nm}/*allyears_median*.tif')
    median_runoff = rxr.open_rasterio(filename_median[0]).squeeze()
    
    filename_dem = glob.glob(f'../output/cascades/{nm}/*dem*.tif')
    dem = rxr.open_rasterio(filename_dem[0]).squeeze().where(median_runoff>-1)
    
    filename_aspect = glob.glob(f'../output/cascades/{nm}/*aspect*.tif')
    aspect = rxr.open_rasterio(filename_aspect[0]).squeeze().where(median_runoff>-1)
    
    filename_slope = glob.glob(f'../output/cascades/{nm}/*slope*.tif')
    slope = rxr.open_rasterio(filename_slope[0]).squeeze().where(median_runoff>-1)
    
    stats_df = pd.DataFrame(columns=['runoff_dates','elevation','aspect','slope'])
    stats_df['runoff_dates'] = median_runoff.data.reshape(-1)
    stats_df['elevation'] = dem.data.reshape(-1)
    stats_df['aspect'] = np.cos(np.deg2rad(aspect.data.reshape(-1)))
    stats_df['slope'] = slope.data.reshape(-1)
    stats_df['dah'] = np.cos(np.deg2rad(202.5)-np.deg2rad(aspect.data.reshape(-1)))*np.arctan(np.deg2rad(slope.data.reshape(-1)))
    stats_df = stats_df.dropna()
                 
    predictors = np.append(np.ones_like([stats_df['runoff_dates']]).T,stats_df.loc[:,['elevation']].to_numpy(),axis=1)
    B,_,_,_ = scipy.linalg.lstsq(predictors, stats_df.loc[:,'runoff_dates'])
    print(f'For {nm} during 2015-2022, the average runoff onset date is predicted {B[1]*100:.2f} days later for every 100m of elevation gain.')
    print(stats_df.corr())
    print(f'---------------------------------------------------------------------------------------------------------------------------')

For MtBaker during 2015-2022, the average runoff onset date is predicted 6.04 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.886714  0.181437  0.022553 -0.146489
elevation         0.886714   1.000000  0.039475  0.113781  0.000104
aspect            0.181437   0.039475  1.000000  0.000335 -0.860765
slope             0.022553   0.113781  0.000335  1.000000  0.012618
dah              -0.146489   0.000104 -0.860765  0.012618  1.000000
---------------------------------------------------------------------------------------------------------------------------
For GlacierPeak during 2015-2022, the average runoff onset date is predicted 5.80 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.823068  0.226649 -0.109422 -0.222444
elevation         0.823068   1.000000 -0.076671  0.018794  0.059215
aspect            0.226649  -0.076671  1.000000 -0.074629 -0.883128
slope            -0.109422   0.018794 -0.074629  1.000000  0.048337
dah              -0.222444   0.059215 -0.883128  0.048337  1.000000
---------------------------------------------------------------------------------------------------------------------------
For MtRainier during 2015-2022, the average runoff onset date is predicted 4.65 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.875413  0.127707  0.082632 -0.157772
elevation         0.875413   1.000000 -0.021132  0.245073 -0.000167
aspect            0.127707  -0.021132  1.000000 -0.057488 -0.854259
slope             0.082632   0.245073 -0.057488  1.000000  0.051093
dah              -0.157772  -0.000167 -0.854259  0.051093  1.000000
---------------------------------------------------------------------------------------------------------------------------
For MtAdams during 2015-2022, the average runoff onset date is predicted 5.18 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.842649  0.198514  0.272018 -0.205100
elevation         0.842649   1.000000  0.058961  0.430279 -0.060014
aspect            0.198514   0.058961  1.000000 -0.056586 -0.830179
slope             0.272018   0.430279 -0.056586  1.000000  0.015045
dah              -0.205100  -0.060014 -0.830179  0.015045  1.000000
---------------------------------------------------------------------------------------------------------------------------
For MtStHelens during 2015-2022, the average runoff onset date is predicted 4.27 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.813078  0.013158  0.309737 -0.052175
elevation         0.813078   1.000000 -0.042051  0.455983  0.033239
aspect            0.013158  -0.042051  1.000000 -0.087717 -0.841057
slope             0.309737   0.455983 -0.087717  1.000000  0.018238
dah              -0.052175   0.033239 -0.841057  0.018238  1.000000
---------------------------------------------------------------------------------------------------------------------------
For MtHood during 2015-2022, the average runoff onset date is predicted 5.91 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.801229  0.101277  0.088885 -0.127618
elevation         0.801229   1.000000 -0.085252  0.204024  0.028095
aspect            0.101277  -0.085252  1.000000  0.007668 -0.876505
slope             0.088885   0.204024  0.007668  1.000000 -0.009968
dah              -0.127618   0.028095 -0.876505 -0.009968  1.000000
---------------------------------------------------------------------------------------------------------------------------
For MtJefferson during 2015-2022, the average runoff onset date is predicted 5.80 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.707183  0.238503 -0.021747 -0.230197
elevation         0.707183   1.000000  0.017097  0.183973 -0.024328
aspect            0.238503   0.017097  1.000000 -0.040026 -0.869908
slope            -0.021747   0.183973 -0.040026  1.000000  0.149889
dah              -0.230197  -0.024328 -0.869908  0.149889  1.000000
---------------------------------------------------------------------------------------------------------------------------
For ThreeSisters during 2015-2022, the average runoff onset date is predicted 4.27 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.601043  0.207990  0.108313 -0.219256
elevation         0.601043   1.000000 -0.059353  0.468607  0.041223
aspect            0.207990  -0.059353  1.000000 -0.045623 -0.822953
slope             0.108313   0.468607 -0.045623  1.000000 -0.005446
dah              -0.219256   0.041223 -0.822953 -0.005446  1.000000
---------------------------------------------------------------------------------------------------------------------------
For MtShasta during 2015-2022, the average runoff onset date is predicted 4.32 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.803060  0.307896  0.194023 -0.294400
elevation         0.803060   1.000000  0.074466  0.381479 -0.007221
aspect            0.307896   0.074466  1.000000 -0.044537 -0.884200
slope             0.194023   0.381479 -0.044537  1.000000  0.067924
dah              -0.294400  -0.007221 -0.884200  0.067924  1.000000
---------------------------------------------------------------------------------------------------------------------------
For LassenPeak during 2015-2022, the average runoff onset date is predicted 2.53 days later for every 100m of elevation gain.
              runoff_dates  elevation    aspect     slope       dah
runoff_dates      1.000000   0.416987  0.350259  0.122971 -0.449301
elevation         0.416987   1.000000 -0.208923  0.479766  0.164401
aspect            0.350259  -0.208923  1.000000  0.000154 -0.855658
slope             0.122971   0.479766  0.000154  1.000000 -0.029238
dah              -0.449301   0.164401 -0.855658 -0.029238  1.000000
---------------------------------------------------------------------------------------------------------------------------

In [9]:

f,ax = plt.subplots(5,2,figsize=(6,10))

for nm,axs in zip(names,ax.ravel()):
    filename_median = glob.glob(f'../output/cascades/{nm}/*allyears_median*.tif')
    median_runoff = rxr.open_rasterio(filename_median[0]).squeeze()
    
    filename_dem = glob.glob(f'../output/cascades/{nm}/*dem*.tif')
    dem = rxr.open_rasterio(filename_dem[0]).squeeze().where(median_runoff>-1)
    
    filename_aspect = glob.glob(f'../output/cascades/{nm}/*aspect*.tif')
    aspect = rxr.open_rasterio(filename_aspect[0]).squeeze().where(median_runoff>-1)
    
    filename_slope = glob.glob(f'../output/cascades/{nm}/*slope*.tif')
    slope = rxr.open_rasterio(filename_slope[0]).squeeze().where(median_runoff>-1)
    
    stats_df = pd.DataFrame(columns=['runoff_dates','elevation','aspect','slope'])
    stats_df['runoff_dates'] = median_runoff.data.reshape(-1)
    stats_df['elevation'] = dem.data.reshape(-1)
    stats_df['aspect'] = np.cos(np.deg2rad(aspect.data.reshape(-1)))
    stats_df['slope'] = slope.data.reshape(-1)
    stats_df['dah'] = np.cos(np.deg2rad(202.5)-np.deg2rad(aspect.data.reshape(-1)))*np.arctan(np.deg2rad(slope.data.reshape(-1)))
    stats_df = stats_df.dropna()
    
    stats_elevation = stats_df.groupby(pd.cut(stats_df["elevation"], 60)).corr()
    #stats_elevation['runoff_dates'].unstack().drop(columns='elevation').drop(columns=['runoff_dates','slope','aspect']).plot(kind='barh',ax=axs,width=1) #.reset_index()
    stats_elevation['runoff_dates'].unstack().drop(columns='elevation').drop(columns=['runoff_dates','slope','dah']).plot(kind='barh',ax=axs,width=1) #.reset_index()
    #stats_elevation['runoff_dates'].unstack().drop(columns='elevation').drop(columns=['runoff_dates','dah','aspect']).plot(kind='barh',ax=axs,width=1) #.reset_index()

    #stats_df.groupby(pd.cut(stats_df["elevation"], 20)).count()['runoff_dates'].plot(ax=axs,kind='barh',width=1,color=cmap(norm(stats_df.groupby(pd.cut(stats_df["elevation"], 20)).median()['runoff_dates'])))
    axs.set_title(f'{nm}')
    #axs.set_xticks([]) 
    axs.set_yticks([])
    axs.set_ylabel('Elevation Bins')
    axs.set_xlabel('Correlation Strength')
    axs.set_xlim([0,0.7])
    
#f.suptitle(f'DAH Correlation Strength at Each Elevation Bin')
f.suptitle(f'Snowmelt Runoff Onset Timing and Aspect Correlation Strength at Each Elevation Bin')

f.tight_layout()

In [10]:

f,ax=plt.subplots(figsize=(8,12))
for nm in names:
    print(f'Processing {nm}...')
    filename_median = glob.glob(f'../output/cascades/{nm}/*allyears_median*.tif')
    median_runoff = rxr.open_rasterio(filename_median[0]).squeeze()
    filename_dem = glob.glob(f'../output/cascades/{nm}/*dem*.tif')
    dem = rxr.open_rasterio(filename_dem[0]).squeeze().where(median_runoff>-1)
    dem_and_runoff_df = pd.DataFrame((dem.to_numpy().ravel(),median_runoff.to_numpy().ravel())).T.dropna()
    dem_and_runoff_df.columns = ['Elevation','Runoff DOY']

    median = dem_and_runoff_df.groupby(pd.qcut(dem_and_runoff_df["Elevation"], 20)).median()
    std = dem_and_runoff_df.groupby(pd.qcut(dem_and_runoff_df["Elevation"], 20)).std()
    median.plot(ax=ax,x='Runoff DOY',y='Elevation',label=f'{nm}')
    ax.fill_betweenx(median['Elevation'],median['Runoff DOY']-std['Runoff DOY'],median['Runoff DOY']+std['Runoff DOY'],alpha=0.1)
    

ax.set_title(f'2015-2022 Median Snowmelt Runoff Onset vs Elevation')
ax.set_xlabel(f'Median Snowmelt Runoff Onset [DOY]')
ax.set_ylabel(f'Elevation [m]')    
ax.grid()
ax.legend(loc='lower right')

Processing MtBaker...
Processing GlacierPeak...
Processing MtRainier...
Processing MtAdams...
Processing MtStHelens...
Processing MtHood...
Processing MtJefferson...
Processing ThreeSisters...
Processing MtShasta...
Processing LassenPeak...

Out[10]:

<matplotlib.legend.Legend at 0x7f6eb5fabaf0>

In [11]:

f,ax = plt.subplots(5,2,figsize=(12,20))

for nm,axs in zip(names,ax.ravel()):
    filename = glob.glob(f'../output/cascades/{nm}/*allyears_std*.tif')
    plots = rxr.open_rasterio(filename[0])
    plots.plot(ax=axs,cmap='Reds',vmin=0,vmax=60,add_colorbar=False)
    axs.set_title(f'{nm}')
    axs.set_aspect('equal')
    axs.set_xlabel('')
    axs.set_ylabel('')
    axs.axis('off')
    ctx.add_basemap(ax=axs, crs=plots.rio.crs, source=ctx.providers.Stamen.Terrain, attribution=False)
f.suptitle(f'2015-2022 Snowmelt Runoff Onset Day of Year Standard Devation')
f.tight_layout()

In [12]:

f,ax = plt.subplots(5,2,figsize=(12,20))

for nm,axs in zip(names,ax.ravel()):
    filename = glob.glob(f'../output/cascades/{nm}/*allyears_polyfit*.tif')
    plots = rxr.open_rasterio(filename[0])
    plots.plot(ax=axs,cmap='RdBu',vmin=-10,vmax=10,add_colorbar=False)
    axs.set_title(f'{nm}')
    axs.set_aspect('equal')
    axs.set_xlabel('')
    axs.set_ylabel('')
    axs.axis('off')
    ctx.add_basemap(ax=axs, crs=plots.rio.crs, source=ctx.providers.Stamen.Terrain, attribution=False)
f.suptitle(f'2015-2022 Snowmelt Runoff Onset Day of Year Trend')
f.tight_layout()

In [13]:

f,ax = plt.subplots(5,2,figsize=(10,16))

for nm,axs in zip(names,ax.ravel()):
    print(f'Processing {nm}...')
    geotiff_list = glob.glob(f'../output/cascades/{nm}/runoff_onset_*_2*.tif')
    year_list = []
    for geotiff in geotiff_list:
        year_list.append(re.search("([0-9]{4})", geotiff).group(0))
    year_list = [int(year) for year in year_list]

    #Create variable used for time axis
    time_var = xr.Variable('time', year_list)

    # Load in and concatenate all individual GeoTIFFs
    runoffs_allyears = xr.concat([rxr.open_rasterio(i) for i in geotiff_list],dim=time_var).squeeze().sortby('time')
    
    filename_dem = glob.glob(f'../output/cascades/{nm}/*dem*.tif')
    dem = rxr.open_rasterio(filename_dem[0]).squeeze().where(runoffs_allyears.isel(time=0)>-1)
    
    color = iter(cm.Blues(np.linspace(0.2, 0.99, 8)))

    year_list.sort()
    for year in year_list:
        median_runoff = runoffs_allyears[runoffs_allyears.time==year]
        dem_and_runoff_df = pd.DataFrame((dem.data.reshape(-1),median_runoff.data.reshape(-1))).T.dropna()
        dem_and_runoff_df.columns = ['Elevation','Runoff DOY']
        c = next(color)
        dem_and_runoff_df.groupby(pd.qcut(dem_and_runoff_df["Elevation"], 10)).median().plot(ax=axs,x='Runoff DOY',y='Elevation',label=f'{year}',c=c)
        axs.set_title(f'{nm}')
f.suptitle(f'Elevation vs. Runoff Date')
f.tight_layout()

Processing MtBaker...
Processing GlacierPeak...
Processing MtRainier...
Processing MtAdams...
Processing MtStHelens...
Processing MtHood...
Processing MtJefferson...
Processing ThreeSisters...
Processing MtShasta...
Processing LassenPeak...

In [ ]: