<jupyter_start><jupyter_text>Load and interpolate data sources<jupyter_code>data_dir = '/home/teisberg/data/common_data/'

nc_measure = Dataset(os.path.join(data_dir, 'antarctica_ice_velocity_450m_v2.nc'))
nc_bedmachine = Dataset(os.path.join(data_dir, 'BedMachineAntarctica_2019-11-05_v01.nc'))
crs_latlon = pyproj.crs.CRS.from_proj4("+proj=latlon")
crs_epsg3031 = pyproj.crs.CRS.from_epsg(3031)
transformer_latlon_to_3031 = pyproj.Transformer.from_crs(crs_latlon, crs_epsg3031)
start_point = [550e3, -0.8e6] # Byrd
X_ctr, Y_ctr = start_point
size_x, size_y = 200e3, 200e3
x = np.array(nc_bedmachine.variables['x'])
y = np.array(nc_bedmachine.variables['y'])
mask_x = (x > (X_ctr - size_x)) & (x < (X_ctr + size_x))
mask_y = (y > (Y_ctr - size_y)) & (y < (Y_ctr + size_y))
bm_X, bm_Y = np.meshgrid(x[mask_x], y[mask_y])
bm_h = np.array(nc_bedmachine.variables['thickness'])[mask_y, :][:, mask_x]
bm_source = np.array(nc_bedmachine.variables['source'])[mask_y, :][:, mask_x]
bm_mask = np.array(nc_bedmachine.variables['mask'])[mask_y, :][:, mask_x]
bm_surface = np.array(nc_bedmachine.variables['surface'])[mask_y, :][:, mask_x]

x = np.array(nc_measure.variables['x'])
y = np.array(nc_measure.variables['y'])
mask_x = (x > (X_ctr - size_x)) & (x < (X_ctr + size_x))
mask_y = (y > (Y_ctr - size_y)) & (y < (Y_ctr + size_y))
measure_X, measure_Y = np.meshgrid(x[mask_x], y[mask_y])
measure_vx = np.array(nc_measure.variables['VX'])[mask_y, :][:, mask_x]
measure_vy = np.array(nc_measure.variables['VY'])[mask_y, :][:, mask_x]
# Interpolate everything to the BedMachine 450 m grid
bm_vx = scipy.interpolate.griddata((measure_X.flatten(), measure_Y.flatten()), measure_vx.flatten(),
                                   (bm_X, bm_Y))
bm_vy = scipy.interpolate.griddata((measure_X.flatten(), measure_Y.flatten()), measure_vy.flatten(),
                                   (bm_X, bm_Y))

vx_interp = scipy.interpolate.LinearNDInterpolator((measure_X.flatten(), measure_Y.flatten()), measure_vx.flatten())
vy_interp = scipy.interpolate.LinearNDInterpolator((measure_X.flatten(), measure_Y.flatten()), measure_vy.flatten())
h_interp = scipy.interpolate.LinearNDInterpolator((bm_X.flatten(), bm_Y.flatten()), bm_h.flatten())
# Load cresis data

cresis_2017_basler = pd.read_csv(os.path.join(data_dir, "cresis/Browse_2017_Antarctica_Basler.csv")).rename(columns = {'UTCTIMESOD': 'TIME'}) # TODO: See note below!
# The name of the time column changed at some point. We don't use the time column, so I'm just renaming it to avoid it being NaN later. If you want to use the time for anything, check on this!
cresis_2011_dc8 = pd.read_csv(os.path.join(data_dir, "cresis/Browse_2011_Antarctica_DC8.csv"))
cresis_2011_to = pd.read_csv(os.path.join(data_dir, "cresis/2011_Antarctica_TO.csv"))
cresis_df = pd.concat([cresis_2011_dc8, cresis_2011_to, cresis_2017_basler])
cresis_df = cresis_df[['LAT', 'LON', 'THICK']]
cresis_df['SOURCE'] = 'cresis'
# Load HiCARS data

headers = "YEAR DOY SOD LON LAT THICK SRF_RNG BED_ELEVATION SURFACE_ELEVATION PARTIAL_BED_REFLECT SRF_RELFECT AIRCRAFT_ROLL".split(' ')
# Note I renamed THK to THICK for consistency with cresis data

all_dfs = []

for f in glob.iglob(data_dir + r'hicars-byrd/**/*icethk.txt', recursive=True):
    df = pd.read_csv(f, sep=' ', comment='#', names=headers, index_col=False)
    df = df[['LON', 'LAT', 'THICK']].dropna()
    all_dfs.append(df)

hicars_df = pd.concat(all_dfs).reset_index(drop=True)
hicars_df['SOURCE'] = 'hicars'
radar_df = pd.concat([cresis_df, hicars_df]).reset_index(drop=True)

radar_proj = [transformer_latlon_to_3031.transform(lon, lat) for lat,lon in zip(radar_df['LAT'], radar_df['LON'])]
radar_df['X'] = [a[0] for a in radar_proj]
radar_df['Y'] = [a[1] for a in radar_proj]


radar_df['BM_DIFF'] = [h_interp((x,y))-h for x,y,h in zip(radar_df['X'], radar_df['Y'], radar_df['THICK'])]

radar_df = radar_df.dropna()
radar_df = radar_df[radar_df['THICK'] >= 0] # No negative ice thicknesses plz

