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)