# need to upgrade odc-stac to fix keyerror issue for some tiles
!pip install -q odc-stac -U
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
import pystac
import pystac_client
import stackstac
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
from rechunker import rechunk
sys.path.append('../sar_snowmelt_timing')
import s1_rtc_bs_utils
import contextily as ctx
import rioxarray as rxr
import skimage
import pathlib
import glob
import re
import time
import fsspec
cluster = dask_gateway.GatewayCluster()
client = cluster.get_client()
cluster.adapt(minimum=90, maximum=100)
print(client.dashboard_link)
https://pccompute.westeurope.cloudapp.azure.com/compute/services/dask-gateway/clusters/prod.ec6eee513a0f48d6bf281eafb5d2a9e4/status
#from dask.distributed import Client, LocalCluster
#cluster = LocalCluster(processes=False, n_workers=4, local_directory='/tmp')
#client = Client(cluster)
#client
url = 'https://github.com/scottyhq/mgrs/raw/main/MGRS_LAND.parquet' # Scott created an MGRS parquet file here https://github.com/scottyhq/mgrs
with fsspec.open(url) as file:
mgrs_gdf = gpd.read_parquet(file)
#tile_names_10T = list(mgrs_gdf[mgrs_gdf.tile.str.match('10T[DEFG]')].tile)
tile_names = open('../input/MGRS_square_list/westernUS_square_list.txt','r').readlines()
tile_names = [i.rstrip('\n') for i in tile_names]
# comment out if want more than 10T
#tile_names = [tile for tile in tile_names if tile in tile_names_10T]
years = [2015,2016,2017,2018,2019,2020,2021,2022]
resolution = 80 #40
snow_classes = [ # https://zenodo.org/record/2626737#.ZDMMf3bMIQ8
# 0, #Little-to-no snow
1, #Indeterminate due to clouds
2, #Ephemeral snow
3, #Seasonal snow
]
classes = [ # page 13 of https://esa-worldcover.s3.amazonaws.com/v100/2020/docs/WorldCover_PUM_V1.0.pdf
# 10, # treecover
20, # shrubland
30, # grassland
40, # cropland
# 50, # built-up
60, #bare / sparse vegetation
70, # snow and ice
# 80, # permanent water bodies
90, # herbaceous wetlands
95, # mangroves
100 # loss and lichen
]
# for nm in tile_names:
# print(f'Processing {nm}:')
# bbox_gdf = mgrs_gdf[mgrs_gdf.tile==nm]
# ts_ds = s1_rtc_bs_utils.get_s1_rtc_stac_odc_pc(bbox_gdf,start_time='2014-01-01',end_time='2022-12-31',resolution=resolution,epsg=f'EPSG:{bbox_gdf.epsg.values[0]}')
# worldcover = s1_rtc_bs_utils.get_worldcover(ts_ds)
# ts_ds = ts_ds.where(worldcover.isin(classes))
# path = f'../output/MGRS/{nm}'
# savedir = pathlib.Path(path)
# savedir.mkdir(parents=True,exist_ok=True)
# for year in years:
# print(f'Processing {year}...')
# tic = time.perf_counter()
# one_year = slice(f'{year}-01-01',f'{year}-12-31')
# ts_ds_oneyear = ts_ds.sel(time=one_year)
# runoffs_median = s1_rtc_bs_utils.get_runoff_onset(ts_ds_oneyear)
# runoffs_median_computed = runoffs_median.compute()
# runoffs_median_computed.rio.write_crs(rio.CRS.from_epsg(bbox_gdf['epsg']),inplace=True)
# runoffs_median_computed.dt.dayofyear.rio.to_raster(f'../output/MGRS/{nm}/runoff_onset_{nm}_{year}_median_{resolution}m.tif')
# toc = time.perf_counter()
# print(f'Processing time: {toc - tic:0.1f} seconds')
# print(f'Done!')
# for nm in tile_names:
# print(f'Processing {nm}:')
# geotiff_list = glob.glob(f'../output/MGRS/{nm}/runoff_onset_*.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')
# runoffs_allyears.median(dim='time').rio.to_raster(f'../output/MGRS/{nm}/runoff_onset_{nm}_allyears_median_{resolution}m.tif')
# runoffs_allyears.std(dim='time').rio.to_raster(f'../output/MGRS/{nm}/runoff_onset_{nm}_allyears_std_{resolution}m.tif')
# polyfits= runoffs_allyears.polyfit('time',deg=1,full=True,cov=True)
# polyfits.polyfit_coefficients[0].rio.set_crs(runoffs_allyears.rio.crs).rio.to_raster(f'../output/MGRS/{nm}/runoff_onset_{nm}_allyears_polyfit_{resolution}m.tif')
# xr.corr(runoffs_allyears.time,runoffs_allyears,dim='time').rio.set_crs(runoffs_allyears.rio.crs).rio.to_raster(f'../output/MGRS/{nm}/runoff_onset_{nm}_allyears_corr_strength_{resolution}m.tif')
#bbox_gdf = gpd.read_file('../input/shapefiles/western_us.geojson')
#ts_ds = s1_rtc_bs_utils.get_s1_rtc_stac_pc(bbox_gdf,start_time='2020-01-01',end_time='2020-12-31',resolution=80)
#ts_ds
import getpass
import azure.storage.blob
import io
#sas_token = getpass.getpass() # prompts for the sas_token
container_client = azure.storage.blob.ContainerClient(
"https://snowmelt.blob.core.windows.net",
container_name="snowmelt",
credential=sas_token,
)
file_list = [item['name'] for item in list(container_client.list_blobs())]
for nm in tile_names:
if len(list(filter(lambda x: nm in x, list(filter(lambda x: str(resolution) in x, file_list))))) < 6:
print(f'Processing {nm}:')
bbox_gdf = mgrs_gdf[mgrs_gdf.tile==nm]
ts_ds = s1_rtc_bs_utils.get_s1_rtc_stac_odc_pc(bbox_gdf,start_time='2015-01-01',end_time='2022-12-31',resolution=resolution,epsg=f'EPSG:{bbox_gdf.epsg.values[0]}')
ts_ds = ts_ds.sel(band='vv') # remove for regular processing ######################################################
worldcover = s1_rtc_bs_utils.get_worldcover(ts_ds)
snow_mask = s1_rtc_bs_utils.get_snowmask(ts_ds)
ts_ds = ts_ds.where(worldcover.isin(classes))
ts_ds = ts_ds.where(snow_mask.isin(snow_classes))
for year in years:
print(f'Processing {year}...')
try:
tic = time.perf_counter()
#one_year = slice(f'{year}-01-01',f'{year}-12-31') # DEFINE MELT PERIOD HERE!!! TRY MARCH1-AUG1???
melt_year = slice(f'{year}-02-01',f'{year}-07-31')
ts_ds_oneyear = ts_ds.sel(time=melt_year)
runoffs_median = s1_rtc_bs_utils.get_runoff_onset(ts_ds_oneyear)
runoffs_median_computed = runoffs_median.compute()
runoffs_median_computed.rio.write_crs(rio.CRS.from_epsg(bbox_gdf['epsg']),inplace=True)
with io.BytesIO() as buffer:
runoffs_median_computed_doy = runoffs_median_computed.dt.dayofyear.where(runoffs_median_computed.dt.dayofyear>0).rio.write_nodata(-32768,encoded=True)
runoffs_median_computed_doy.rio.to_raster(buffer, driver="COG", dtype='int16')#dtype='int16' #compress='LZW' tiled=True, dytpe=''
buffer.seek(0)
blob_client = container_client.get_blob_client(f"eric/MGRS/{nm}/runoff_onset_{nm}_{year}_median_{resolution}m.tif")
blob_client.upload_blob(buffer, overwrite=True)
toc = time.perf_counter()
print(f'Processing time: {toc - tic:0.1f} seconds')
except:
print(f'Failed.')
print(f'Done!')
Done!
# only for azure blob storage
file_list = [item['name'] for item in list(container_client.list_blobs())]
url = "https://snowmelt.blob.core.windows.net/snowmelt/"
for nm in tile_names:
if len(list(filter(lambda x: nm in x, list(filter(lambda x: str(resolution) in x, file_list))))) < 12:
print(f'Processing {nm}:')
geotiff_list = list(filter(lambda x: nm in x, list(filter(lambda x: str(resolution) in x, file_list))))
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(f"{url}{i}") for i in geotiff_list],dim=time_var).squeeze().sortby('time')
runoffs_allyears = runoffs_allyears.where(runoffs_allyears>-32768)
#runoffs_allyears.plot(col='time')
with io.BytesIO() as buffer:
runoffs_allyears_doy = runoffs_allyears.median(dim='time').where(runoffs_allyears.median(dim='time')>0).rio.write_nodata(-32768,encoded=True)
runoffs_allyears_doy.rio.to_raster(buffer, driver="COG",dtype='int16')
buffer.seek(0)
blob_client = container_client.get_blob_client(f"eric/MGRS/{nm}/runoff_onset_{nm}_allyears_median_{resolution}m.tif")
blob_client.upload_blob(buffer, overwrite=True)
with io.BytesIO() as buffer:
runoffs_allyears.std(dim='time').astype('float32').rio.to_raster(buffer, driver="COG",dtype='float32')
buffer.seek(0)
blob_client = container_client.get_blob_client(f"eric/MGRS/{nm}/runoff_onset_{nm}_allyears_std_{resolution}m.tif")
blob_client.upload_blob(buffer, overwrite=True)
polyfits= runoffs_allyears.polyfit('time',deg=1,full=True,cov=True)
with io.BytesIO() as buffer:
polyfits.polyfit_coefficients[0].astype('float32').rio.set_crs(runoffs_allyears.rio.crs).rio.to_raster(buffer, driver="COG",dtype='float32')
buffer.seek(0)
blob_client = container_client.get_blob_client(f"eric/MGRS/{nm}/runoff_onset_{nm}_allyears_polyfit_{resolution}m.tif")
blob_client.upload_blob(buffer, overwrite=True)
with io.BytesIO() as buffer:
xr.corr(runoffs_allyears.time,runoffs_allyears,dim='time').astype('float32').rio.set_crs(runoffs_allyears.rio.crs).rio.to_raster(buffer, driver="COG", dtype='float32')
buffer.seek(0)
blob_client = container_client.get_blob_client(f"eric/MGRS/{nm}/runoff_onset_{nm}_allyears_corr_strength_{resolution}m.tif")
blob_client.upload_blob(buffer, overwrite=True)
print(f'Done!')
Processing 13UDP: Processing 12UVU: Processing 12UYU: Processing 13UEP: Processing 13TDN: Processing 13TEM: Processing 13TEN: Processing 13TFM: Processing 12UWU: Processing 13UCP: Processing 13TFL: Processing 13TEL: Processing 13TDM: Processing 12UXU: Processing 13TCN: Processing 13UFP: Processing 12TYN: Processing 13TFN: Processing 13TDK: Processing 12TXM: Processing 12TYT: Processing 13TBG: Processing 12TYP: Processing 12TYM: Processing 13TDL: Processing 13TCH: Processing 13TCM: Processing 12TUP: Processing 13TDJ: Processing 12TYL: Processing 13TBF: Processing 12TXT: Processing 13TCG: Processing 13TCJ: Processing 13TCL: Processing 12TWM: Processing 12TWT: Processing 11TNG: Processing 12TVT: Processing 12TYR: Processing 12TXS: Processing 12TXQ: Processing 12TXR: Processing 12TXP: Processing 11TPF: Processing 11TQJ: Processing 12TUQ: Processing 13TCK: Processing 12TYS: Processing 13TDH: Processing 12TWS: Processing 12TXL: Processing 12TXN: Processing 13TEK: Processing 12TYQ: Processing 11TPG: Processing 13TDG: Processing 13TFK: Processing 12UUU: Processing 11TQK: Processing 12TTM: Processing 11TNF: Processing 11TQG: Processing 12TVP: Processing 12TVS: Processing 12TWN: Processing 12TVR: Processing 12TWL:
2023-04-10 17:31:37,273 - distributed.client - ERROR - Failed to reconnect to scheduler after 30.00 seconds, closing client Exception in callback None() handle: <Handle cancelled> Traceback (most recent call last): File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 1391, in _do_ssl_handshake self.socket.do_handshake() File "/srv/conda/envs/notebook/lib/python3.10/ssl.py", line 1342, in do_handshake self._sslobj.do_handshake() ssl.SSLEOFError: EOF occurred in violation of protocol (_ssl.c:997) During handling of the above exception, another exception occurred: Traceback (most recent call last): File "/srv/conda/envs/notebook/lib/python3.10/asyncio/events.py", line 80, in _run self._context.run(self._callback, *self._args) File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/platform/asyncio.py", line 189, in _handle_events handler_func(fileobj, events) File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 696, in _handle_events self._handle_read() File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 1478, in _handle_read self._do_ssl_handshake() File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 1400, in _do_ssl_handshake return self.close(exc_info=err) File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 611, in close self._signal_closed() File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 641, in _signal_closed self._ssl_connect_future.exception() asyncio.exceptions.CancelledError Exception in callback None() handle: <Handle cancelled> Traceback (most recent call last): File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 1391, in _do_ssl_handshake self.socket.do_handshake() File "/srv/conda/envs/notebook/lib/python3.10/ssl.py", line 1342, in do_handshake self._sslobj.do_handshake() ssl.SSLEOFError: EOF occurred in violation of protocol (_ssl.c:997) During handling of the above exception, another exception occurred: Traceback (most recent call last): File "/srv/conda/envs/notebook/lib/python3.10/asyncio/events.py", line 80, in _run self._context.run(self._callback, *self._args) File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/platform/asyncio.py", line 189, in _handle_events handler_func(fileobj, events) File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 696, in _handle_events self._handle_read() File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 1478, in _handle_read self._do_ssl_handshake() File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 1400, in _do_ssl_handshake return self.close(exc_info=err) File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 611, in close self._signal_closed() File "/srv/conda/envs/notebook/lib/python3.10/site-packages/tornado/iostream.py", line 641, in _signal_closed self._ssl_connect_future.exception() asyncio.exceptions.CancelledError
Processing 13TDF: Processing 12TWR: Processing 12TYK: Processing 12TVM: Processing 12TUR: Processing 13TBE: Processing 12TVN: Processing 13SDB: Processing 11TPE: Processing 13TEJ: Processing 12UYV: Processing 11TPJ: Processing 12UXV: Processing 13UCQ: Processing 11TMG: Processing 12TUN: Processing 11TNE: Processing 13UFQ: Processing 10TGM: Processing 12UWV: Processing 13TCF: Processing 13UDQ: Processing 12UVV: Processing 13UEQ: Processing 11TKG: Processing 12TXK: Processing 12TVQ: Processing 13TCE: Processing 11SND: Processing 12TUT: Processing 12UUV: Processing 12TWP: Processing 13SCC: Processing 12TWK: Processing 13SCB: Processing 11TLG: Processing 12TVL: Processing 12TUS: Processing 11TPK: Processing 11SPD: Processing 11TQL: Processing 13SBC: Processing 13SDD: Processing 13TGL: Processing 12TUM: Processing 13SDC: Processing 11UQQ: Processing 12TWQ: Processing 12SVJ: Processing 11TMK: Processing 12TVK: Processing 12SYH: Processing 12SYJ: Processing 13TGM: Processing 11UQP: Processing 11TQM: Processing 13TGK: Processing 13SCD: Processing 13SBD: Processing 10TGL: Processing 11TKF: Processing 10UGU: Processing 10TGP: Processing 11TNH: Processing 13TEH: Processing 13SBB: Processing 11TML: Processing 11SLB: Processing 13TGN: Processing 11TQN: Processing 12SYG: Processing 12SVH: Processing 11TNK: Processing 12SXJ: Processing 13UGP: Processing 11TQH: Processing 11TPL: Processing 11TLJ: Processing 13TDE: Processing 13SDA: Processing 11TNL: Processing 10TGN: Processing 11SMD: Processing 10UFU: Processing 11ULP: Processing 13SCA: Processing 11SKC: Processing 12SWJ: Processing 11TLH: Processing 11TLK: Processing 12SUH: Processing 13TFJ: Processing 11SLC: Processing 13SED: Processing 11TQF: Processing 12SUG: Processing 11TPH: Processing 10TGQ: Processing 11TPM: Processing 11TMF: Processing 11UPP: Processing 11SLA: Processing 10UGV: Processing 10TFP: Processing 10TFT: Processing 11TME: Processing 11SQD: Processing 12TTL: Processing 13UGQ: Processing 10SGH: Processing 10TFN: Processing 10UFV: Processing 11TQE: Processing 11TMJ: Processing 11TPN: Processing 11TNJ: Processing 11TMH: Processing 11ULQ: Processing 10TGT: Processing 12SXH: Processing 12STJ: Processing 12TUK: Processing 10TFS: Processing 11SNC: Processing 12TTK: Processing 11TLF: Processing 11TNM: Processing 13TEF: Processing 10TFM: Processing 12SXG: Processing 11SQC: Processing 10SGJ: Processing 11UPQ: Processing 11TLL: Processing 11SMC: Processing 11SKB: Processing 11SKD: Processing 12STH: Processing 12SUJ: Processing 11SPC: Processing 13TGJ: Processing 10TGK: Processing 11TKE: Processing 10TFL: Processing 11TNN: Processing 11UNP: Processing 11UMP: Processing 13SEB: Processing 12SVG: Processing 10TES: Processing 13SCV: Processing 13TEG: Processing 13SEA: Processing 12TUL: Processing 11UNQ: Processing 10TEQ: Processing 10UEV: Processing 12SXC: Processing 10TDT: Processing 10TEP: Processing 10TEN: Processing 11TMM: Processing 10UDU: Processing 13SBA: Processing 11UMQ: Processing 10TFK: Processing 10TEL: Processing 11TLE: Processing 13SDV: Processing 10TDL: Processing 12SWH: Processing 10TGR: Processing 11SMB: Processing 10TFQ: Processing 10TET: Processing 13TFH: Processing 10UEU: Processing 10TDM: Processing 12SWG: Processing 11SMA: Processing 10TGS: Processing 10TEM: Processing 10SGG: Processing 10TFR: Processing 12STG: Processing 10TER: Processing 11TLN: Processing 11SLD: Processing 12SXF: Processing 12SYF: Processing 13SEC: Processing 12SUF: Processing 13TGH: Processing 11SLV: Processing 12SYB: Processing 10SFJ: Processing 10UCV: Processing 12SXE: Processing 11SPA: Processing 13SBV: Processing 11SNT: Processing 12SYC: Processing 11TMN: Processing 10UDV: Processing 12SVF: Done!