# User Manual for Slush Limit Detection in Greenland
## 1 Introduction
The programs in this repository are used to find the slush limit in MODIS scenes.
The approach is described in *Machguth et al.*, Journal of Glaciology (2022) and is partly based
on *Greuell and Knap* (2000) who originally developed their algorithm for
application to AVHRR scenes.
The approach is split in four steps that are performed successively. The steps
are summarized briefly in the first section of this manual and explained in
detail in the third and fourth section. The second section provides basic
information about the required Python setup to run the code.
## 2 Brief Summary of the Projects
The code for the approach is split into the following four projects:
1. get_modis
2. modis_tools
3. MODIS_Greenland
4. MODIS_Greenland_analysis
The first downloads the files, the second pre-processes the data, the
third performs the actual computation of the slush limits and the fourth
step is about analysing the files. Output is created at the end of
steps 1, 2, 3 and 4. The output is always the input for the subsequent step.
This step-wise procedure was initially chosen for development but was then kept.
An overview of the inputs and outputs for all the steps is given in Table 1.
*Table 1:* Overview of the processing.
| Step | Code | Input | Output |
|------|--------------------------------------|--------------------------------------|-------------------------------------------------------|
| 1 | `get_modis_batch.py` | NASA EarthData MOD10A1 (and MOD09GA) | MODIS data |
| 2 | `stitch_warp_modis_l2.py` | MODIS data & configuration template | reprojects, stitched and cropped GeoTIFFs |
| 3a | `MODIS_NDWI.py` | processed MOD09GA from step 2 | NDWI grids (.nc) |
| 3b | `MODIS_array_filter_parallel.py` | processed MOD10A1 from step 2 | filtered grids (.nc) |
| 3c | `MODIS_stddev_spatial_parallel.py` | grids from step 3b | albedo and stddev. albedo grids (.nc) |
| 3d | `MODIS_find_slush_limit_parellel.py` | grids from step 3c and DEM grid | Excel and .csv sheets Ys candidates, plots |
| 4a | `MODIS_A_SL_max.py` | Excel sheet from step 3d | Maximum annual slush limits (Excel) and various plots |
| 4b | `MODIS_A_SLanalysis.py` | Excel sheets from step 4a | Maximum annual slush limits (Excel) and various plots |
### 2.1 get_modis
The code downloads MODIS MOD10A1 and MOD09GA images in the form of *.hdf* files
by using `get_modis_batch.py`. ice
### 2.2 modis_tools
In this project the files are stitched together and warped to a specific
projection using `stitch_warp_modis_l2.py`. The project also includes templates
to specify the input. The output are files in the form of GeoTiff.
Alternatively, for example Google Earth Engine can be used for download of MODIS
MOD10A1 and MOD09GA data (making get_modis and modis_tools obsolete).
### 2.3 MODIS_Greenland
In this project the files are processed in four steps. (1) The MODIS files
are filtered according to the steps explained by *Machguth et al.* (2022) by using
`MODIS_array_filter_parallel.py`. (2) NDWI files are created using `MODIS_NDWI.py`.
The NDWI files can also be acquired using Google Earth Engine. These can be
converted from .tiff to .nc files using `MODIS_NDWI_tiff_to_nc.py`. (3) The
spatial standard deviation of surface albedo is calculated using
`MODIS_stddev_spatial_parallel.py`. (4) The slush limit is detected by running
the code `MODIS_find_slush_limit_parallel.py`.
This fourth step creates the following outputs: (i) an Excel sheet with the
metadata of all slush limit candidates. (ii) plots for each individual year
containing all the detected slush limits over all the elevation bands.
(iii) plots for each individual year and each individual elevation band
containing all the detected slush limits (see Fig. 2 in *Machguth et al.,* 2022).
These plots also show detected outliers among the slush limits candidates.
(iv) Folders for each elevation band containing daily plots for all days with
slush limit candidates (see Fig. 1 in *Machguth et al.,* 2022).
### 2.4 MODIS_Greenland_analysis
In this project the files are analysed by using `MODIS_A_SLmax.py`. This
program creates a number of plots and Excel sheets that help analyse the
results generated in the previous project.
The program `MODIS_A_SLanalysis.py` creates a number of plots (such as
Fig. 6 in *Machguth et al.,* 2022) and also calculates PDH statistics for the
K-Transect as shown in *Machguth et al.,* (2022). Note that in the code of
`MODIS_A_SLanalysis.py` one can choose whether only the plots of Ys progression
(Fig. 6 in *Machguth et al.,* 2022), only the analysis of the forcing of progression
of Ys (K-Transect analysis) or both should be done.
## 3 Preparations for Using the Code
The code is in Python. In the given example below the code was executed inside
PyCharm. Packages used in the code were downloaded from conda-forge and
installed using Anaconda. The code can be either forked from http://github.com/machguth
or downloaded here. Note that the newest versions on Github differ substantially
from the code here. The following steps explain the setup.
### 3.1 Packages required
To be able to do all the steps for the slush limit processing and analysis
install the following packages: *bs4, numpy, pandas, rasterio, xarray, joblib,
gdal, netcdf4, matplotlib, geopandas, shapely, scipy, pyproj,
lxml, openpyxl.*
## 4 Calculation and Analysis of the Slush Limits
The following steps explain the download, the processing and the analysis of the
data. Modify the variables mentioned in the following according to your area of
interest. Consider Table 1 for an overview of input and output for the programs.
Note that the steps 2 to 4 are carried out in parallel computing, trying to
make optimal use of the available CPUs (more precisely of the available threads)
on a computer. The number of available threads needs to be specified in the
settings of the code.
### 4.1 Raw Data: `get_modis_batch.py`
Below explained is the download using `get_modis_batch.py`. Below only basic functionality
is explained, more details on the functioning and settings can be found in the code of
`get_modis.py` and in the readme of the *get_modis* project.
> The MODIS data can also be downloaded using other methods such as for example Google Earth Engine.
1. For the download of MODIS data an account at the NASA EarthData website is necessary. The username and the password will be used in the code, therefore use a different password than elsewhere.
2. `yr_start`, `yr_end`: Period that will be downloaded, last year is included
3. `user`, `pw`: Enter your username and password
4. `outdir`: Modify the path: This is where your files will be saved at, i.e. your output files
5. `tiles`: Modify the tiles where your area of interest is located in
6. `codepath`: Enter the path to the `get_modis_batch.py` program
7. Run the program to start the download
**Output:** This step downloads the files from the NASA EarthData website.
The files are images from the MODIS sensor on the Terra Satellite. The output
are the raw files in the form of *.hdf* files.
### 4.2 Processing to Level 1: `stich_warp_modis_l2.py`
Below only basic functionality is explained, more details on the functioning and settings
can be found in the code of `stich_warp_modis_l2.py` and in the readme of
the *modis_tools* project.
> Any other tool can be used to stitch, reproject and crop the MODIS data. The output needs to be
in geoTIFF format.
**Input:** The input for the first processing step are MOD10A1 *.hdf* files
(see Section 4.1). Furthermore, you need to create a template for your area
of interest, best start with one of the existing templates, such as
`MODIS_reprojection_MOD10A1_EPSG3413.template`.
1. Modify `year_start`, `year_end`, `day_start` and `day_end` according to the files you want to download. This should be the same as in the previous Section 4.1.
2. `input_path`, `output_path`: Modify the paths
3. Define your area of interest in the template file.
4. The output folder where the files are saved (`output_path`) needs to be created manually before the program is run, otherwise there will be an error
5. Run `stich_warp_modis_l2.py` using the parameter file as input.
**Output:** The outputs are the satellite images stitched together, reprojected
and cropped to the study area. They are in GeoTiff format.
### 4.3 Processing to Level 2: `MODIS_array_filter_parallel.py`
Performs filtering to remove erroneous values as described in *Machguth et al.,* (2022).
**Input:** The input for the second processing step are the files form
Section 4.2.
1. `infolder`, `outfolder`: Modify the path
2. `n_cores`: Adjust to the number of CPU cores that you want to use for the processing. Here and in the following steps: You can specify the number of threads instead of physical cores (the latter often being twice as high as the former). If you enter a value that exceeds the actual number of threads, the code will reset the value to the actual maximum.
3. Run `MODIS_array_filter_parallel.py` to process the files
**Output:** The files are now filtered for non-albedo values, clouds and
erroneous albedo values. Output files are netCDF.
### 4.4 Processing to Level 3: `MODIS_stddev_spatial_parallel.py`
Calculates spatial standard deviation of surface albedo.
**Input:** The input for the third processing step are the files form Section 4.3.
1. `infolder`, `outfolder`: Modify the path
2. `n_cores`: Set number of CPU cores used
3. Run `MODIS_stddev_spatial_parallel.py` to process the files
**Output:** Daily netCDF files containing albedo and spatial standard deviation of albedo.
### 4.5 Calculating of the slush limits Ys: `MODIS_find_slush_limit_parallel.py`
**Input:** The input for the fourth processing step are the files form Section 4.4
and optionally the NDWI-files. In addition, a digital elevation model (DEM)
is required. The DEM has to fulfil the following conditions: (i) Format GeoTIFF,
(ii) same projection as the pre-processed MODIS grids, (ii) same extent and
spatial resolution as the MODIS grids. (iv) The DEM serves also as glacier or
ice sheet mask. This means in the DEM all grid cells outside the glacier or
ice sheet perimeters need to be set to NaN. Otherwise the algorithm will attempt
to detect slush limits also for non-glaciated areas.
1. `infolder`, `infolder_ndwi`, `outfolder`, `demfile`: Modify the paths or file names
2. `num_cores`: Set number of CPU cores used
3. Specify further options as explained below
4. Specify output as explained below
5. Run the program `MODIS_array_filter_parallel.py` for the final processing of the files
**Options:** The various options and parameter settings are explained in the
code. Key options are detailed here.
`batch_stripes`: If `batch_stripes`, then the domain is divided in latitudinal stripes that
run from the
very west to the very east of the domain and are of north-south extent as
specified in `ycw`. The algorithm automatically analyses how many of these
stripes can be fitted into the north-south extent of the domain. Scanning always
starts from the south. Note that batch_stripes are only suitable if the ice
sheet surface slopes approximately from east to west (or west to east).
`yc`: If not `batch_stripes`, then one single stripe is analysed
(running east-west or west-east), of north-south extent `ycw`. This option is
useful for fast model runs and debugging. Note that useful results are only
obtained if the ice sheet surface slopes approximately from east to west
(or west to east).
`NDWI`: If `False`, then no NDWI input grids are being read and the slush limit
candidates are identified solely based on spatial standard deviation of
surface albedo and surface albedo. This option can be used to more closely mimic
the original approach by *Greuell and Knap* (2000). Their approach was based on
AVHRR data that do not allow the calculation of NDWI (it lacks the blue channel
needed for the calculation of NDWIice, cf. *Yang and Smith,* 2013).
`q_index`: This option should typically set to `False` and while still functional,
is currently poorly maintained. The option is related to the filtering for
outliers and is obsolete when MODIS data are used and `NDWI` is `True`. If
`q_index == True`, then for each detected slush limit candidate a quality index
is calculated. The calculation is done in the function
`MODIS_find_slush_limit_funcs.qindex_SL`. The quality index attempts to
evaluate the reliability of each detected slush limit candidate and uses the
calculated quality index in the filtering for outliers (the basic assumption
being that candidates with a lower quality index are more likely outliers).
The `q_index` option is still included in the code as it could become relevant
again if ever attempting to use, e.g., AVHRR data.
`ziw`: Specifies the elevation extent of the elevation bins that are used in
detecting the slush limit candidates (see *Machguth et al.,* 2022).
Typically, `ziw` is set to 20 m. The parameter should be chosen in a way that the
number of grid cells falling into an elevation bin is big enough to calculate
reliable statistic parameters over each elevation bin. The number of pixels per
elevation bin is also controlled by the width of the stripes (width of the
polygons or parameter `ycw`), the resolution of the input grids and the steepness
of the terrain.
`start_year`, `end_year`: Specifies the range of years over which the algorithm
should be ran. `end_year` is included in the range.
**Output:** This step concludes the calculation of slush limits and creates a
number of outputs (see Table 1 and Section 2.3). Mandatory output is (i) an
Excel and *.csv* file that contains comprehensive information on every detected
slush limit candidate, and (ii) one plot per year and stripe which shows the
evolution of the detected slush limits over the course of the melt season.
Detected outliers are shown but marked in these plots. Additional
output can be requested as follows:
`daily_plots`: set to `True` to obtain a plot for every
detected slush limit candidate. This option will substantially slow down
processing and use disk space. Option is recommended to visualize and/or analyse
algorithm performance.
`all_daily_plots`: set to `True` to obtain a plot for each time
step where cloudiness is low enough to theoretically detect a slush limit
candidate. Option is recommended to analyse algorithm performance and for debugging. If
set to `True`, then the option `daily_plots` is automatically also selected. This
option will very substantially slow down processing and use disk space.
### 4.6 Analysis: `MODIS_A_SLmax.py`
Here annual maximum slush limits are calculated based on the list of valid slush limit
candidates calculated in Section 4.5.
**Input:** Excel sheet that was created in Section 4.5.
1. `data`, `outfolder`: Modify the file and path
2. Run `MODIS_A_SLmax.py` to create plots for analysing the results
**Output:** The output of the analysis creates two new Excel sheets with metadata
about the valid slush limit candidates as well as a simple overview sheet of all the
annual maximum slush limits over all the years. Furthermore, a number of plots are
created that help analyse the results generated in the processing steps.
**Options:** `Reeh` and `Greuell` specify two Excel tables that contain the data from earlier
estimates of the West Greenland slush limit (*Reeh,* 1991; *Greuell and Knap,* 2000).
The two tables are used to create figures that compare the slush limits derived with
the present algorithm to these earlier estimates. The two Excel sheets are contained in the
*Docs* folder of this project. `MODIS_A_SLmax.py` should work without the two tables, in
that case all the other output will still be written. To avoid error messages the user
can comment out the call to the functions `boxplot_max_slush_lim_comp_reeh`
and `boxplot_max_slush_lim_comp_reeh_vertical`.
## References
Greuell, W. B. and W. H. Knap (2000). Remote sensing of the albedo and detection of the
slush line on the Greenland ice sheet, *J. Geophys. Res.*, **105**(D2), 15567–15576
(doi: 10.1029/1999JD901162).
Machguth, H., Tedstone, A., Mattea, E. (2022). Daily Variations in Western
Greenland Slush Limits, 2000 to 2021. *Journal of Glaciology*.
Reeh, N. (1991): Parameterizations of melt rate and surface temperature on the
Greenland ice sheet, Polarforschung, **59,** 113-128.
Yang, K. and L. C. Smith (2013): Supraglacial streams on the Greenland Ice Sheet
delineated from combined spectral-shape information in high-resolution satellite
imagery, *IEEE Geosci. Remote Sens. Lett.*, **10**, 4, 801-805.