|
32 | 32 | merge-multiple: true
|
33 | 33 |
|
34 | 34 | - name: Compute summary statistics
|
35 |
| - shell: bash -el -c "python -u {0}" |
36 | 35 | run: |
|
37 |
| - import xarray as xr |
38 |
| - import rasterio as rio |
39 |
| - import rioxarray |
40 |
| - from glob import glob |
41 |
| - from datetime import datetime |
42 |
| - import seaborn as sns |
43 |
| - import matplotlib.pyplot as plt |
44 |
| -
|
45 |
| - # functions to load tifs to xarray |
46 |
| - def xr_read_geotif(geotif_file_path, masked=True): |
47 |
| - da = rioxarray.open_rasterio(geotif_file_path, masked=True) |
48 |
| - # Extract bands and assign as variables in xr.Dataset() |
49 |
| - ds = xr.Dataset() |
50 |
| - for i, v in enumerate(da.band): |
51 |
| - da_tmp = da.sel(band=v) |
52 |
| - da_tmp.name = "band" + str(i + 1) |
53 |
| - ds[da_tmp.name] = da_tmp |
54 |
| - # Delete empty band coordinates. |
55 |
| - del ds.coords["band"] |
56 |
| - |
57 |
| - return ds |
58 |
| - |
59 |
| - def combine_ds(data_dir, file_type='horizontal_velocity'): |
60 |
| - datasets = [] |
61 |
| - tif_list = glob(f'{data_dir}/S2*{file_type}.tif') |
62 |
| - for tif_path in tif_list: |
63 |
| - dates = tif_list[0].split('/')[-1][3:20] #parse filename for dates |
64 |
| - start_date = datetime.strptime(dates[:8], '%Y%m%d') |
65 |
| - end_date = datetime.strptime(dates[-8:], '%Y%m%d') |
66 |
| - t_baseline = end_date - start_date |
67 |
| - |
68 |
| - src = xr_read_geotif(tif_path, masked=False) #read product to xarray ds |
69 |
| - src = src.assign_coords({"dates": dates}) |
70 |
| - src = src.expand_dims("dates") |
71 |
| - |
72 |
| - src = src.assign_coords(start_date = ('dates', [start_date])) |
73 |
| - src = src.assign_coords(end_date = ('dates', [end_date])) |
74 |
| - src = src.assign_coords(t_baseline = ('dates', [t_baseline])) |
75 |
| - |
76 |
| - src = src.rename({'band1':file_type}) |
77 |
| - |
78 |
| - datasets.append(src) |
79 |
| - |
80 |
| - ds = xr.concat(datasets, dim="dates", combine_attrs="no_conflicts") #create dataset |
81 |
| - ds = ds.sortby('dates') |
82 |
| - |
83 |
| - return ds |
84 |
| -
|
85 |
| - # read in tifs |
86 |
| - veloc_ds = combine_ds(data_dir='glacier_image_correlation', file_type='horizontal_velocity') |
87 |
| - # calculate and save median velocity |
88 |
| - veloc_da_median = veloc_ds.horizontal_velocity.median(dim='dates') |
89 |
| - # veloc_da_median.rio.to_raster('glacier_image_correlation/median_horizontal_velocity.tif') |
90 |
| - # save standard deviation of velocity |
91 |
| - veloc_da_stdev = veloc_ds.horizontal_velocity.std(dim='dates') |
92 |
| - # veloc_da_stdev.rio.to_raster('glacier_image_correlation/stdev_horizontal_velocity.tif') |
93 |
| - # save valid velocity pixel count |
94 |
| - veloc_da_count = veloc_ds.horizontal_velocity.count(dim='dates') |
95 |
| - # veloc_da_count.rio.to_raster('glacier_image_correlation/count_horizontal_velocity.tif') |
96 |
| -
|
97 |
| - # plot summary statistics |
98 |
| - sns.set_theme() |
99 |
| - f, ax = plt.subplots(1, 3, figsize=(15, 5), sharex=True, sharey=True) |
100 |
| - veloc_da_median.plot(ax=ax[0], vmin=0, vmax=600, cmap='inferno', cbar_kwargs= {'shrink':0.7, 'label':'velocity (m/yr)'}) |
101 |
| - veloc_da_stdev.plot(ax=ax[1], vmin=0, vmax=300, cmap='cividis', cbar_kwargs= {'shrink':0.7, 'label':'standard deviation (m/yr)'}) |
102 |
| - veloc_da_count.plot(ax=ax[2], vmin=0, cmap='Blues', cbar_kwargs= {'shrink':0.7, 'label':'pixel count'}) |
103 |
| - ax[0].set_aspect('equal') |
104 |
| - ax[1].set_aspect('equal') |
105 |
| - ax[2].set_aspect('equal') |
106 |
| - ax[0].set_title(f'median velocity') |
107 |
| - ax[1].set_title(f'velocity standard deviation') |
108 |
| - ax[2].set_title(f'valid pixel count') |
109 |
| - ax[0].set_xticks([]) |
110 |
| - ax[0].set_yticks([]) |
111 |
| - ax[1].set_xticks([]) |
112 |
| - ax[1].set_yticks([]) |
113 |
| - ax[2].set_xticks([]) |
114 |
| - ax[2].set_yticks([]) |
115 |
| - ax[0].set_xlabel('') |
116 |
| - ax[0].set_ylabel('') |
117 |
| - ax[1].set_xlabel('') |
118 |
| - ax[1].set_ylabel('') |
119 |
| - ax[2].set_xlabel('') |
120 |
| - ax[2].set_ylabel('') |
121 |
| - f.tight_layout() |
122 |
| - |
123 |
| - f.savefig('glacier_image_correlation/velocity_summary_statistics.png', dpi=300) |
| 36 | + python glacier_image_correlation/summary_statistics.py |
124 | 37 |
|
125 | 38 | - name: Upload summary statistics
|
126 | 39 | uses: actions/upload-artifact@v4
|
|
0 commit comments