<jupyter_start><jupyter_text>Explore Grand Mesa data Calculating the Rolling Window Standard Deviation of the Elevation, and Exploratory Statistics of this Output<jupyter_code># Import necessary packages
import os
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
# Use geopandas for vector data and rasterio for raster data
import geopandas as gpd
import rasterio as rio
import rioxarray as riox
import xarray as xr
# Plotting extent is used to plot raster & vector data together
from rasterio.plot import plotting_extent
from scipy import stats

# Prettier plotting with seaborn
sns.set(font_scale=1.5, style="white")

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Set the Working Directory and Point to the Data<jupyter_code>#set working directory?
os.chdir(os.path.join('/Users',
                      'Hannah',
                    'Documents', 
                      'Documents_Grad', 
                      'Research', 
                      'ICESat_2',
                     'data'))
fn = os.path.join("ASO",
                        "snow_off",
                        "GrandMesa",
                       "ASO_3M_PCDTM_USCOGM_20160926_20160926.tif")
!pwd<jupyter_output>/Users/Hannah/Documents/Documents_Grad/Research/ICESat_2/data<jupyter_text>Open the Dataset using Xarray, Set -9999 values to nan, and Plot<jupyter_code>im = xr.open_rasterio(fn)
im = im.squeeze()
im = im.where(im != (-9999.))
#gm_std_low = gm_std_2d.where(gm_std_2d < 4)
plt.imshow(im)
plt.colorbar()
plt.axis('off');<jupyter_output><empty_output><jupyter_text>Introduce the Function I wrote, changed slightly to avoid creating a tuple for simplicity<jupyter_code># This function computes the rolling stadard deviation, taking as input an xarray DataArray

def rolling_stdev(im, size):                # im: a DataArray loaded in using xarray. size: define the size (dimension of one side) of your moving window
  """
  Short Summary
  --------------------
  This function calculates the rolling standard deviation within a window of
  specified dimensions.

  Extended Summary
  ------------------
  This function utilizes xarray's rolling window operation to first create a 
  Rolling object, and then to apply that Rolling object to the DataArray
  to calculate standard deviation within each (size, size) window. The output is
  a new DataArray with each pixel populated with the standard deviation of all 
  pixels within the window for which the pixel was at the center. This function
  also outputs a plot of the output file for quick and easy data visualization.

  Parameters
  ------------
  im:   DataArray
        Generally a .tif file that has been loaded in using xarray.open_rasterio()
  size 
        The dimension of the square rolling window to use. If one dimension is supplied, 
        it will be propagated to the proper number of dimensions. For example, a size
        set to 3 will generate a rolling window square of dimensions (3,3).

  Returns
  ------------
  tuple
        consisting of the new DataArray of standard deviation values
        and the figure and axes of the plot of this new DataArray
  
  Example
  -------------
  Implement by saving the output as a variable. 
  >>>output = rolling_stdev(im=im_subset, size=3)
  >>>output[0]
  xarray.DataArray

    band: 1y: 511x: 511

    nan nan nan nan nan ... 0.06020909 0.1668852 0.025760788 0.05106023
    Coordinates:
        band
        (band)
        int64
        1
        y
        (y)
        float64
        4.2e+06 4.2e+06 ... 4.198e+06
        x
        (x)
        float64
        2.272e+05 2.273e+05 ... 2.288e+05
    Attributes: (0)

  """
  
  r = im.rolling(y=size)                    # this creates the xarray Rolling object 
  im_stdev = r.reduce(np.std)               # use the Rolling object to calculate the standard deviation within each window. Could substitute other functions like mean or median. 
  
  return im_stdev<jupyter_output><empty_output><jupyter_text>Implement the Function<jupyter_code>gm_std = rolling_stdev(im, 3)
gm_std<jupyter_output><empty_output><jupyter_text>Plot the Standard Deviation output, and Exclude Outliers<jupyter_code>fig, ax = plt.subplots(ncols=2, figsize=(20,10))

c1 = ax[0].imshow(gm_std)
c2 = ax[1].imshow(gm_std, vmin=0, vmax=5)
plt.colorbar(c1, ax=ax[0])
plt.colorbar(c2, ax=ax[1]);
gm_std.shape<jupyter_output><empty_output><jupyter_text>Output Statistics About the Standard Deviation Output, potentially for use as thresholding values outside of Grand Mesa<jupyter_code>print(np.nanmedian(gm_std))
print(np.nanmean(gm_std))
print(np.nanmax(gm_std))
print(np.nanmin(gm_std))
#print(stats.mode(gm_std.values))<jupyter_output>0.18067864
0.34063563
66.755646
0.0<jupyter_text>Subset to Standard Deviations Below 5, to Remove Outliers<jupyter_code>gm_std_low = gm_std.where(gm_std < 5)
plt.imshow(gm_std_low)
plt.colorbar();
print(np.nanmean(gm_std_low.values))

plt.hist(np.ravel(gm_std_low), bins=100);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09242364141213805, 'fraction_numerical': 0.03351844506148354, 'mean_word_length': 3.209515859766277, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '555818', 'license': 'MIT', 'n_tokens_mistral': 1757, 'n_tokens_neox': 1573, 'n_words': 503, 'path': 'FirstLookGM.ipynb', 'repo_name': 'bessoh2/IA_DEM_Roughness', 'size': '5042'}