#!/usr/bin/env python # coding: utf-8 # # Snowmelt timing with Sentinel-1 SAR at Grand Mesa, CO # #### Applying my [SAR snowmelt timing toolbox](https://github.com/egagli/sar_snow_melt_timing) to Grand Mesa, CO for snowex hackweek. Snowex hackweek repo can be found [here](https://github.com/snowex-hackweek/snowmelt-timing) and presentation can be found [here](https://docs.google.com/presentation/d/1D-t0BT5zrmGV-8o6O4DEFhmdDeUwdOxDfkmsYYA0vOg/edit#slide=id.g13d0e1ed0eb_4_14). # #### Our timeframe of interest is December 1st 2019 - August 1st 2020. # In[1]: import numpy as np import pandas as pd import geopandas as gpd from dask.distributed import Client import py3dep import geopandas as gpd import rasterio as rio import pystac import pystac_client import stackstac import math import shapely import matplotlib.pyplot as plt import os import xarray as xr from datetime import datetime import ulmo from datetime import datetime import sys import contextily as ctx import rioxarray as rxr from matplotlib_scalebar.scalebar import ScaleBar import matplotlib import matplotlib.transforms as mtransforms import logging import geopandas as gpd import requests sys.path.append('../sar_snowmelt_timing') import s1_rtc_bs_utils # In[2]: #### CODE FOR DEVELOPMENT-- THIS RELOADS s1_rtc_bs_utils MODULE sys.modules.pop('s1_rtc_bs_utils') import s1_rtc_bs_utils # In[3]: # GDAL environment variables for better performance os.environ['AWS_REGION']='us-west-2' os.environ['GDAL_DISABLE_READDIR_ON_OPEN']='EMPTY_DIR' os.environ['AWS_NO_SIGN_REQUEST']='YES' # xr.set_options(keep_attrs=True) # Paste /proxy/localhost:8787 for cluster diagnostics client = Client(local_directory='/tmp', processes=False,silence_logs=logging.ERROR) client # ## Let's set the location via geojson obtained from geojson.io. # In[4]: grandmesa_bbox_gdf = gpd.read_file('../input/shapefiles/grandmesa.geojson') # In[5]: ts_ds = s1_rtc_bs_utils.get_s1_rtc_stac(grandmesa_bbox_gdf,start_time='2015-01-01',end_time=datetime.today().strftime('%Y-%m-%d'),orbit_direction='all',polarization='gamma0_vv',collection='../input/sentinel1-rtc-aws/mycollection_grandmesa.json') ts_ds = ts_ds.compute() ts_ds = ts_ds.dropna('time',how='all') # ## Let's check out the xarray dataset. # In[6]: ts_ds # ## Let's see the timings of the Sentinel-1 acquisitions. # In[7]: f,ax=plt.subplots(figsize=(30,7)) s1_rtc_bs_utils.plot_sentinel1_acquisitons(ts_ds,ax=ax,textsize=6) f.suptitle(f'Grand Mesa, CO') plt.tight_layout() # In[8]: f,ax=plt.subplots(figsize=(30,7)) s1_rtc_bs_utils.plot_sentinel1_acquisitons(ts_ds,ax=ax,start_date='2019-12-01',end_date='2020-08-01',textsize=14) f.suptitle(f'Grand Mesa, CO') plt.tight_layout() # ### Looks like 6pm as expected, sweet. # ## Let's clip to the melt season. # In[9]: ts_ds = s1_rtc_bs_utils.get_s1_rtc_stac(grandmesa_bbox_gdf,start_time='2019-12-01',end_time='2020-08-01',orbit_direction='all',polarization='gamma0_vv',collection='../input/sentinel1-rtc-aws/mycollection_grandmesa.json') ts_ds = ts_ds.compute() ts_ds = ts_ds.dropna('time',how='all') # ## Let's take a preliminary look at RGB and NDSI data from Sentinel-2 throughout the melt season. # In[10]: rgb = s1_rtc_bs_utils.get_s2_rgb(ts_ds) ndsi = s1_rtc_bs_utils.get_s2_ndsi(ts_ds) #ndwi = s1_rtc_bs_utils.get_s2_ndwi(ts_ds) new_rgb = np.clip(rgb/np.amax(rgb),0,1) # In[11]: f,ax=plt.subplots(min(len(ndsi.time),len(rgb.time)),2,figsize=(10,40)) for i in range(min(len(ndsi.time),len(rgb.time))): new_rgb.isel(time=i).plot.imshow(ax=ax[i,0],rgb='band',robust=True) ndsi.isel(time=i).plot.imshow(ax=ax[i,1],add_colorbar=False,cmap='RdBu',vmin=-1,vmax=1) #ndwi.isel(time=i).plot.imshow(ax=ax[i,2],add_colorbar=False,vmin=-1,vmax=1) ax[i,0].set_aspect('equal') ax[i,1].set_aspect('equal') #ax[i,2].set_aspect('equal') ax[i,0].axis('off') ax[i,1].axis('off') #ax[i,2].axis('off') ax[i,0].set_title(f'{pd.to_datetime(rgb.time.values[i]).strftime("%Y-%m-%d")} S2 RGB') ax[i,1].set_title(f'{pd.to_datetime(rgb.time.values[i]).strftime("%Y-%m-%d")} S2 NDSI') #ax[i,2].set_title(f'{pd.to_datetime(rgb.time.values[i]).strftime("%Y-%m-%d")} S2 NDWI') plt.tight_layout() # ## Let's take a peak at the SAR scenes. # In[12]: ts_ds_db = 10.0*np.log10(ts_ds) ts_ds_db.plot(col='time',col_wrap=5,robust=False,cmap='gray',vmin=-15,vmax=-6) # ## Let's check the SWE measurements at nearby SNOTEL sites. # In[13]: f,ax=plt.subplots(1,2,figsize=(20,6),gridspec_kw={'width_ratios': [1, 2]}) s1_rtc_bs_utils.plot_closest_snotel(ts_ds,distance_cutoff=30,ax=ax[0]) snotel_data = s1_rtc_bs_utils.get_closest_snotel_data(ts_ds,variable_code='SNOTEL:WTEQ_D',distance_cutoff=30,start_date='1900-01-01', end_date=datetime.today().strftime('%Y-%m-%d')) snotel_data_cm = snotel_data*2.54 snotel_data_cm.plot(ax=ax[1]) ax[1].set_xlabel('Time') ax[1].set_ylabel('Snow Water Equivalent [cm]') ax[1].set_title('SWE for SNOTEL Stations Close in Proximety to Dataset') ax[1].set_xlim(['2019-12-01','2020-08-01']) plt.tight_layout() # ## Let's try out a plot that combines S1 backscatter, S2 NDSI, and the nearest SNOTEL station. # In[14]: s1_rtc_bs_utils.plot_bs_ndsi_swe_precip_with_context(ts_ds,start_date='2019-12-01',end_date='2020-08-01') # ## Looks like we need to do some masking. Let's first try masking by NDVI. # In[15]: summer_ndvi_ds = s1_rtc_bs_utils.get_median_ndvi(ts_ds) f,ax=plt.subplots() summer_ndvi_ds.plot(ax=ax) # ## Masking by NDVI might be a bit too noisy here. Maybe let's just mask by an elevation height. Let's test a couple thresholds: # In[16]: dem = s1_rtc_bs_utils.get_py3dep_dem(ts_ds) f,ax=plt.subplots(1,5,figsize=(30,6)) dem.plot(ax=ax[0]) dem.where(dem>2800).plot(ax=ax[1]) dem.where(dem>3000).plot(ax=ax[2]) dem.where(dem>3100).plot(ax=ax[3]) dem.where(dem>3200).plot(ax=ax[4]) ax[0].set_title('No Elevation Cutoff') ax[1].set_title('2800m Elevation Cutoff') ax[2].set_title('3000m Elevation Cutoff') ax[3].set_title('3100m Elevation Cutoff') ax[4].set_title('3200m Elevation Cutoff') for axs in ax: axs.set_aspect('equal') # ## 3000m looks like a good cutoff. Let's rerun our analysis with this elevation mask. # In[17]: ts_ds_dem_masked = ts_ds.where(dem>3000) s1_rtc_bs_utils.plot_bs_ndsi_swe_precip_with_context(ts_ds_dem_masked,start_date='2019-12-01',end_date='2020-08-01') # ## Great! Now let's pull in the snowpit location and swe data. (Thank you Megan Mason!) # In[18]: url = "https://gist.githubusercontent.com/meganmason/dde13c46a90875e364b1c25a31bff1d8/raw/d361e1a56eeb867bd89c87b4dc13ed431ffff6f6/SNEX_TS_SP_avg_prelim_Hackweek.geojson" snowpits = gpd.read_file(url) # In[19]: # keep columns we want in the dataframe url = "https://gist.githubusercontent.com/meganmason/40a87840b52c5516efc716236a67ccf3/raw/9767072779ac6f22d13e5e4dbf2c7e85f3e6e15e/SNEX20_TS_SP_cogm_swe.geojson" swe = gpd.read_file(url) trunc = lambda x: x[:6] swe.PitID = swe.PitID.apply(trunc) swe = swe[['PitID','Site', 'Date/Local Standard Time', 'SWE (mm)']] swe = swe.rename({'Date/Local Standard Time': 'date', 'SWE (mm)': 'swe (mm)'}, axis=1) swe.date = pd.to_datetime(swe.date) # ## Let's select all of the Grand Mesa sites and reproject to UTM. # In[20]: snowpits = snowpits[snowpits['Location']=='Grand Mesa'] snowpits = snowpits.to_crs(ts_ds.rio.crs) # In[21]: snowpits # ## Now let's plot our snowpit sites on top of our estimated DOY runoff onset Map. # In[22]: f,ax=plt.subplots(figsize=(12,7)) runoff_map = s1_rtc_bs_utils.get_runoff_onset(ts_ds_dem_masked).dt.dayofyear runoff_map.plot(ax=ax,cmap='twilight') snowpits.plot(ax=ax,categorical=True,column='Site',legend=True,legend_kwds={'bbox_to_anchor':(1.0, .35),'fontsize':14,'frameon':False}) ctx.add_basemap(ax=ax, crs=ts_ds.rio.crs, source=ctx.providers.Stamen.Terrain) ax.set_title('DOY Estimates for Snowmelt Runoff Onset',fontsize=22) ax.set_xlim(left=218000) ax.set_aspect('equal') ax.add_artist(ScaleBar(1.0)) plt.tight_layout() plt.savefig('../output/grandmesa/plots/runoff_map.png') # ## Let's zoom in a bit more. # In[23]: f,ax=plt.subplots(figsize=(20,4)) runoff_map.plot(ax=ax,cmap='twilight') snowpits.plot(ax=ax,categorical=True,column='Site',legend=True,legend_kwds={'bbox_to_anchor':(0.8, -0.1),'fontsize':10,'frameon':False,'ncol':6}) ctx.add_basemap(ax=ax, crs=ts_ds.rio.crs, source=ctx.providers.Stamen.Terrain) ax.set_title('DOY Estimates for Snowmelt Runoff Onset',fontsize=22) ax.set_xlim([220000,240000]) ax.set_ylim([4.324E6,4.327E6]) ax.add_artist(ScaleBar(1.0)) plt.tight_layout() plt.savefig('../output/grandmesa/plots/runoff_map_zoomed.png') # ## Let's check out the runoff estimates around the snowpits. Let's use a 100m buffer. # In[24]: snowpits_buffer = snowpits.copy() snowpits_buffer['geometry'] = snowpits_buffer.geometry.buffer(100) # In[25]: PitIDs = ['COGMCO','COGMCT','COGMWO','COGMWT','COGMSO','COGMST'] f,ax=plt.subplots(3,2,figsize=(8,10)) for pitid, axs in zip(PitIDs,ax.ravel()): snowpit = snowpits_buffer[snowpits_buffer['PitID']==pitid] runoff_map_clipped = runoff_map.rio.clip(snowpit.geometry) runoff_map_clipped.plot(ax=axs,cmap='twilight',vmin=0,vmax=365) #ctx.add_basemap(ax=axs, crs=runoff_map_clipped.rio.crs, source=ctx.providers.Stamen.Terrain) snowpits[snowpits['PitID']==pitid].plot(ax=axs,marker='*',color='green') axs.set_title(f'{snowpit.Site.values[0]} \n Mean Snowmelt Runoff Onset DOY = {runoff_map_clipped.mean().values:.2f}') axs.set_aspect('equal') axs.axis('off') axs.add_artist(ScaleBar(1.0)) plt.tight_layout() plt.savefig('../output/grandmesa/plots/runoff_map_snowpits.png') # In[26]: PitIDs = ['COGMCO','COGMCT','COGMWO','COGMWT','COGMSO','COGMST'] f,ax=plt.subplots(3,2,figsize=(8,10)) for pitid, axs in zip(PitIDs,ax.ravel()): snowpit = snowpits_buffer[snowpits_buffer['PitID']==pitid] runoff_map_clipped = new_rgb.isel(time=5).rio.clip(snowpit.geometry) runoff_map_clipped.plot.imshow(ax=axs,rgb='band',vmax=0.7) #ctx.add_basemap(ax=axs, crs=runoff_map_clipped.rio.crs, source=ctx.providers.Stamen.Terrain) snowpits[snowpits['PitID']==pitid].plot(ax=axs,marker='*',color='green') axs.set_title(f'{snowpit.Site.values[0]}') axs.set_aspect('equal') axs.axis('off') axs.add_artist(ScaleBar(1.0)) plt.tight_layout() plt.savefig('../output/grandmesa/plots/runoff_map_snowpits_ndsi.png') # #### The open snowpit sites are a lot more spatially uniform in DOY estimates as expected. C-band radar will lose coherence in vegetated areas, so the backscatter minima at tree covered pixels will not necessarily be reflecting snowmelt properties. # ## And let's extract the backscatter timeseries over the buffered pixels for each site. # In[27]: PitIDs = ['COGMCO','COGMCT','COGMWO','COGMWT','COGMSO','COGMST'] f,ax=plt.subplots(3,2,figsize=(12,10)) for pitid, axs in zip(PitIDs,ax.ravel()): snowpit = snowpits_buffer[snowpits_buffer['PitID']==pitid] ts_ds_clipped_site = ts_ds_dem_masked.rio.clip(snowpit.geometry) axs.plot(ts_ds_clipped_site.time,10*np.log10(ts_ds_clipped_site).mean(dim=['x','y']),color='black',label='Mean Backscatter [Watts]') start_date = '2019-12-01' end_date = '2020-08-01' axs.set_xlim([datetime.strptime(start_date,'%Y-%m-%d'),datetime.strptime(end_date,'%Y-%m-%d')]) axs.set_ylabel('Backscatter [dB]') #runoff_map_clipped.plot(ax=axs,cmap='twilight',vmin=0,vmax=365) axs.set_title(f'{snowpit.Site.values[0]} \n Sentinel-1 Backscatter Timeseries') # to csv timeseries_time = ts_ds_clipped_site.time.values timeseries_values = 10*np.log10(ts_ds_clipped_site).mean(dim=['x','y']).values df = pd.DataFrame(timeseries_values,timeseries_time) df = df.reset_index() df.columns = ['datime','backcsatter_db'] df.to_csv(f'../output/grandmesa/snowpit_backscatter_means_csv/{snowpit.Site.values[0]}.csv') plt.tight_layout() plt.savefig('../output/grandmesa/plots/bs_snowpits.png') # #### We can clearly see backscatter minima in the open sites. # ## Let's add Sentinel-2 NSDI and NDWI. # In[28]: ndsi_dem_masked = s1_rtc_bs_utils.get_s2_ndsi(ts_ds_dem_masked) ndwi_dem_masked = s1_rtc_bs_utils.get_s2_ndwi(ts_ds_dem_masked) # In[29]: PitIDs = ['COGMCO','COGMCT','COGMWO','COGMWT','COGMSO','COGMST'] f,ax=plt.subplots(3,2,figsize=(12,10)) for pitid, axs in zip(PitIDs,ax.ravel()): snowpit = snowpits_buffer[snowpits_buffer['PitID']==pitid] ts_ds_clipped_site = ts_ds_dem_masked.rio.clip(snowpit.geometry) backscatter_line, = axs.plot(ts_ds_clipped_site.time,10*np.log10(ts_ds_clipped_site).mean(dim=['x','y']),color='black',label='S1 Mean Backscatter [dB]') start_date = '2019-12-01' end_date = '2020-08-01' axs.set_xlim([datetime.strptime(start_date,'%Y-%m-%d'),datetime.strptime(end_date,'%Y-%m-%d')]) axs.set_ylabel('Backscatter [dB]') #runoff_map_clipped.plot(ax=axs,cmap='twilight',vmin=0,vmax=365) axs.set_title(f'{snowpit.Site.values[0]}') ndsi_ax = axs.twinx() ndsi_clipped_site = ndsi_dem_masked.rio.clip(snowpit.geometry) ndsi_clipped_site = ndsi_clipped_site.where(ndsi_clipped_site<1) ndsi_line, = ndsi_ax.plot(ndsi_clipped_site.time,ndsi_clipped_site.mean(dim=['x','y']),color='green',label='S2 Mean NDSI') ndwi_clipped_site = ndwi_dem_masked.rio.clip(snowpit.geometry) ndwi_clipped_site = ndwi_clipped_site.where(ndsi_clipped_site<1) ndwi_line, = ndsi_ax.plot(ndwi_clipped_site.time,ndwi_clipped_site.mean(dim=['x','y']),color='blue',label='S2 Mean NWSI') #ndsi_ax.set_ylim([0,1]) axs.legend(handles=[backscatter_line,ndsi_line,ndwi_line]) plt.tight_layout() plt.savefig('../output/grandmesa/plots/bs_snowpits_ndsi_ndwi.png') # #### Interesting interplay between NDSI and NDWI. # ## Let's look at S1 Backscatter + S2 NDSI + Snowpit Data. # In[30]: PitIDs = ['COGMCO','COGMCT','COGMWO','COGMWT','COGMSO','COGMST'] f,ax=plt.subplots(3,2,figsize=(14,10)) for pitid, axs in zip(PitIDs,ax.ravel()): snowpit = snowpits_buffer[snowpits_buffer['PitID']==pitid] ts_ds_clipped_site = ts_ds_dem_masked.rio.clip(snowpit.geometry) backscatter_line, = axs.plot(ts_ds_clipped_site.time,10*np.log10(ts_ds_clipped_site).mean(dim=['x','y']),color='black',label='S1 Mean Backscatter') ts_log = 10*np.log10(ts_ds_clipped_site) axs.fill_between(ts_ds_clipped_site.time,ts_log.mean(dim=['x','y'])-ts_log.std(dim=['x','y']),ts_log.mean(dim=['x','y'])+ts_log.std(dim=['x','y']),alpha=0.1,color='black') start_date = '2019-12-01' end_date = '2020-08-01' axs.set_xlim([datetime.strptime(start_date,'%Y-%m-%d'),datetime.strptime(end_date,'%Y-%m-%d')]) axs.set_ylabel('Backscatter [dB]') #runoff_map_clipped.plot(ax=axs,cmap='twilight',vmin=0,vmax=365) axs.set_title(f'{snowpit.Site.values[0]}') ndsi_ax = axs.twinx() swe_ax = axs.twinx() ndsi_clipped_site = ndsi_dem_masked.rio.clip(snowpit.geometry) ndsi_clipped_site = ndsi_clipped_site.where(ndsi_clipped_site<1) ndsi_line, = ndsi_ax.plot(ndsi_clipped_site.time,ndsi_clipped_site.mean(dim=['x','y']),color='green',label='S2 Mean NDSI') swe_site = swe[swe['PitID'] == pitid] swe_scatter = swe_ax.scatter(x=swe_site['date'],y=swe_site['swe (mm)'],label='Snowpit SWE',c='blue') ndsi_ax.spines['right'].set_position(('outward', 40)) ndsi_ax.set_ylim([-1,1]) swe_ax.set_ylim([0,600]) ndsi_ax.spines["right"].set_visible(True) ndsi_ax.yaxis.set_label_position('right') ndsi_ax.yaxis.set_ticks_position('right') ndsi_ax.set_ylabel('NDSI') swe_ax.set_ylabel('SWE [mm]') swe_ax.yaxis.label.set_color('blue') ndsi_ax.yaxis.label.set_color(ndsi_line.get_color()) axs.tick_params(colors='black', which='both') swe_ax.tick_params(colors='blue', which='both') ndsi_ax.tick_params(colors='green', which='both') axs.legend(handles=[backscatter_line,ndsi_line,swe_scatter],loc='upper left') plt.tight_layout() plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.5, hspace=None) plt.savefig('../output/grandmesa/plots/bs_snowpits_ndsi_swe.png') # #### Everything seems to sync up nicely. We see minima in Sentinel-1 backscatter right as SWE seems to hit an inflection point. NDSI stays plateaued up until this point and drops significantly soon after. Something strange going on in the Mesa West Open data? Might want to check that out, maybe data processing error? # In[ ]: