Snow Depth Applications

import numpy as np
import pandas as pd
import geopandas as gpd
import rioxarray as rxr
import matplotlib.pyplot as plt
from sliderule import icesat2, sliderule
from scipy.interpolate import RectBivariateSpline

icesat2.init("slideruleearth.io")

# Define region of interest over Toolik, Alaska
region = [{"lon":-149.5992624418217, "lat":68.63358948385529}, 
          {"lon":-149.5954459662985, "lat":68.60200878043223}, 
          {"lon":-149.2821268688734, "lat":68.60675802967609}, 
          {"lon":-149.2855031235162, "lat":68.63834638180673},
          {"lon":-149.5992624418217, "lat":68.63358948385529}]

print(region)

## Build SlideRule request
# Define parameters (described below)
parms = {
    "poly": region,
    "srt": icesat2.SRT_LAND,
    "cnf": icesat2.CNF_SURFACE_HIGH,
    "atl08_class": ["atl08_ground"],
    "ats": 5.0,
    "len": 20.0,
    "res": 10.0,
    "maxi": 5
}

# Calculated ATL06 dataframe
is2_df = icesat2.atl06p(parms)

# Print SlideRule output
print(is2_df.head())

import contextily as ctx
from shapely.geometry import Polygon

# Convert region to a Polygon
coords = [(point["lon"], point["lat"]) for point in region]
polygon = Polygon(coords)
region_gdf = gpd.GeoDataFrame([1], geometry=[polygon], crs="EPSG:4326")

# Reproject to Web Mercator for contextily
is2_df_mercator = is2_df.to_crs(epsg=3857)
region_mercator = region_gdf.to_crs(epsg=3857)


fig, ax = plt.subplots(figsize=(12, 8))
# Plot surface height
is2_df_mercator.plot(column='h_mean', 
                  ax=ax, 
                  cmap='viridis',
                  legend=True,
                  markersize=10,
                  alpha=0.8)

# Plot the region bounding box
region_mercator.plot(ax=ax, 
                     facecolor='none', 
                     edgecolor='red', 
                     linewidth=2)

# Add ESRI World Imagery basemap
ctx.add_basemap(ax, 
                crs=is2_df_mercator.crs, 
                source=ctx.providers.Esri.WorldImagery)
plt.tight_layout()
plt.show()

# Subset ICESat-2 data to single RGT, time of year
is2_df_subset = is2_df[is2_df['rgt']==152]
is2_df_subset = is2_df_subset.loc['2023-03-31']

# Display top of dataframe
print(is2_df_subset.head())

import earthaccess
earthaccess.login(strategy='interactive', persist=True)
auth = earthaccess.login()

# Coordinates for SW/NE corners
lon_min = min([coord[0] for coord in coords])
lat_min = min([coord[1] for coord in coords])
lon_max = max([coord[0] for coord in coords])
lat_max = max([coord[1] for coord in coords])

# Data query for lidar point clouds over Fairbanks, AK
results = earthaccess.search_data(
    short_name='SNEX23_Lidar',
    bounding_box = (lon_min, lat_min, lon_max, lat_max),
    temporal = ('2023-03-13', '2023-03-14')
)

# File paths for UAF rasters (TODO: Add lidar files, and update names)
tifpath = "/your/path/here"
f_snow_off = f"{tifpath}/uaf_lidar_snowoff.tif"
f_snow_on = f"{tifpath}/uaf_lidar_snowon.tif"

# Load files as rioxarray datasets
lidar_snow_off = rxr.open_rasterio(f_snow_off)
lidar_snow_on = rxr.open_rasterio(f_snow_on)

# Initialize ICESat-2 coordinate projection
is2_df_subset = is2_df_subset.set_crs("EPSG:4326")

# Change to EPSG:32606
is2_df_subset = is2_df_subset.to_crs("EPSG:32606")

# Display top of dataframe
print(is2_df_subset.head())

