Add some error trapping to avoid report failures.

why?

currently using SD < 1e-5 as threshold
* clhs() breaks when sd == 0
* masking 5-95 pctile interval results in not enough data for MDS
* figure out reasonable heuristic (multivariate CV?)

Should account for specific low-variance map units and filter those out instead of barfing.

It would be nice if the output file names were specific to the map units named in config.R.

http://ncss-tech.github.io/soil-pit/projects/MLRA111B_Glynwood_B-slope_Erosion_Northeastern_IN.html

display at bottom of report

https://rmgsc.cr.usgs.gov/outgoing/ecosystems/USdata/

1. tried other MDS methods:
   • MASS::sammon() fails when there are duplicates in the initial configuration
   • vegan::monoMDS() gives similar results as MASS::isoMDS()
   • MASS::isoMDS() is the fastest, most stable algorithm I have tried
   • tried tsne algorithm (https://lvdmaaten.github.io/tsne/), slow and runs out of memory with large datasets

Ideas:

• https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html

perform adjustment based on:

• number of map units
• number of rasters

From here: http://ncss-tech.github.io/AQP/sharpshootR/GIS-summary-by-MU.html

Specifically:

1. flagging polygons with values outside of the 5-95 interval, all rasters?
2. clustering of polygons based on raster values

These operations would have to be done a on a subset of the samples, possibly in conjunction with the multivariate summary.

Currently the only "value added" components are median raster values for each polygon.

Likely in sharpshootR functions.

Ideas: http://datascienceplus.com/parallel-vectorized-operations/

There is no clean way to get devtools and soilReports into the "correct" library paths when:

1. an .Rprofile file is missing
2. HOME is set to H:/

Possible solution: source installRprofile() from GH and run before installing anything.

TODO: test!

1. slope class defs in config.R
2. classify slope map (if there is one) and present as proportions
3. sanity checks

Sometimes and inefficient bandwidth parameter is selected for smoothing in density plots:

New process:

1. load rasters into memory if possible
2. perform cursory grid-based sampling to determine Moran's I of each raster
3. set sampling intensity based on I of each raster

This is all implemented in:

• sharpshootR::sampleRasterStackByMU(..., estimateEffectiveSampleSize=TRUE)
• sharpshootR::Moran_I_ByRaster()
• sharpshootR::ESS_by_Moran_I()

These functions have not been extensively documented or tested.

In the 'region2' folder under branch 'AGB' you will find a skeleton of a pedon summary report that employs Shiny and flexdashboard for user interaction/visualization of pedon data. A static report can be knitted after tweaking parameters of interest.

There are a few different moving parts here that extend beyond the typical case we have for soilReports.

The main workhorse is shiny.Rmd. This is the file that would be run from within Rstudio.

shiny.Rmd sources config.R, utility_functions.R, and main.R.

Upon running shiny.Rmd, an interface showing data from the user's NASIS selected set will appear. config.R allows a few different settings (linework data source, caching, generalized horizons, as well as loading of raster data sources) that are pre-requisites for creating the Shiny interface and its contents.

From the Shiny interface, the user can utilize regex patterns to filter pedon data by MUSYM, taxon name, pedonID or can explicitly specify a list of pedons. Also, they can select a 'modal' pedon from their subset for comparison.

The tabs to the right of the input panel correspond to sections of the pedon summary report with some new additions that display modal pedon data, horizon generalization, etc.

Once a user is happy with the selected set of pedons and is ready to create a permanent record of their work, they hit the "Export" button which will knit a static HTML file containing the same information displayed in the Shiny interface. The R objects that are built in the Shiny interface are passed on to the knit session by creating an R environment containing the relevant dependencies for the template "report.Rmd" to run. "Export" also dumps out copies of the R objects (associated component table, soil profile collections) and tabular output that uniquely identifies the pedons used.

The issue here is how can this style of report be integrated into soilReports?

Let me know what you think if you get a chance to try this out. You'll need pedons and some DMU components in your selected set, update the config to point to your linework, and then once the dashboard loads will need to adjust the default MUSYM pattern for parsing the NASIS selected set.

Caching samples between report runs would speed report development and tweaking of parameters. However, it could cause confusion when a report is re-run after changes have been made to the linework.

Consider setting this flag to FALSE as a default. Also, this may be related to debugging output such as chunk timing.

strip non alpha-numeric chars

Graphical and formalized comparison are difficult due to the VERY high spatial autocorrelation. How can we down-weight the DF accordingly?

Ideas

1. http://www.inside-r.org/packages/cran/SpatialPack/docs/modified.ttest
2. use `clhs` output for bwplots and distributional tests... realistic?
3. weighting via local Moran's I
4. https://github.com/jebyrnes/spatial_correction_lavaan
5. [comparison of methods](http://www.petrkeil.com/?p=1050)
6. [faster ? calculation via ape::Moran.I and distance matrix](http://www.ats.ucla.edu/stat/r/faq/morans_i.htm)

Once we figure out an accurate and scale-able approach, we still need to figure out how to get these adjust sample size numbers into custom.bwplot() via bwplot(). This may require a custom panel function that performs a look-up against the spatial stats summary table.

Ideally, most reports should be configurable via:

• config.R
• interactive report parameters
• parameters specified in a call to render() or in the report YAML header

There are cases where each has benefits. Converting the "region 2" reports to use rmarkdown report parameter names shouldn't take all that long.

When there are MU with low sd, they are filtered before subsetting via clhs. This causes color syncing issues with univariate summaries. Does this also cause syncing issues with labels in the ordination?

Start with something like this:

raster.list <- list(
  continuous=list(
    `Mean Annual Air Temperature (degrees C)`='E:/gis_data/prism/final_MAAT_800m.tif',
    `Mean Annual Precipitation (mm)`='E:/gis_data/prism/final_MAP_mm_800m.tif',
    `Effective Precipitation (mm)`='E:/gis_data/prism/effective_precipitation_800m.tif',
    `Frost-Free Days`='E:/gis_data/prism/ffd_mean_800m.tif',
    `Growing Degree Days (degrees C)`='E:/gis_data/prism/gdd_mean_800m.tif',
    `Elevation (m)`='E:/gis_data/region-2-mu-analysis/elev_30.tif',
    `Slope Gradient (%)`='E:/gis_data/region-2-mu-analysis/slope_30.tif',
    `Annual Beam Radiance (MJ/sq.m)`='E:/gis_data/ca630/beam_rad_sum_mj_30m.tif',
    `(Estimated) MAST (degrees C)`='E:/gis_data/ca630/mast-model.tif',
    `Compound Topographic Index`='E:/gis_data/ca630/tci30.tif',
    `MRVBF`='E:/gis_data/ca630/mrvbf_10.tif',
    `SAGA TWI`='E:/gis_data/ca630/saga_twi_10.tif'
  ),
  categorical=list(
    `Geomorphon Landforms`='L:/Geodata/DEM_derived/forms10.tif',
    `Curvature Classes`='E:/gis_data/ca630/curvature_classes_15.tif'
  ),
  circular=list(`Slope Aspect (degrees)`='E:/gis_data/region-2-mu-analysis/aspect_30.tif')
)

Once #28 is closed it should be possible to iterate over "categorical" type variables and generate a summary, using a template derived from the current geomorphons section.

Might be possible to test for univariate or multivariate separation between map units via supervised classification.

This would make it possible to iterate over items in the "categorical" slice of raster.list and generate a narrative.

this would replace a .Rmd with the latest version, maybe even show differences.

.ipkCRAN() should update existing packages if they are already installed. The current approach will skip any installed packages that match list of packages to install.

Need to fix resetting the factor levels with the nlcd classes.

Add an additional .SHP file to output that includes information regarding polygons with above-threshold proportion of sample values outside 5-95% percentile range for the MU.

Currently, statistics shapefile output contains median values for each raster and the "toCheck" flag which is a ranking based on the amount of samples outside range. This polygon is useful for symbolizing where problematic polygons occur, and also the distribution of abiotic factors (vis-à-vis the median) across a MU extent... but does not tell why a particular polygon is flagged. This can rarely be determined from the median values unless it is a very extreme case.

The new shapefile will follow same format as the stats shapefile (1 column per raster) only instead of medians, contains the "proportion of samples outside range" for each raster data source.

This will reduce the iterative process of looking up polygon IDs to see why they were flagged. Currently that information is only available in the tabular output at the end of the report HTML file.

In addition to reducing iteration between the report and the shapefile display in e.g. ArcMap this will allow symbolizing of MUs based on proportion outside range for INDIVIDUAL data sources rather than an aggregate of all data sources.

This would require a pre-made slab-style database of aggregate data for all MLRA. New configuration options in config.R could allow for specification of 4-5 soil properties from 10-20 possible properties.

Sampling / Aggregating:

library(soilDB)
# iterate over MLRA in current shapefile or local list of codes
x <- fetchKSSL(mlra = '136')

# aggregate, property selection important
a <- slab(x, mlra ~ clay + estimated_ph_h2o + caco3 + bs82)

# save to file

Combine into database:

# assemble samples into single file and distribute

Symbolize in report

library(RColorBrewer)
library(lattice)

# adjust factor labels for MLRA to include number of pedons
pedons.per.mlra <- tapply(site(x)$mlra, site(x)$mlra, length)
a$mlra <- factor(a$mlra, levels=names(pedons.per.mlra), labels=paste(names(pedons.per.mlra), ' (', pedons.per.mlra, ' profiles)', sep=''))

# re-name variables
a$variable <- factor(a$variable, labels=c('Clay %', 'pH 1:1 Water', 'CaCO3 Equiv. (%)', 'Base Sat. pH 8.2'))

# make some nice colors
cols <- brewer.pal('Set1', n=7)

xyplot(
  top ~ p.q50 | variable, groups=mlra, data=a, lower=a$p.q25, upper=a$p.q75, 
  ylim=c(170,-5), alpha=0.25, scales=list(y=list(tick.num=7, alternating=3), x=list(relation='free',alternating=1)),
  panel=panel.depth_function, prepanel=prepanel.depth_function, sync.colors=TRUE, asp=1.5,
  ylab='Depth (cm)', xlab='median bounded by 25th and 75th percentiles', strip=strip.custom(bg=grey(0.85)),
  par.settings=list(superpose.line=list(col=cols, lty=c(1,2,3), lwd=2)),
  auto.key=list(columns=3, title='MLRA', points=FALSE, lines=TRUE),
  sub=paste(length(x), 'profiles')
)

This is typically only a problem when someone has created their input features (mu) with many tables joined to it. Some options:

• filter out all columns except mu.col
• smarter abbreviateNames()

The context of this error is:

poly.check.wide <- dcast(polygons.to.check, pID ~ variable, value.var = 'prop.outside.range')
poly.check.wide[is.na(poly.check.wide)] = 0 #replace NAs with zero (no samples outside 5-95% percentile range)

mu.check <- merge(mu, poly.check.wide, by='pID', all.x=TRUE)
names(mu.check)[-1] <- abbreviateNames(mu.check)

# fix names for printing
names(polygons.to.check)[1] <- mu.col

# print table (removed from report now that shapefile with proportions outside range is generated)
#kable(polygons.to.check[display.idx,], row.names = FALSE) #only shows polys with p.crit > 0.15 in report tabular output

#save a SHP file with prop.outside.range for each polygon and raster data source combination
if(nrow(polygons.to.check) > 0) {
  shp.fname <- paste0('poly-qc-', paste(mu.set, collapse='_'))
  writeOGR(mu.check, dsn='output', layer=shp.fname, driver='ESRI Shapefile', overwrite_layer=TRUE)
  write.csv(mu.check,file=paste0("output\\",shp.fname,".csv")) 
}

NEWS file, or some other update of changes / by report. Perhaps a new list or vector in the report-metadata chunk.

Currently, we assume that every report is comprised of report.Rmd, config.R and setup.R files, but sometimes these are not the only documents that would need to be copied by reportInit().

Examples of this case are:

1. region11 lab summary reports (using custom.R)
2. Shiny interactive pedon summary report (currently under region2 in branch 'AGB')

A few questions for your consideration:

• Should setup.R be extended to include a variable (e.g. files_to_copy; vector of directories and files to copy) that lists additional items to be copied when creating a new report instance?
• Should we distinguish reports that have an interactive/dashboard component (i.e. 'has_dashboard' boolean that indicates presence of shiny.Rmd or equivalent)?
• Should the report metadata (currently stored at top of report.Rmd) be moved to a set of metadata variables in setup.R?

These sections should be more robust in case of NA. This is typically caused by sampling rasters that do not extend beyond the extent of MU polygons.

This is related to #11.

Some Region 2 folks are not convinced that the sampling approach is appropriate--they want to use all pixels within a suite of polygons.

Demonstrate appropriate sampling density via Moran's I and stability of the median as a function of sample size.

Figure out a better way to filter out "id" columns from multivariate analysis.

Building off of #28...

  categorical=list(
    `Geomorphon Landforms`=list(
       data='L:/Geodata/DEM_derived/forms10.tif',
        legend=list(
          levels=1:10,
          labels=c('flat', 'summit', 'ridge', 'shoulder', 'spur', 'slope', 'hollow', 'footslope', 'valley', 'depression')
        )
     ),
  ...
 )

This is likely associated with sharpshootR::sampleRasterStackByMU(), and lower-level raster access as specified in the traceback:

 rgdal::getRasterData(con, offset = offs, region.dim = c(1, nc), 
    band = layers) 
13 .readCellsGDAL(x, uniquecells, layers) 
12 .readCells(x, cells, 1) 
11 .cellValues(object, cells, layer = layer, nl = nl) 
10 .xyValues(x, coordinates(y), ..., df = df) 
9 .local(x, y, ...) 
8 raster::extract(r, s) 
7 raster::extract(r, s) 
6 data.frame(value = raster::extract(r, s), pID = s$pID, sid = s$sid) 
5 (function (r) 
{
    res <- data.frame(value = raster::extract(r, s), pID = s$pID, 
        sid = s$sid) ... 
4 rapply(raster.list, how = "replace", f = function(r) {
    res <- data.frame(value = raster::extract(r, s), pID = s$pID, 
        sid = s$sid)
    return(res) ... 
3 sampleRasterStackByMU(mu, mu.set, mu.col, raster.list, pts.per.acre, 
    estimateEffectiveSampleSize = correct.sample.size) 
2 withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning")) 
1 suppressWarnings(sampleRasterStackByMU(mu, mu.set, mu.col, raster.list, 
    pts.per.acre, estimateEffectiveSampleSize = correct.sample.siz

Some ideas:

• http://r-sig-geo.2731867.n2.nabble.com/Failure-during-raster-IO-td7582421.html
• crop() before extract()

Ideas, based on the WRS package:

• https://garstats.wordpress.com/2016/07/12/shift-function/
• https://freakonometrics.hypotheses.org/4199
• http://www.nicebread.de/comparing-all-quantiles-of-two-distributions-simultaneously/

Citation Here

The ordination figure would be a lot more useful if there were some way to link points in the figure with points on the ground. This could be as simple as preserving coordinates / IDs and outputting sampling locations (cLHS sub-sample). Some degree of interaction is needed when plots are cluttered:

This will be a lot faster.

graphical display of extents