Skip to content
Lemos and Sanso' (2012) model.
R
Branch: master
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
R
figs
man
.Rbuildignore
.gitignore
DESCRIPTION
NAMESPACE
README.md
Scallops.Rproj

README.md

title author date output
Scallops
Ricardo T. Lemos
7/21/2018
html_document

Description

Scallops is an R package for geostatistical modeling. It is inspired by the paper Lemos and Sanso' (2012).

Installation

You can install the github development version via

devtools::install_github('rtlemos/Scallops')

Example

Let us load the famous "scallop" data and plot the georeferenced log-catches, along with a grid that we will use to interpolate those values.

library(SemiPar)
library(ggmap)
data(scallop)
scallop$logcatch = log(scallop$tot.catch + 1)
dpc_grid = get_grid(c(-73.75, -71.25), c(38.5, 41), 0.5)

map.ny = get_map(location = c(-72.5, 39.75), zoom = 8)
ggmap(map.ny) + 
  geom_point(data = scallop, aes(x = longitude, y = latitude, size = logcatch), shape = 1) +
  geom_point(data = dpc_grid$coord, aes(x = lon, y = lat), shape = 3, color = 'red')

We will be placing one Gaussian variable over each gridpoint. To have an idea of its influence on the interpolation, let us pick one in the middle of the plot, (39.5N, -72.75E), and depict how its weight changes across space.

get_influence_plot(dpc_grid, lat = 39.5, lon = -72.75) +
    geom_point(data = scallop, aes(x = longitude, y = latitude, size = logcatch), shape = 1)

The plot above is for an "isotropic" kernel (actually, the kernel is not exactly isotropic, because we are using a latitude-longitude coordinate system). By default, the kernel's range corresponds to twice the grid spacing.

Let us now fit the model and depict the spatial variation of the posterior weight for the same Gaussian variable.

fit = get_mcmc(nburn = 10, nsample = 1000,
               s = data.frame(lon = scallop$longitude, lat = scallop$latitude), 
               dpc_grid = dpc_grid, y = scallop$logcatch, 
               priors = get_priors(dpc_grid = dpc_grid, precision_diagonal_value = 0.01),
               seed = 1)
get_influence_plot(dpc_grid, lat = 39.5, lon = -72.75, fit = fit) +
    geom_point(data = scallop, aes(x = longitude, y = latitude, size = logcatch), shape = 1)

We can also look at the shape of the kernels at all gridpoints.

get_ellipses_plot(dpc_grid, fit)

And finally, let us look at the interpolation.

get_interpolation_plot(obs_coord = scallop, dpc_grid = dpc_grid, fit = fit, contour_binwidth = 1) +
                       geom_text(data = scallop, aes(x = longitude, y = latitude, label = round(logcatch)))

You can’t perform that action at this time.