radar_kdtree = scipy.spatial.KDTree(np.array(list(zip(radar_df['X'], radar_df['Y']))))
fig, axs = plt.subplots(2, 3, figsize=(18,8), facecolor='white', sharex=True, sharey=True)

h_norm = matplotlib.colors.Normalize(vmin=0, vmax=np.max(bm_h))

pcm0 = axs[0,0].pcolormesh(bm_X, bm_Y, bm_h, norm=h_norm, shading='nearest')
clb0 = fig.colorbar(pcm0, ax=axs[0,0])
clb0.set_label('BedMachine Ice Thickness [m]')
axs[0,0].set_aspect('equal')

pcm1 = axs[0,1].pcolormesh(measure_X, measure_Y, np.sqrt(measure_vx**2 + measure_vy**2), shading='nearest')
arrow_subsample = 20
axs[0,1].quiver(measure_X[::arrow_subsample,::arrow_subsample],
              measure_Y[::arrow_subsample,::arrow_subsample],
              measure_vx[::arrow_subsample,::arrow_subsample],
              measure_vy[::arrow_subsample,::arrow_subsample],
              angles='xy', color='white', scale=5e4)
clb1 = fig.colorbar(pcm1, ax=axs[0,1])
clb1.set_label('MEaSURE Ice Velocity [m/yr]')
axs[0,1].set_aspect('equal')

bedmachine_source_key = {k: v for k, v in zip(nc_bedmachine.variables['source'].flag_values,
                                              nc_bedmachine.variables['source'].flag_meanings.split())}
pcm2 = axs[1,0].pcolormesh(bm_X, bm_Y, bm_source, vmin=1, vmax=10, cmap=plt.cm.get_cmap('terrain', 11), shading='nearest')
formatter2 = plt.FuncFormatter(lambda val, loc: bedmachine_source_key.get(val, ""))
clb2 = fig.colorbar(pcm2, ax=axs[1,0], format=formatter2)
clb2.set_label('BedMachine Ice Thickness Source')
axs[1,0].set_aspect('equal')

bedmachine_mask_key = {k: v for k, v in zip(nc_bedmachine.variables['mask'].flag_values,
                                              nc_bedmachine.variables['mask'].flag_meanings.split())}
pcm3 = axs[1,1].pcolormesh(bm_X, bm_Y, bm_mask, vmin=1, vmax=4, cmap=plt.cm.get_cmap('terrain', 11), shading='nearest')
formatter3 = plt.FuncFormatter(lambda val, loc: bedmachine_mask_key.get(val, ""))
clb3 = fig.colorbar(pcm3, ax=axs[1,1], format=formatter3)
clb3.set_label('BedMachine Mask')
axs[1,1].set_aspect('equal')

pcm4 = axs[0,2].scatter(radar_df['X'], radar_df['Y'], c=radar_df['THICK'], norm=h_norm, s=0.1)
clb4 = fig.colorbar(pcm4, ax=axs[0,2])
clb4.set_label('Radar Ice Thickness [m]')
axs[0,2].set_aspect('equal')

sc = axs[1,2].scatter(radar_df['X'], radar_df['Y'], c=radar_df['BM_DIFF'], s=0.1, cmap=matplotlib.cm.get_cmap('coolwarm'), vmin=-200, vmax=200)
clb5 = fig.colorbar(sc, ax=axs[1,2])
axs[1,2].set_aspect('equal')
clb5.set_label('Bed Machine Thickness - Radar Thickness')

axs[1,0].tick_params(axis='x', labelrotation=45) 
axs[1,1].tick_params(axis='x', labelrotation=45)
axs[1,2].tick_params(axis='x', labelrotation=45)

fig.tight_layout()

plt.show()
ice_free_mask = bm_mask == 1

data_dict = {
        'radar': {
            'x': np.concatenate([np.array(radar_df['X']), bm_X[ice_free_mask]]),
            'y': np.concatenate([np.array(radar_df['Y']), bm_Y[ice_free_mask]]),
            'h': np.concatenate([np.array(radar_df['THICK']), 0*bm_Y[ice_free_mask]])
        },
        'bedmachine': {
            'x': bm_X,
            'y': bm_Y,
            'h': bm_h
        },
        'velocity': {
            'x': measure_X,
            'y': measure_Y,
            'vx': measure_vx,
            'vy': measure_vy
        }
    }

# Add in interpolated velocity to radar data
    
tree = scipy.spatial.KDTree(list(zip(measure_X.flatten(), measure_Y.flatten())))
    
vals_out = []
for x, y in zip(data_dict['radar']['x'], data_dict['radar']['y']):
    dist, idx = tree.query((x,y), k=1)
    vals_out.append(np.sqrt(measure_vx.flatten()[idx]**2 + measure_vy.flatten()[idx]**2))

data_dict['radar']['v_nn'] = np.array(vals_out)
with open('byrd-data.pickle', 'wb') as f:
    pickle.dump(data_dict, f)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17590797466718366, 'fraction_numerical': 0.026754556029468787, 'mean_word_length': 5.1025236593059935, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '242880', 'license': 'MIT', 'n_tokens_mistral': 3346, 'n_tokens_neox': 3197, 'n_words': 495, 'path': 'data_preprocess.ipynb', 'repo_name': 'thomasteisberg/igarss2021', 'size': '7737'}