scbrown86

# Use the Cumulative Relative Frequency Curve hotspot method
# to identify hotspots of high values

# Ref:
## Bartolino, V., Maiorano, L., & Colloca, F. (2011).
## A frequency distribution approach to hotspot identification.
## Population Ecology, 53(2), 351-359. doi:10.1007/s10144-010-0229-2

# example array data
# could do with raster/terra, but wanted example in base r

scbrown86

# Use the Cumulative Relative Frequency Curve hotspot method
# to identify hotspots of high values

# Ref:
## Bartolino, V., Maiorano, L., & Colloca, F. (2011).
## A frequency distribution approach to hotspot identification.
## Population Ecology, 53(2), 351-359. doi:10.1007/s10144-010-0229-2


library(raster)

scbrown86

# Convert raster to polygon using GDAL
## edited 2019-03-27 to use sf::st_read
polygonise <- function(x, outshape=NULL, gdalformat = 'ESRI Shapefile', 
                        pypath=NULL, readpoly=TRUE, quietish=TRUE) {
  ## x: an R Raster layer, or the file path to a raster file recognised by GDAL
  ## outshape: the path to the output shapefile (if NULL, a temporary file will be created)
  ## gdalformat: the desired OGR vector format. Currently only supports ESRI Shapefile!
  ## pypath: the path to gdal_polygonize.py (if NULL, an attempt will be made to determine the location
  ## readpoly: should the polygon shapefile be read back into R, and returned by this function? (logical)
  ## quietish: should (some) messages be suppressed? (logical)

scbrown86

## Function for returning weighted centroid location
## currently hard-coded to return mean, q10, q20, q80 & q90.
weightedCentre <- function(x, y, z) {
  require(matrixStats); require(Hmisc)
  if (anyNA(c(x, y, z))) {
    stop("There are missing values present in x, y, or z")
  }
  if (length(z) != length(x)) {
    stop("Number of weights supplied not equal to number of coordinates")
  }

scbrown86

## GISFrag metric
## https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19950017676.pdf
## 1) Produce a proximity (distance) map between 'patches'
## 2) GISFrag == mean of all values on the proximity map
## 3) Large GISFrag values reflect low forest fragmentation, low values == high fragmentation

gisFrag <- function(x, ...) {
  ## x needs to be a raster where cells with suitable habitat are coded as 1
  ## unsuitable cells coded with 0
  ## extract cell numbers for suitable cells

scbrown86

x<- seq (from = 1, to = 100, by=1)
y<-cumsum(rnorm(100, mean = 1, sd = 20))
normalized = (y-min(y))/(max(y)-min(y))
y.mod <- (normalized*2)+10

scbrown86

# somewhat hackish solution to:
# https://twitter.com/EamonCaddigan/status/646759751242620928
# based mostly on copy/pasting from ggplot2 geom_violin source:
# https://github.com/hadley/ggplot2/blob/master/R/geom-violin.r

library(ggplot2)
library(dplyr)


"%||%" <- function(a, b) {

scbrown86

library(tidyverse)
library(scales)
data(diamonds)

diamonds %>% 
  filter(str_detect(cut, "Fair|Ideal")) %>% 
  ggplot(aes(price, carat)) +
  geom_point(color = "skyblue", alpha = 0.5) +
  facet_wrap(~cut, strip.position = "bottom") + 
  scale_x_continuous(labels = comma) +

scbrown86

cd /mnt/c/Users/Stu/Desktop/CMIP5/

#CCSM4
cd CCSM4
cdo -ensmean tas_Amon_CCSM4_rcp26_r* CCSM4_rcp26_ensAvg.nc #6 seconds

#CSIRO
cd ../CSIRO-Mk3.6.0/
cdo -ensmean tas_Amon_CSIRO-Mk3-6-0_rcp26_r* CSIRO-Mk3-6-0_rcp26_ensAvg.nc #2.5 seconds

scbrown86

m.stat <- function(mod, obs, ...) {
  # model field = x; observed field = y
  # mse the mean square error between x and y,
  # V and G are spatial variance and domain mean of the respective fields
  # https://doi.org/10.1002/(SICI)1097-0088(199604)16:4%3C379::AID-JOC18%3E3.0.CO;2-U
  obs = na.omit(obs)
  mod = na.omit(mod)
  stopifnot(length(obs) == length(mod))
  se = (obs - mod)^2
  mse = mean(se)