ITS_LIVE: Extract velocity map data at GNSS locations
ITS_LIVE: Extract velocity map data at GNSS locationsΒΆ
This script samples velocity at GNSS locations and updates all pt* fields in notebooks/results_ITSLIVE.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_ITSLIVE.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
from datetime import datetime
from pyproj import Transformer
df = pd.read_csv('../results_ITSLIVE.csv', dtype=str)
# df
# All GNSS coordinates are in EPSG 32607 and have to be reprojected to EPSG3413 before sampling the geotiffs.
transformer = Transformer.from_crs("epsg:32607", "epsg:3413")
The cell below is the main procedure.
GPS_root = '/home/jovyan/Projects/PX_comparison/GPS/'
for idx, row in df.iterrows():
startdate = datetime.strptime(row['Start date'], '%Y%m%d')
enddate = datetime.strptime(row['End date'], '%Y%m%d')
timedel = enddate - startdate
duration = timedel.days / 365 # in yrs
startdate_gpsstr = startdate.strftime('%Y-%m-%d')
enddate_gpsstr = enddate.strftime('%Y-%m-%d')
gps_file = GPS_root + 'Kaskawulsh_{}_to_{}_GPS'.format(enddate_gpsstr, startdate_gpsstr)
gps = pd.read_csv(gps_file)
######## Prep corrdinates in EPSG3413
for idx2, row2 in gps.iterrows():
x32607 = row2.start_easting
y32607 = row2.start_northing
x3413, y3413 = transformer.transform(x32607, y32607)
gps.loc[idx2, 'start_easting_3413'] = x3413
gps.loc[idx2, 'start_northing_3413'] = y3413
x32607 = row2.end_easting
y32607 = row2.end_northing
x3413, y3413 = transformer.transform(x32607, y32607)
gps.loc[idx2, 'end_easting_3413'] = x3413
gps.loc[idx2, 'end_northing_3413'] = y3413
#########
# This is beginning coordinates
gps = gpd.GeoDataFrame(gps, geometry=gpd.points_from_xy(gps['end_easting_3413'], gps['end_northing_3413']), crs='EPSG:3413')
gps_xy = list(gps[['end_easting_3413', 'end_northing_3413']].to_records(index=False))
gps['vx (m/yr)'] = (gps['start_easting_3413'] - gps['end_easting_3413']) / duration
gps['vy (m/yr)'] = (gps['start_northing_3413'] - gps['end_northing_3413']) / duration
gps['v (m/yr)'] = np.abs(gps['velocity (m/d)'] * 365)
vx_grid = Raster(row.Vx)
vy_grid = Raster(row.Vy)
v_grid = Raster(row.Vx.replace('vx', 'v'))
verr_grid = Raster(row.Vx.replace('vx', 'v_error'))
sampled = []
sampled2 = []
for x, y in gps_xy:
# print(gps_file, x, y)
if np.isnan(x) or np.isnan(y) or np.isinf(x) or np.isinf(y):
sampled.append([np.nan, np.nan, np.nan, np.nan])
sampled2.append([np.nan, np.nan])
else:
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 # ITSLIVE nodata = -32767
vy_3by3[vy_3by3 < -9998] = np.nan # ITSLIVE nodata = -32767
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
sampled.append([vx_nn, vx_avg, vy_nn, vy_avg])
v_avg, v_3by3 = v_grid.value_at_coords(x, y, window=3, return_window=True)
v_3by3[v_3by3 < -9998] = np.nan # ITSLIVE nodata = -32767
if np.any(~np.isnan(v_3by3)):
v_avg = np.nanmean(v_3by3)
else:
v_avg = np.nan
verr_avg, verr_3by3 = verr_grid.value_at_coords(x, y, window=3, return_window=True)
verr_3by3[verr_3by3 < -9998] = np.nan # ITSLIVE nodata = -32767
if np.any(~np.isnan(verr_3by3)):
verr_avg = np.nanmean(verr_3by3)
else:
verr_avg = np.nan
sampled2.append([v_avg, verr_avg])
sampled = np.array(sampled)
sampled2 = np.array(sampled2)
# print(sampled)
# print(row.Vx, float(df.loc[idx, 'SAV-peak-x']), float(df.loc[idx, 'SAV-peak-y']), sampled)
df.loc[idx, 'pt0_vxavg'] = sampled[0, 1]
df.loc[idx, 'pt0_vxgps'] = np.abs(gps.loc[0, 'vx (m/yr)'])
df.loc[idx, 'pt0_vyavg'] = sampled[0, 3]
df.loc[idx, 'pt0_vygps'] = np.abs(gps.loc[0, 'vy (m/yr)'])
df.loc[idx, 'pt1_vxavg'] = sampled[1, 1]
df.loc[idx, 'pt1_vxgps'] = np.abs(gps.loc[1, 'vx (m/yr)'])
df.loc[idx, 'pt1_vyavg'] = sampled[1, 3]
df.loc[idx, 'pt1_vygps'] = np.abs(gps.loc[1, 'vy (m/yr)'])
df.loc[idx, 'pt2_vxavg'] = sampled[2, 1]
df.loc[idx, 'pt2_vxgps'] = np.abs(gps.loc[2, 'vx (m/yr)'])
df.loc[idx, 'pt2_vyavg'] = sampled[2, 3]
df.loc[idx, 'pt2_vygps'] = np.abs(gps.loc[2, 'vy (m/yr)'])
df.loc[idx, 'pt0_vxdiff'] = sampled[0, 0] - gps.loc[0, 'vx (m/yr)']
df.loc[idx, 'pt0_vxavgdiff'] = sampled[0, 1] - gps.loc[0, 'vx (m/yr)']
df.loc[idx, 'pt0_vydiff'] = sampled[0, 2] - gps.loc[0, 'vy (m/yr)']
df.loc[idx, 'pt0_vyavgdiff'] = sampled[0, 3] - gps.loc[0, 'vy (m/yr)']
df.loc[idx, 'pt1_vxdiff'] = sampled[1, 0] - gps.loc[1, 'vx (m/yr)']
df.loc[idx, 'pt1_vxavgdiff'] = sampled[1, 1] - gps.loc[1, 'vx (m/yr)']
df.loc[idx, 'pt1_vydiff'] = sampled[1, 2] - gps.loc[1, 'vy (m/yr)']
df.loc[idx, 'pt1_vyavgdiff'] = sampled[1, 3] - gps.loc[1, 'vy (m/yr)']
df.loc[idx, 'pt2_vxdiff'] = sampled[2, 0] - gps.loc[2, 'vx (m/yr)']
df.loc[idx, 'pt2_vxavgdiff'] = sampled[2, 1] - gps.loc[2, 'vx (m/yr)']
df.loc[idx, 'pt2_vydiff'] = sampled[2, 2] - gps.loc[2, 'vy (m/yr)']
df.loc[idx, 'pt2_vyavgdiff'] = sampled[2, 3] - gps.loc[2, 'vy (m/yr)']
df.loc[idx, 'pt0_vavg'] = sampled2[0, 0]
df.loc[idx, 'pt0_vgps'] = gps.loc[0, 'v (m/yr)']
df.loc[idx, 'pt0_vdiff'] = sampled2[0, 0] - gps.loc[0, 'v (m/yr)']
df.loc[idx, 'pt0_verr'] = sampled2[0, 1]
df.loc[idx, 'pt1_vavg'] = sampled2[1, 0]
df.loc[idx, 'pt1_vgps'] = gps.loc[1, 'v (m/yr)']
df.loc[idx, 'pt1_vdiff'] = sampled2[1, 0] - gps.loc[1, 'v (m/yr)']
df.loc[idx, 'pt1_verr'] = sampled2[1, 1]
df.loc[idx, 'pt2_vavg'] = sampled2[2, 0]
df.loc[idx, 'pt2_vgps'] = gps.loc[2, 'v (m/yr)']
df.loc[idx, 'pt2_vdiff'] = sampled2[2, 0] - gps.loc[2, 'v (m/yr)']
df.loc[idx, 'pt2_verr'] = sampled2[2, 1]
You can comment/uncomment these lines to examine the data/results.
gps
# gps_xy
# df
| Unnamed: 0 | date1 | date2 | start_easting | start_northing | end_easting | end_northing | distance_traveled (m) | velocity (m/d) | start_easting_3413 | start_northing_3413 | end_easting_3413 | end_northing_3413 | geometry | vx (m/yr) | vy (m/yr) | v (m/yr) | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 2018-10-05 | 2018-09-30 | 621383.084841 | 6.738920e+06 | 621381.738315 | 6.738918e+06 | 2.015470 | -0.403094 | -3.227459e+06 | 212767.842312 | -3.227460e+06 | 212769.385111 | POINT (-3227460.401 212769.385) | 102.707969 | -112.624285 | 147.129311 |
| 1 | 1 | 2018-10-05 | 2018-09-30 | 610531.630118 | 6.737073e+06 | 610530.195434 | 6.737074e+06 | 1.703303 | -0.340661 | -3.228240e+06 | 224145.144740 | -3.228239e+06 | 224146.528571 | POINT (-3228239.026 224146.529) | -79.947760 | -101.019600 | 124.341145 |
| 2 | 2 | 2018-10-05 | 2018-09-30 | 601810.429736 | 6.733774e+06 | 601809.124834 | 6.733773e+06 | 1.881326 | -0.376265 | -3.230736e+06 | 233478.821992 | -3.230737e+06 | 233480.308331 | POINT (-3230737.265 233480.308) | 92.088135 | -108.502758 | 137.336804 |
df = df.replace(-np.inf, np.nan)
df = df.replace(np.inf, np.nan)
df.to_csv('../results_ITSLIVE.csv', index=False)