# Make coregistration function
def coregister_is2(lidar_snow_off, lidar_snow_on, is2_df):
    """
    Co-registers UAF data with ICESat-2 data with a rectangular 
    bivariate spline.

    Parameters
    ------------
    lidar_snow_off: rioxarray dataset
        Lidar DEM/DTM in rioxarray format.
    lidar_snow_on: rioxarray dataset
        Lidar-derived snow depth in rioxarray format.
    is2_df: GeodataFrame
        GeoDataFrame for the ICESat-2 data generated with SlideRule.

    Returns
    ------------
    is2_uaf_df: GeoDataFrame
        Contains the coordinate and elevation data that matches best
        between ICESat-2 and UAF.
    """

    # Helper function to prepare lidar data
    def prepare_lidar_data(raster):
        coords_x = np.array(raster.x)
        coords_y = np.array(raster.y)
        values = np.array(raster.sel(band=1))[::-1, :]
        values[np.isnan(values)] = -9999
        return coords_x, coords_y, values

    # Get coordinates and height/depth values from lidar data
    x0, y0, dem_heights = prepare_lidar_data(lidar_snow_off)
    xs, ys, dem_depths = prepare_lidar_data(lidar_snow_on)

    # Generate interpolators
    interp_height = RectBivariateSpline(y0[::-1], x0, dem_heights)
    interp_height = RectBivariateSpline(ys[::-1], xs, dem_depths)

    # Pre-filter IS2 data to bounds (apply once instead of per beam)
    x_bounds = (is2_df.geometry.x > np.min(x0)) & (is2_df.geometry.x < np.max(x0))
    y_bounds = (is2_df.geometry.y > np.min(y0)) & (is2_df.geometry.y < np.max(y0))
    is2_filtered = is2_df[x_bounds & y_bounds].copy()

    if is2_filtered.empty:
        print('Error with GeoDataFrame or raster bounds.')
        return gpd.GeoDataFrame()

    # Extract coordinates once
    xn = is2_filtered.geometry.x.values
    yn = is2_filtered.geometry.y.values

    # Estimate lidar height and snow depth at the ICESat-2 coordinates
    lidar_heights = interp_height(yn, xn, grid=False)
    lidar_snow_depths = interp_depth(yn, xn, grid=False)

    # Create result DataFrame in one operation
    is2_uaf_df = gpd.GeoDataFrame({
        'x': xn,
        'y': yn,
        'time': is2_filtered.index.values,
        'beam': is2_filtered['gt'].values,
        'lidar_height': lidar_heights,
        'lidar_snow_depth': lidar_snow_depths,
        'is2_height': is2_filtered['h_mean'].values,
        'h_sigma': is2_filtered['h_sigma'].values,
        'dh_fit_dx': is2_filtered['dh_fit_dx'].values
    })

    # Add coordinate transformation
    transformer = Transformer.from_crs("EPSG:32606", "EPSG:4326", always_xy=True)
    is2_uaf_df['lon'], is2_uaf_df['lat'] = transformer.transform(
        is2_uaf_df['x'], is2_uaf_df['y']
    )

    return is2_uaf_df

# Co-locate ICESat-2 and UAF using the above function
is2_uaf_df = coregister_is2(lidar_snow_off, 
                            lidar_snow_on, 
                            is2_df_subset
                           )
# Convert to a GeoDataFrame
geom = gpd.points_from_xy(is2_uaf_df.lon, is2_uaf_df.lat)
is2_uaf_gdf = gpd.GeoDataFrame(is2_uaf_df,
                               geometry=geom,
                               crs="EPSG:4326"
                              )
# Print head of GeoDataFrame
is2_uaf_gdf.head()

# Derive snow depth using snow-on/snow-off differencing
is2_uaf_gdf['is2_snow_depth'] = is2_uaf_gdf['is2_height'] - is2_uaf_gdf['lidar_height']

# Estimate the residual (bias) between IS-2 and UAF depths
is2_uaf_gdf['snow_depth_residual'] = is2_uaf_gdf['is2_snow_depth'] - is2_uaf_gdf['lidar_snow_depth']

# Create plot
fig, ax = plt.subplots(figsize=(12, 8))
# Plot surface height
is2_uaf_gdf.to_crs("EPSG:3857").plot(column='is2_snow_depth', 
                  ax=ax, 
                  cmap='viridis',
                  legend=True,
                  markersize=10,
                  alpha=0.8)

# Plot the region bounding box
region_mercator.plot(ax=ax, 
                     facecolor='none', 
                     edgecolor='red', 
                     linewidth=2)

# Add ESRI World Imagery basemap
ctx.add_basemap(ax, 
                crs="EPSG:3857", 
                source=ctx.providers.Esri.WorldImagery)
plt.tight_layout()
plt.show()

# Plot snow depths for the three strong beams
fig,axs = plt.subplots(3)

#Left strong beam
tmp_df = is2_uaf_gdf[is2_uaf_gdf['beam']==10]
axs[0].plot(tmp_df['lat'], tmp_df['is2_snow_depth'], label='ICESat-2')
axs[0].plot(tmp_df['lat'], tmp_df['lidar_snow_depth'], label='UAF')
axs[0].set_title('gt1l')
axs[0].legend()

# Central strong beam
tmp_df = is2_uaf_gdf[is2_uaf_gdf['beam']==30]
axs[1].plot(tmp_df['lat'], tmp_df['is2_snow_depth'])
axs[1].plot(tmp_df['lat'], tmp_df['lidar_snow_depth'])
axs[1].set_ylabel('Snow depth [m]', fontsize=18)
axs[1].set_title('gt2l')

# Right strong beam
tmp_df = is2_uaf_gdf[is2_uaf_gdf['beam']==50]
axs[2].plot(tmp_df['lat'], tmp_df['is2_snow_depth'])
axs[2].plot(tmp_df['lat'], tmp_df['lidar_snow_depth'])
axs[2].set_xlabel('Latitude [m]', fontsize=18)
axs[2].set_title('gt3l')
plt.tight_layout()

# Only include outer axis labels
for ax in axs:
    ax.label_outer()
    ax.set_ylim([0, 1.5])

fig, ax = plt.subplots(figsize=(12,6))
s = is2_uaf_gpd.plot.scatter(ax=ax,
                             x='lidar_snow_depth',
                             y='is2_snow_depth',
                             c='snow_depth_residual',
                             vmin=-0.5, vmax=0.5
                            )
ax.set_xlabel("UAF snow depth [m]", fontsize=18)
ax.set_ylabel("ICESat-2 snow depth [m]", fontsize=18)
ax.set_xlim([0, 1.5])
ax.set_ylim([0, 1.5])
cbar = fig.colorbar(s, ax=ax)
cbar.set_label("Snow depth residual [m]")
plt.tight_layout()