## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
library(intkrige)
data(ohtemp)
head(ohtemp)

## ----ohioMap, eval = FALSE, echo = FALSE, fig.width = 7, fig.height = 5-------
# 
# library(ggmap)
# # # Load the Ohio River Basin Shapefile
# # ohMap <- rgdal::readOGR(dsn = "ohMap", layer = "ohMap")
# #
# # # Simplify the shapefile
# # ohMap <- rgeos::gSimplify(ohMap, tol = .05, topologyPreserve = TRUE)
# data(ohMap)
# 
# # key has been removed for security reasons
# # register_google(key = "insert key here")
# register_google(key = "<my key>")
# center = apply(ohMap@bbox, 1, mean)
# 
# oh.google <- ggmap::get_googlemap(center = center, zoom = 6,
#                                     maptype = "terrain",
#                                     color = "bw",
#                                   style = "feature:all|element:labels|visibility:off|
#                            &style=feature:road|element:all|visibility:off|
#                            &style=feature:landscape.man_made|element:all|visibility:off|")
# 
# map <- ggmap(oh.google) +
#   xlim(ohMap@bbox[1, 1], ohMap@bbox[1, 2]) +
#   xlab("Longitude (easting)") +
#   ylim(ohMap@bbox[2, 1], ohMap@bbox[2, 2]) +
#   ylab("Latitude (northing)") +
#   geom_polygon(data = ohMap, aes(x = long, y = lat),
#                colour = "gray50", fill = NA, lwd = 1.1, lty = 1) +
#   geom_point(data = as.data.frame(ohtemp), aes(x = LONGITUDE, y = LATITUDE),
#              pch = 16, alpha = 0.9) +
#   theme(legend.title = element_blank(),
#         legend.text = element_text(size = 16),
#         axis.text = element_text(size = 16),
#         axis.title = element_text(size = 16))
# 
# pdf("vignettes/ohmap.png", fig.width = 4, fig.height = 4)
# map
# dev.off()

## -----------------------------------------------------------------------------
# First, create a SpatialPointsDataFrame in the usual way
tproj <- sp::CRS(SRS_string = "EPSG:4326")
sp::coordinates(ohtemp) <- c("LONGITUDE", "LATITUDE")
sp::proj4string(ohtemp) <- tproj
interval(ohtemp) <- c("minm", "maxm")

head(ohtemp)

## -----------------------------------------------------------------------------
interval(ohtemp) <- log(interval(ohtemp))
head(ohtemp)

## ----fig.align='center', fig.width = 7, fig.height = 5------------------------
# Revert back to the standard interval
interval(ohtemp) <- exp(interval(ohtemp))
varios <- intvariogram(ohtemp, cutoff = 500)

plot(varios)

## ----fig.align='center', fig.width = 7, fig.height = 5------------------------
varioFit <- fit.intvariogram(varios, models = gstat::vgm(c("Lin", "Sph", "Sph")))
varioFit
intvCheck(varios, varioFit)

## ----fig.align='center', fig.width = 7, fig.height = 5------------------------
# Replace non-convergent variogram fit with a surrogate that 
# contains a reasonable range and sill. 
varioFit[[1]] <- gstat::vgm(psill = 350, nugget = 4.608542, 
                            range = 600, model = "Sph")

intvCheck(varios, varioFit)

## ----fig.align = 'center', fig.width = 7, fig.height = 5----------------------
# Include the Ohio river basin shapefile
data(ohMap)

# New location data preparation
lon <- seq(-89.26637, -77.83937, length.out = 10)
lat <- seq(35.31332, 42.44983, length.out = 10)
newlocations <- expand.grid(lon, lat)
colnames(newlocations) <- c("lon", "lat")
sp::coordinates(newlocations) <- c("lon", "lat")
sp::proj4string(newlocations) <- tproj
sp::gridded(newlocations) <- TRUE

# Adjust r and theta to ensure answers remain in feasible region.
preds <- intkrige(ohtemp, newlocations, varioFit, 
                  A = c(1, 1, 0.5), r = 200, eta = 0.9, maxp = 225)

plot(preds, beside = FALSE, circleCol = "gray") + 
  latticeExtra::layer(sp.lines(ohMap, col = "white"))
plot(preds, beside = TRUE)