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:

  1. Get and print the UTM coordinates of three GNSS statations for each scene pair.

  2. Sample every velocity maps.

  3. 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)