Kevin

.gitignore

@@ -1,11 +1,9 @@
 # Ignore everything
 /*
 !.gitignore
-!its_live_web_app_tutorial.ipynb
-!/geoms 
-!/images
-!environment.yml
-!store_its_live_velocities_2022.json
-!ITS_LIVE API_tutorial.ipynb
-!its_live_datacube_term_poly_tutorial.ipynb
-!its_live_datacube_tutorial.ipynb
+!term_poly_gpkgs
+!/csv_outputs
+!/poly_sample
+!/term_poly_data_viz
+!/figs
+!gdf_concat.py

gdf_concat.py

@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Feb 13 17:40:43 2024
+
+@author: laserglaciers
+"""
+
+import geopandas as gpd
+import os
+
+glacier_names = ['KUJ', 'AVA','KAN','EQI']
+
+for name in glacier_names:
+    shp_file_path = f'/media/laserglaciers/upernavik/its_live_tutorial/{name}_polygons/'
+    os.chdir(shp_file_path)
+    gdfs_dates = [(gpd.read_file(fn), gpd.pd.to_datetime(fn[:10].replace('_','-'))) for fn in sorted(os.listdir(shp_file_path)) if fn.endswith('shp')] # 
+
+    gdfs, dates = zip(*gdfs_dates)
+
+    rdf = gpd.GeoDataFrame(gpd.pd.concat(gdfs, ignore_index=True), crs=gdfs[0].crs)
+    rdf['date'] = dates
+
+    if 'FID' in rdf.columns:
+        rdf.drop(columns=['FID'], inplace=True)
+        rdf.to_crs(crs=None,epsg=3413,inplace=True)
+    elif 'col1' in rdf.columns:
+        rdf.drop(columns=['col1'], inplace=True)
+        rdf.to_crs(crs=None,epsg=3413,inplace=True)
+
+    op = '/media/laserglaciers/upernavik/its_live_tutorial/merged_geoms/'
+    rdf.to_file(f'{op}{name}_merged.gpkg',driver='GPKG',crs="EPSG:3413")

poly_sample/datacube_poly_sample_multiprocess.py

@@ -0,0 +1,95 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Apr  4 15:45:50 2024
+
+@author: laserglaciers
+"""
+import os
+import json
+import geopandas as gpd
+from shapely.geometry import Point
+from shapely.geometry import Polygon
+import pandas as pd
+import numpy as np
+import xarray as xr
+import rioxarray as rxr
+import geoviews as gv
+import geoviews.feature as gf
+# import hvplot.pandas
+import s3fs
+import matplotlib.pyplot as plt
+import matplotlib.dates as mdates
+from datetime import datetime, timedelta
+from multiprocessing import Pool, cpu_count
+import dask
+from dask.distributed import Client, progress, LocalCluster
+import dask.array as da
+import geopandas as gpd
+from shapely import box
+from datacube_strain_rates_serial_function import *
+# %%
+
+def grab_cube(its_live_catalog, term_poly_path):
+
+    term_poly_path_abs = './term_poly_gpkgs/'
+    term_poly = gpd.read_file(f'{term_poly_path_abs}{term_poly_path}')
+    term_centroid = term_poly.centroid.to_crs(crs='EPSG:4326')
+    term_centroid = gpd.GeoDataFrame(geometry=term_centroid,crs=4326)
+    granule = its_live_catalog.sjoin(term_centroid, how='inner')
+    url = granule['zarr_url'].values[0]
+
+
+    cube = read_in_s3(url)
+    cube = cube.chunk({"mid_date": len(cube.mid_date)}) # setting the mid_date chunksize to the time dimensions helps A LOT
+    cube = cube.sortby('mid_date')
+    cube.rio.write_crs('EPSG:3413',inplace=True)
+
+    total_bounds = term_poly.total_bounds
+    bounding_box = box(*total_bounds)
+    cube = cube.rio.clip([bounding_box],term_poly.crs, drop=True,invert=False,all_touched=True)
+
+
+
+    cube_dates = cube.mid_date.dt.date.data # get the mid_dates
+    term_dates_pd = term_poly['date'].dt.date.values # get the polygon dates
+    matched_idx = np.where(np.isin(cube_dates,term_dates_pd))[0] # returns the indices of the matched dates so we can just use an isel (integer select) on the data cube
+    cube = cube.isel(mid_date=matched_idx) # use isel to grab the data we want
+
+    days = cube.date_dt.compute() # get the dt_array
+    days = days.data.astype('timedelta64[D]') # convert to days using .astype
+    day_filter = np.timedelta64(100, 'D') # the day separation we filter by
+    matched_date_dt_idx = np.where(days<=day_filter)[0] # get this indices 
+    cube = cube.isel(mid_date=matched_date_dt_idx) # use isel to grab the data we want
+    term_poly.to_crs(crs="EPSG:3413",inplace=True)
+
+    return cube, term_poly, term_poly_path[:3]
+
+
+catalog = gpd.read_file('https://its-live-data.s3.amazonaws.com/datacubes/catalog_v02.json') # ITS_LIVE data catalog where all the paths are defined
+term_poly_path = './term_poly_gpkgs/'
+
+terminus_poly_list = [file for file in os.listdir(term_poly_path) if file.endswith('gpkg')]
+
+cube_list = []
+for term_poly in terminus_poly_list:
+    filtered_cube = grab_cube(catalog,term_poly)
+    cube_list.append(filtered_cube)
+
+# %%
+
+for data in cube_list:
+
+    poly_sample(*data)
+
+# %%
+def worker(cube, term_poly,glacier_abbr):
+    proc = poly_sample(cube, term_poly,glacier_abbr)
+    proc.wait()
+
+
+if __name__ == "__main__":
+    with Pool(processes=cpu_count()-6) as pool:
+          results = pool.starmap(worker, cube_list[:2])
+
+