1 Introduction
2 Data & Methods
These methods and data layers are preliminary and subject to change
2.1 Elevation
2.1.1 DEM
Elevation and elevation derivatives from 30-m NASA SRTM. USGS 3-DEM (10m) product not suitable for full study area analysis due to (1) the large area of missing data in Mexico, and (2), the excessively high spatial resolution (massively increasing computational requirements).
SRTM elevation sinks filled prior to calculating slope and aspect.
Should elevation be directly used in the suitability analysis?
2.1.2 Slope
Slope derived from hydrologically conditioned (filled) 30-m SRTM layer using quadratic surface function and a fixed 30-m neighborhood. Slope measured in °.
Higher slopes are less suitable because thinning is both more expensive and more precipitation will end up as runoff.
Slope classified from 1-10 using a continuous function in ArcPro Suitability Mapper.
| Parameter | Setting |
|---|---|
| Function | MSSSmall |
| Mean multiplier | 1 |
| Sddv multiplier | 2 |
| Lower threshold | 0 |
| Value below threshold | 0 |
| Upper threshold | 90 |
| Value above threshold | 0 |
| Invert function | FALSE |
| Save transformed dataset | TRUE |
| Output | Transformed_SRTM_slope |
2.1.3 Aspect
Aspect calculated as with slope. Aspect reference point at N. Pole.
Aspect has a large impact on solar radiation.
Closer to 0 or 360 is desired, low suitability scores for closeness.
Aspect classified from 1-10 using a continuous function in ArcPro Suitability Mapper.
| Parameter | Setting |
|---|---|
| Function | Near |
| Mid Point | 180 |
| Point spread | 0.0011049638968393428 (default) |
| Lower threshold | -1 (flat) |
| Value below threshold | 0 |
| Upper threshold | 360 |
| Value above threshold | 0 |
| Invert function | TRUE |
| Save transformed dataset | TRUE |
| Output | Transformed_SRTM_aspect |
2.2 Precipitation
PRISM normals, 800m resolution. Annual precipitation.
Mean annual precipitation must be higher than 500mm 1990 - 2020
Precipitation classified from 1-10 using a continuous function in ArcPro Suitability Mapper.
NOTE: The logistic growth function may also be a good choice for this dataset. See Logistic Growth function
| Parameter | Setting |
|---|---|
| Function | MSLarge |
| Mean multiplier | 1.68 (approximates 500mm at x-intercept) |
| Sddv multiplier | 1 |
| Lower threshold | 67.33789825439453 (default, minimum) |
| Value below threshold | 0 |
| Upper threshold | 1214.5689697265625 (default, maximum) |
| Value above threshold | 0 |
| Invert function | FALSE |
| Save transformed dataset | TRUE |
| Output | Transformed_PRISM_ppt_30yrnormal_800m |
2.3 Vegetation Characteristics
2.3.1 NLCD 2021 Total Canopy Cover
2.3.2 Landfire
2.4 Soil Hydrology
AZ_Soil_Hydric_Group data layer
Classification Schema
| Class | Count (pixels) | Text | Value |
|---|---|---|---|
| A | 62559472 | Group A soils consist of deep, well drained sands or gravelly sands with high infiltration and low runoff rates. | 10 |
| B | 76665198 | Group B soils consist of deep well drained soils with a moderately fine to moderately coarse texture and a moderate rate of infiltration and runoff. | 8 |
| C | 88491710 | Group C consists of soils with a layer that impedes the downward movement of water or fine textured soils and a slow rate of infiltration. | 5 |
| D | 155095790 | Group D consists of soils with a very slow infiltration rate and high runoff potential. This group is composed of clays that have a high shrink-swell potential, soils with a high water table, soils that have a clay pan or clay layer at or near the surface, and soils that are shallow over nearly impervious material. | 2 |
| A/D | 43192 | Group A/D soils naturally have a very slow infiltration rate due to a high water table but will have high infiltration and low runoff rates if drained. | 7 |
| B/D | 18456 | Group B/D soils naturally have a very slow infiltration rate due to a high water table but will have a moderate rate of infiltration and runoff if drained. | 6 |
| C/D | 217771 | Group C/D soils naturally have a very slow infiltration rate due to a high water table but will have a slow rate of infiltration if drained. | 3 |
Transformed dataset Transformed_AZ_Soil_Hydric_Group
2.5 Depth to Bedrock
There are 2 data layers which represent depth to bedrock and it is not clear which data layer is preferred!
- AZ_BedrockDepth_cm.tif
- 218 m resolution
- UTM 12N, NAVD88 depth (m) positive down
- 0 – 108,273 cm
- Depth to Bedrock WTA
- Classified
- 30 m resolution
- UTM 12N, NAVD88 height (m) positive up
- vertical datum is incorrect. Should be depth (m) positive down
- 0 – 269 cm
- Extremely skewed distribution clustering around 200 cm
2.5.1 Soil vs. Subsurface Geology Weighting Layers
To quantify the differential importance of soils vs. subsurface geology layers for determining suitability two related data layers had to be calculated.
The logic assumes that there are two uniform subsurface layers, soil, and subsurface geology (i.e. geology). However, the weighted importance of these layers is not uniform across space. Where the bedrock is close to the surface, we assume that the soil is the most important layer for ground water storage. Inversely, when the bedrock is extremely deep, we assume that the geology is the more important layer. Our soil layer is measured at a depth of 200cm (2m), and we assume a uniform soil depth across the state. Therefore, the depth to bedrock was divided by 200 to get a depth to bedrock (dtb) in soil units. The first “soil depth” was ascribed to the soil layer, and varies from 0 to 1, while the remaining “soil depth” were attributed to the geology layer, with a range from 0 to 541. Ergo, where the bedrock is deepest, the geology layer is 541 time more influential than the soils layer.
These layers were created in a custom R script using the following raster math, with their resulting outputs.
2.5.1.1 Soils
```{r}
# Where depth to bedrock (dtb) = 0cm, soil multiplier = 0 (no soil)
# Where depth to dtb >= 200cm, soil multiplier = 1 (Full depth of soil)
# Intermediate depths = linear
soilMultiplier = masked
soilMultiplier[soilMultiplier > 200] = 200 # Fix upper limit of soil depth = 200 cm
soilMultiplier = soilMultiplier/200
```
2.5.1.2 Geology
```{r}
# Where dtb < 200cm, geology multiplier = 0 (soil only)
# Where dtb >= 200cm, geology multiplier = dtb/200 (in units of relative soil depth)
geologyMultiplier = masked
geologyMultiplier[geologyMultiplier < 200] = 0
geologyMultiplier = geologyMultiplier/200
```
2.6 Other Data Layers for Consideration
2.6.1 Global Hydrologic Curve Number(GCN250)
https://gee-community-catalog.org/projects/gcn250/?h=hydrologic
The GCN250 is a globally consistent, gridded dataset defining CNs at the 250 m spatial resolution from new global land cover (300 m) and soils data (250 m). GCN250 represents runoff for a combination of the European space agency global land cover dataset for 2015 (ESA CCI-LC) resampled to 250 m and geo-registered with the hydrologic soil group global data product (HYSOGs250m) released in 2018. The potential application of this data includes hydrologic design, land management applications, flood risk assessment, and groundwater recharge modeling. The CN values vary depending on antecedent runoff conditions (ARC), which is affected by the rainfall intensity and duration, total rainfall, soil moisture conditions, cover density, stage of growth, and temperature[.] emphasis mine
2.6.2 Soil Properties 800m
https://gee-community-catalog.org/projects/soilprop/?h=
The data shown here were obtained by aggregating current USDA-NCSS soil survey data (SSURGO back-filled with STATSGO where SSURGO is not available) within 800m² grid cells. This data aggregation technique results in maps that may not match the original data at any given point, and is intended to depict regional trends in soil properties at the statewide scale. emphasis mine
- Pros:
- Lots of relevant data layers, such as:
- Avail. Water Holding Capacity
- Drainage Class
- Sat. Hyd. Conductivity
- Depth to Restrictive Layer
- Hydrologic Group
- Soil Depth
- etc.
- Lots of relevant data layers, such as:
- Cons:
- 800m resolution
- Large data gaps (layer dependent)
2.6.2.1 Alternative layers
gNATSGO (gridded National Soil Survey Geographic Database)
- Pros:
- Authoritative
- Source layer for value added products (including Soil Properties 800m)
- 10m resolution
- Cons:
- Large data gaps across AZ
- 10m resolution
Polaris 30m Probabilistic Soil Properties US
- Pros:
- Continuous data availability (no gaps)
- 30m resolution
- Cons:
- Fewer data layers
- Probabilistic model (increased uncertainty)
2.6.3 Leaf Area Index (LAI)
Leaf Area Index (LAI) can be calculated from Landsat data (30 m resolution), as a proxy of land cover. LAI is a unitless index value which is calculated as a function of the Enhanced Vegetation Index (ENVI) and typically ranges from 0 to 3.5 (citation needed). LAI can be calculated efficiently over a large spatial scale in Google Earth Engine (GEE) using the Javascript API, however LAI will vary seasonally, with updated Landsat data available every 8-days. Additionally, due to our large study area, as well as the logistics of Landsat orbital paths, it takes 6 days to photograph the entire study area. Therefore, any LAI analysis on this scale will necessarily be a mosaic image, over roughly a week. For these reasons, a single LAI image should be viewed with some skepticism, and a seasonal mean or median LAI mosaic image may be more desireable as a proxy for land cover.
Code is available for calculating LAI from Landsat 8/9 imagery in my personal GEE scripts folder, which could be easily modified for purpose. Due to the large number of Landsat images involved the study period must be narrowed prior to analysis (i.e. July LAI, Winter LAI, etc.). Additional post-processing of LAI images may be required outside of GEE.
2.6.3.1 Google Earth Engine LAI Javascript Code
var studyArea = ee.FeatureCollection("projects/ee-travisz09/assets/ATUR/WBDHU8_OuterBoundary_Project"),
landsat9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_TOA");
// Map setup
Map.centerObject(studyArea, 6);
Map.addLayer(studyArea, {}, 'Study Area', false);
// Paint study area outline
var empty = ee.Image().byte(); // An empty layer to paint on
var studyAreaOutline = empty.paint(studyArea, 'black', 2);
Map.addLayer(studyAreaOutline, {}, 'Study Area (outline)');
// Study period
var startDate = ee.Date('2023-01-01');
var endDate = ee.Date('2024-01-01'); // Exclusive end date
var timeDif = endDate.difference(startDate, 'day');
// var interval = 8; // days
// print(timeDif);
// Landsat 8/9
var cloudFilter = 10 // 10% max cloud cover
var landsat = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2')) // Combine Lst-9 and Lst-8
.filterDate(startDate, endDate)
.filterBounds(studyArea)
.sort('system:time_start') // Sort by date
.filter(ee.Filter.lt('CLOUD_COVER', cloudFilter));
// Lst scaling factor function
var scaleLst = function(col) {
// Map over images in collection
var scaled = col.map(function(img) {
var opticalBands = img.select('SR_B.').multiply(0.0000275).add(-0.2);
var thermalBands = img.select('ST_B.*').multiply(0.00341802).add(149.0);
return img.addBands(opticalBands, null, true) // Overwrite unscaled bands
.addBands(thermalBands, null, true);
});
return scaled;
};
// print(landsat.first()); // Check the Lst data before scaling
landsat = scaleLst(landsat); // Apply scaling factor
// print(landsat.first()); // Check Lst data after scaling
// NDVI, EVI, LAI, index function
var calcIndices = function(col) {
var indices = col.map(function(img) {
var ndvi = img.normalizedDifference(['SR_B5', 'SR_B4'])
.rename('ndvi');
var evi = img.expression(
'2.5 * ((nir - red) / (nir + 6 * red - 7.5 * blue + 1))',
{
'nir': img.select('SR_B5'),
'red': img.select('SR_B4'),
'blue': img.select('SR_B2')
}).rename('evi');
var lai = evi.expression(
'3.618 * evi - 0.118',
{
'evi': evi.select('evi')
}).rename('lai');
return img.addBands(ndvi)
.addBands(evi)
.addBands(lai);
});
return indices;
};
// Apply indices
landsat = calcIndices(landsat);
// print(landsat.first()); // sanity check
// Isolate lai bands
var lai = landsat.map(function(img) {
return img.select('lai');
});
Map.addLayer(lai.limit(20), {
min: 0,
max: 3.5,
palette: ['red', 'white', 'green']
}, 'LAI');
print(lai);