scbrown86

get_point_at_distance <- function(lon, lat, d, bearing, R = 6378137) {
  # lat: initial latitude, in degrees
  # lon: initial longitude, in degrees
  # d: target distance from initial point (in m)
  # bearing: (true) heading in degrees
  # R: mean radius of earth (in m)
  # Returns new lat/lon coordinate {d} m from initial, in degrees
  stopifnot("lon value not in range {-180, 180}" = !any(lon < -180 | lon > 180))
  stopifnot("lat value not in range {-90, 90}" = !any(lat < -90 | lat > 90))
  stopifnot("bearing not in range {0, 360}" = !any(bearing < 0 | bearing > 360))

scbrown86

# Small function for correctly plotting a global map whith a rotated projection (i.e. non standard centre longitude)
## from here - https://stackoverflow.com/a/68313482
rotate_prj <- function(x, crs) {
  stopifnot(inherits(x, what = "sf"))
  stopifnot(inherits(st_crs(crs), "crs"))
  # make the data valid before doing anything else!
  x <- st_make_valid(x)
  # determine the rotated/centre longitude from crs
  lon <- sapply(strsplit(as.character(st_crs(crs)[2]), "\n"), trimws)
  lon <- lon[which(grepl(pattern = "Longitude of natural origin", x = lon))]

scbrown86

##%######################################################%##
#                                                          #
####           Extracts probability estimates           ####
####         from calibration/validation points         ####
####          for n-dimensional hypervolumes,           ####
####               and compares these to                ####
####    background values from the same hypervolume.    ####
####         Values are then rescaled to {0, 1}         ####
####          based on 10% omission rate, with          ####
####        upper values capped at 90%. AUC and         ####

scbrown86

ortho_plot <- function(poly, lat = 55, lon = 20, lwd = 0.25, bgc = "#bfd7ea") {
  require(sf); require(mapview)
  sf_use_s2(FALSE)
  # Define the orthographic projection
  # Choose lat_0 with -90 <= lat_0 <= 90 and lon_0 with -180 <= lon_0 <= 180
  ortho <- sprintf('+proj=ortho +lat_0=%s +lon_0=%s +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs', lat, lon)
  assign(x = "ortho_crs", value = ortho, envir = .GlobalEnv)
  # Define the polygon that will help you finding the "blade"
  # to split what lies within and without your projection
  # buffer is == Semimajor radius of the ellipsoid axis

scbrown86

# https://doi.org/10.1023/A:1015508322413
# The refet() function assumes that the rasters are ordered Jan-Dec
# if a different order is required, make sure a z-index is provided for
# the Tmin raster. The refet() function will then calculate
# julian dates and month lengths off this z-index.

library(raster)
library(sirad)

Tmin <- resample(crop(x = getData("worldclim", res = 10, var="tmin"),

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)

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

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

# 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

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