Extract velocity map data at GNSS locations
Extract velocity map data at GNSS locations¶
This script samples velocity at GNSS locations and updates all pt* fields in notebooks/results_2022.csv.
To reproduce this workflow, make sure you have downloaded all necessary input files (velocity maps and static terrain geometries) from https://doi.org/10.17605/OSF.IO/HE7YR and have updated the Vx and Vy columns in notebooks/results_2022.csv with the downloaded file paths before starting the analysis.
from glaft.georaster import Raster
import rasterio
import pandas as pd
import geopandas as gpd
import numpy as np
df = pd.read_csv('../results_2022.csv', dtype=str)
# df
The cell below provides a sanity check for glaft.georaster.Raster’s value_at_coords method.
tmp = Raster(df.loc[0, 'Vx'])
a, b = tmp.value_at_coords(621306.41954208, 6738829.50233354, window=3, return_window=True)
vx_grid = rasterio.open(df.loc[0, 'Vx'])
sample_gen_vx = vx_grid.sample([(621306.41954208, 6738829.50233354)])
vx_sampled = np.array([float(record) for record in sample_gen_vx])
Now let’s start the analysis. If you download these files from https://doi.org/10.17605/OSF.IO/HE7YR, make sure to change the paths to the correct file locations on your local machine.
gps_files = ['/home/jovyan/Projects/PX_comparison/GPS/Kaskawulsh_2018-04-05_to_2018-03-04_GPS',
'/home/jovyan/Projects/PX_comparison/GPS/Kaskawulsh_2018-08-18_to_2018-08-02_GPS',
'/home/jovyan/Projects/PX_comparison/GPS/Kaskawulsh_2018-03-14_to_2018-03-04_GPS',
'/home/jovyan/Projects/PX_comparison/GPS/Kaskawulsh_2018-06-27_to_2018-05-08_GPS']
datestrs = ['LS8-20180304-20180405', 'LS8-20180802-20180818', 'Sen2-20180304-20180314', 'Sen2-20180508-20180627']
datenums = [32, 16, 10, 50]
Steps here:
Get and print the UTM coordinates of three GNSS statations for each scene pair.
Sample every velocity maps.
Create additional fields and calculate the difference between GNSS and feature tracked measurements.
for gps_file, datestr, datenum in zip(gps_files, datestrs, datenums):
gps = pd.read_csv(gps_file)
# Additional treatment for Sen2-20180508-20180627
# many of the points here should be (nan, nan), (nan, nan) but nans does not work with rio.sample
if datestr == 'Sen2-20180508-20180627':
gps.loc[1, 'end_easting'] = 610481.2868266493
gps.loc[1, 'end_northing'] = 6737102.953712379
gps.loc[2, 'end_easting'] = 601790.4387747
gps.loc[2, 'end_northing'] = 6733753.77267354
gps = gps.loc[0:2]
gps = gpd.GeoDataFrame(gps, geometry=gpd.points_from_xy(gps['end_easting'], gps['end_northing']), crs='EPSG:32607')
# This is beginning coordinates
gps_xy = list(gps[['end_easting', 'end_northing']].to_records(index=False))
print(datestr, gps_xy)
gps['vx (m/d)'] = (gps['start_easting'] - gps['end_easting']) / datenum
gps['vy (m/d)'] = (gps['start_northing'] - gps['end_northing']) / datenum
df_s = df.loc[df['Date'] == datestr]
for idx, row in df_s.iterrows():
vx_grid = Raster(row.Vx)
vy_grid = Raster(row.Vy)
sampled = []
for x, y in gps_xy:
vx_avg, vx_3by3 = vx_grid.value_at_coords(x, y, window=3, return_window=True)
vy_avg, vy_3by3 = vy_grid.value_at_coords(x, y, window=3, return_window=True)
vx_3by3[vx_3by3 < -9998] = np.nan
vy_3by3[vy_3by3 < -9998] = np.nan
vx_3by3[vx_3by3 == 0.0] = np.nan #Vmap
vy_3by3[vy_3by3 == 0.0] = np.nan #Vmap
vx_nn = vx_3by3[0, 1, 1] # nearest neighbor value
vy_nn = vy_3by3[0, 1, 1]
if np.any(~np.isnan(vx_3by3)):
vx_avg = np.nanmean(vx_3by3)
else:
vx_avg = np.nan
if np.any(~np.isnan(vy_3by3)):
vy_avg = np.nanmean(vy_3by3)
else:
vy_avg = np.nan
vx_nn -= float(df.loc[idx, 'SAV-peak-x'])
vx_avg -= float(df.loc[idx, 'SAV-peak-x'])
vy_nn -= float(df.loc[idx, 'SAV-peak-y'])
vy_avg -= float(df.loc[idx, 'SAV-peak-y'])
sampled.append([vx_nn, vx_avg, vy_nn, vy_avg])
sampled = np.array(sampled)
# print(row.Vx, float(df.loc[idx, 'SAV-peak-x']), float(df.loc[idx, 'SAV-peak-y']), sampled)
df.loc[idx, 'pt0_vxdiff'] = sampled[0, 0] - gps.loc[0, 'vx (m/d)']
df.loc[idx, 'pt0_vxavgdiff'] = sampled[0, 1] - gps.loc[0, 'vx (m/d)']
df.loc[idx, 'pt0_vydiff'] = sampled[0, 2] - gps.loc[0, 'vy (m/d)']
df.loc[idx, 'pt0_vyavgdiff'] = sampled[0, 3] - gps.loc[0, 'vy (m/d)']
df.loc[idx, 'pt1_vxdiff'] = sampled[1, 0] - gps.loc[1, 'vx (m/d)']
df.loc[idx, 'pt1_vxavgdiff'] = sampled[1, 1] - gps.loc[1, 'vx (m/d)']
df.loc[idx, 'pt1_vydiff'] = sampled[1, 2] - gps.loc[1, 'vy (m/d)']
df.loc[idx, 'pt1_vyavgdiff'] = sampled[1, 3] - gps.loc[1, 'vy (m/d)']
df.loc[idx, 'pt2_vxdiff'] = sampled[2, 0] - gps.loc[2, 'vx (m/d)']
df.loc[idx, 'pt2_vxavgdiff'] = sampled[2, 1] - gps.loc[2, 'vx (m/d)']
df.loc[idx, 'pt2_vydiff'] = sampled[2, 2] - gps.loc[2, 'vy (m/d)']
df.loc[idx, 'pt2_vyavgdiff'] = sampled[2, 3] - gps.loc[2, 'vy (m/d)']
LS8-20180304-20180405 [(621306.41954208, 6738829.50233354), (610435.5249175, 6737129.57698521), (601733.22946583, 6733710.66504834)]
LS8-20180802-20180818 [(621363.01607688, 6738895.12164604), (610506.52739125, 6737089.56006354), (601790.43877479, 6733753.77267354)]
Sen2-20180304-20180314 [(621306.41954208, 6738829.50233354), (610435.5249175, 6737129.57698521), (601733.22946583, 6733710.66504834)]
Sen2-20180508-20180627 [(621324.96198502, 6738852.60218059), (610481.28682665, 6737102.95371238), (601790.4387747, 6733753.77267354)]
You can comment/uncomment these lines to examine the data/results.
gps
# df
| Unnamed: 0 | date1 | date2 | start_easting | start_northing | end_easting | end_northing | distance_traveled (m) | velocity (m/d) | geometry | vx (m/d) | vy (m/d) | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 2018-06-27 | 2018-05-08 | 621348.115449 | 6.738880e+06 | 621324.961985 | 6.738853e+06 | 36.045525 | 0.720837 | POINT (621324.962 6738852.602) | 0.463069 | 0.55252 |
| 1 | 1 | 2018-06-27 | 2018-05-08 | 610481.286827 | 6.737103e+06 | 610481.286827 | 6.737103e+06 | NaN | NaN | POINT (610481.287 6737102.954) | 0.000000 | 0.00000 |
| 2 | 2 | 2018-06-27 | 2018-05-08 | NaN | NaN | 601790.438775 | 6.733754e+06 | NaN | NaN | POINT (601790.439 6733753.773) | NaN | NaN |
df.to_csv('../results_2022.csv', index=False)