README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "##",
  fig.path = "man/figures/README-",
  out.width = "100%",
  warning = FALSE,
  message = FALSE
)
```

# spEcula <img src="man/figures/logo.png" align="right" height="90"/>

<!-- badges: start -->

[![CRAN](https://www.r-pkg.org/badges/version/spEcula)](https://CRAN.R-project.org/package=spEcula)
[![r-universe](https://spatlyu.r-universe.dev/badges/spEcula)](https://spatlyu.r-universe.dev/spEcula)

<!-- badges: end -->

The goal of **spEcula** is to make it easier to use **R** for **spatial prediction** based on **various spatial relationships** (e.g. spatial dependence, spatial heterogeneity and geographical similarity).

## Overview

Full document of the most recent release of **spEcula** is online:
<https://spatlyu.github.io/spEcula/>

Current models and functions provided by **spEcula** are:

|   **spatial prediction method**   |   **spEcula function**   |   **support status**   |
|-----------------------------------|--------------------------|------------------------|
|  Geographically Optimal Similarity  |        `gos()`         |          ✔️        |
|  Sandwich Mapping Model  |        `sandwich()`         |          ✔️        |


## Installation

You can install the development version of *spEcula* from [*github*](https://github.com/SpatLyu/spEcula):

``` r
# install.packages("devtools")
devtools::install_github("SpatLyu/spEcula",
                         build_vignettes = T,
                         dep = T)
```

or install *spEcula* from [*r-universe*](https://spatlyu.r-universe.dev/spEcula):

```r
install.packages('spEcula',
                 repos = c("https://spatlyu.r-universe.dev",
                           "https://cran.rstudio.com/"),
                 dep = TRUE)
```

## Example

### Geographically Optimal Similarity (GOS) model

```{r example_gos}
library(spEcula)
data(zn)
data(grid)

zn$Zn = log(zn$Zn)
tictoc::tic()
g1 = gos(Zn ~ Slope + Water + NDVI  + SOC + pH + Road + Mine,
         data = zn, newdata = grid, kappa = 0.08,cores = 6)
tictoc::toc()
g1$pred = exp(g1$pred)
grid$pred = g1$pred
grid$uc99 = g1$`uncertainty99`
g1
```

```{r}
library(ggplot2)
library(cowplot)
library(viridis)
```


```{r gos_result,fig.width=9.5,fig.height=3.5}
f1 = ggplot(grid, aes(x = Lon, y = Lat, fill = pred)) +
  geom_tile() +
  scale_fill_viridis(option="magma", direction = -1) + 
  coord_equal() +
  labs(fill='Prediction') +
  theme_bw() 
f2 = ggplot(grid, aes(x = Lon, y = Lat, fill = uc99)) +
  geom_tile() +
  scale_fill_viridis(option="mako", direction = -1) + 
  coord_equal() +
  labs(fill=bquote(Uncertainty~(zeta==0.99))) +
  theme_bw() 

plot_grid(f1,f2,nrow = 1,label_fontfamily = 'serif',
          labels = paste0('(',letters[1:2],')'),
          label_fontface = 'plain',label_size = 10,
          hjust = -1.5,align = 'hv')  -> p
p
```

### Sandwich Mapping Model

```{r sandwich,fig.width=5.5,fig.height=3.5}
library(sf)
library(tidyverse)
library(spEcula)
simpath = system.file("extdata", "sim.gpkg", package="spEcula")
sampling = read_sf(simpath,layer = 'sim_sampling')
ssh = read_sf(simpath,layer = 'sim_ssh')
reporting = read_sf(simpath,layer = 'sim_reporting')

sampling_zone = sampling %>%
    st_join(ssh['X']) %>%
    st_drop_geometry()

library(ggpubr)

ggerrorplot(sampling_zone, x = "X", y = "Value",
            desc_stat = "mean_sd", color = "black",
            add = "violin", add.params = list(color = "darkgray")) +
  geom_text(data = summarise(sampling_zone,vmean = mean(Value),.by = X), 
            aes(x = X, y = vmean, label = round(vmean,2)), 
            vjust = -0.5, hjust = -0.15, color = "black",size = 3) +
  scale_x_discrete(labels = LETTERS[1:4]) +
  theme(axis.title.x = element_blank())

sim_est = sandwich(sampling = sampling,stratification = ssh,reporting = reporting,
                   sampling_attr = 'Value',ssh_zone = 'X',reporting_id = 'Y',
                   weight_type = 'area')
sim_est
```

```{r sandwich_sim_est,fig.width=9.75, fig.height=3.5}
library(cowplot) 

f1 = ggplot(data = sim_est, aes(fill = sandwichest_mean), 
            color = "darkgray") +
  geom_sf() + 
  labs(fill='mean') +
  scale_fill_gradient(low = "#f0bc9c", high = "red",
                      breaks = range(sim_est$sandwichest_mean)) +
  theme_bw() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    legend.position = 'right',
    legend.background = element_rect(fill = 'transparent',color = NA)
  )

f2 = ggplot(data = sim_est, aes(fill = sandwichest_standarderror), 
            color = "darkgray") +
  geom_sf() + 
  labs(fill='se') +
  scale_fill_gradient(low = "#b6edf0", high = "blue",
                      breaks = range(sim_est$sandwichest_standarderror)) +
  theme_bw() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    legend.position = 'right',
    legend.background = element_rect(fill = 'transparent',color = NA)
  )

plot_grid(f1, f2, nrow = 1,label_fontfamily = 'serif',
          labels = paste0('(',letters[1:4],')'),
          label_fontface = 'plain',label_size = 10,
          hjust = 0.05,align = 'hv')  -> p
p
```