---
title: "Point-referenced spatial data"
subtitle: "Spatial Statistics - BIOSTAT 696/896"
author: "Michele Peruzzi"
institute: "University of Michigan"
execute:
  cache: true
format:
  revealjs:
    auto-stretch: true
    center: true
    html-math-method: katex
    transition: none
    slide-number: true
    callout-appearance: minimal
    code-line-numbers: true
    fontsize: 16pt
    output-file: l03_pointreferenced.html
    theme: [simple, custom.css] #, web/pulse_colors.css] # alternative themes (subset): default, night, dark
    embed-resources: false
    echo: true
    fig-dpi: 150
    # incremental: true  # Remove comment if you like incremental bullet points
    footer: "[Home page](index.html)"
knitr:
  opts_chunk:
    out.width: "60%"
editor: 
  markdown: 
    wrap: 72
editor_options: 
  chunk_output_type: console
---

## Recap: point-referenced spatial data

-   We observe $Y(s)$ at location $s$
-   [What is $s$ in this case?]{.orange}
-   [What are some examples of point-referenced data?]{.orange}

## Point-referenced spatial data

-   Observe $Y(s)$ at $s \in \cal L = \{s_1, \dots, s_n\}$
-   The set $\cal L$ is assumed fixed and known, ie the locations are
    *not random*

### Goals

-   Learn/estimate how $Y(\cdot)$ changes in space. Model the spatial
    dependence
-   Generalize linear regression: drop i.i.d. assumption in favor of spatial dependence
-   Predict $Y(\cdot)$ at locations that were *not* observed

### Example

-   How does air pollution spread geographically?
-   Considering the impact of air pollution on health, how should we
    model the spread of respiratory disease?
-   Using a fixed set of observations of air pollution, how can we
    estimate the level of air pollution *everywhere else*?

## Exploratory analysis and simulation

-   Simulating a fixed set of spatial locations
-   This is just for simulation. Remember with point-referenced data the
    coordinates are given and non-random
-   Remember to set the seed for reproducible analyses!!

```{r echo=FALSE}
library(scico)
```

```{r}
set.seed(2024) # set a seed

n_obs <- 200
x_coord <- runif(n_obs, 0, 1)
y_coord <- runif(n_obs, 0, 1)
coords <- data.frame(x_coord, y_coord)
```

Plotting

```{r}
library(tidyverse)

ggplot(coords, aes(x_coord, y_coord)) +
  geom_point(size=.7) +
  theme_minimal()
```

## Exploratory analysis and simulation

-   We just have locations now, there is no data. Let's generate some
    data for those locations
-   No need to set the seed here. The script will run from the beginning
    and do all RNGs in sequence
-   You *can* set the seed here if you expect to run chunks of code in
    different sequences

```{r}
#set.seed(1234) 
y <- rnorm(n_obs, 0, 1)
```

-   [How many random variables are we generating?]{.orange}
-   [What do you expect to see when we plot `y` on the coordinates we
    had before? Why?]{.orange}

## Exploratory analysis and simulation

```{r}
df <- coords %>% mutate(y=y)

ggplot(df, aes(x_coord, y_coord, color=y)) +
  geom_point() +
  theme_minimal() +
  scale_color_viridis_c()
```

## Exploratory analysis and simulation

Make it easier to see there is no spatial dependence by

-   increasing sample size

```{r}
n_obs <- 20000
x_coord <- runif(n_obs, 0, 1)
y_coord <- runif(n_obs, 0, 1)
coords <- data.frame(x_coord, y_coord)
y <- rnorm(n_obs, 0, 1)
df <- coords %>% mutate(y=y)

ggplot(df, aes(x_coord, y_coord, color=y)) +
  geom_point(size=.5) +
  theme_minimal() +
  scale_color_viridis_c()
```

## Exploratory analysis and simulation

Make it easier to see there is no spatial dependence by

-   plotting on a regular grid of spatial locations
-   regular grid = each data point can be plotted as a pixel on a raster
    image

```{r}
coords <- expand.grid(xlocs <- seq(0,1,length.out=100), xlocs) 
n_obs <- nrow(coords)
y <- rnorm(n_obs, 0, 1)
df <- coords %>% mutate(y=y)

ggplot(df, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  theme_minimal() +
  scale_fill_scico(palette="batlowK")
```

## Exploratory analysis and simulation

-   We have seen simulated spatial data without spatial dependence
-   [What would we expect to see *if* there was spatial
    dependence?]{.orange}
-   [Can you think of ways we could simulate data with spatial
    dependence?]{.orange}

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 1

```{r echo=F}
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- coords %>% mutate(y = Var1)

plotlist <- list()
i <- 1
plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 2

```{r echo=F}
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- coords %>% mutate(y = (Var1-.2)^3 + (2*Var2-1)^2/5)

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 3

```{r echo=F}
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- coords %>% mutate(y = sin(Var1*2*pi) + cos(Var2*2*pi - 1))

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 4

```{r echo=F}
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- coords %>% mutate(y = sin(Var1*10*pi) + cos(Var2*20*pi))

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 5

```{r echo=F}
set.seed(2024)
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- meshed::rmeshedgp(coords, c(10,1.5,1), block_size=200) %>% dplyr::rename(y=w)

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 6

```{r echo=F}
set.seed(2024)
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- meshed::rmeshedgp(coords, c(1,0.5,1), block_size=200) %>% dplyr::rename(y=w)

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 7

```{r echo=F}
set.seed(2024)
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- meshed::rmeshedgp(coords, c(70,1.5,5), block_size=200) %>% dplyr::rename(y=w)

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 8

```{r echo=F}
set.seed(2024)
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- meshed::rmeshedgp(coords, c(3,1.5,5), block_size=200) %>% dplyr::rename(y=w)

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 9

```{r echo=F}
set.seed(2024)
coords <- expand.grid(x <- seq(0,1,length.out=200), x)

warp1 <- meshed::rmeshedgp(coords, c(3,1.5,.01), block_size=100)
warp2 <- meshed::rmeshedgp(coords, c(4,1.5,.01), block_size=100)

newcoords <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(Var1 = Var1 + warp1$w, Var2 = Var2 + warp2$w)

w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   What does spatial dependence "look like"? Example 10

```{r echo=F}
set.seed(2024)
coords <- expand.grid(x <- seq(0,1,length.out=200), x)

warp1 <- meshed::rmeshedgp(coords, c(10,1.5,.02), block_size=100)
warp2 <- meshed::rmeshedgp(coords, c(2,1.5,.02), block_size=100)

newcoords <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(Var1 = Var1 + warp1$w, Var2 = Var2 + warp2$w)

w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)

plotlist[[i]] <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
plotlist[[i]]
i <- i+1
```

## Exploratory analysis and simulation

-   [What common feature can we recognize?]{.orange}

```{r echo=FALSE}
library(gridExtra)

grid.arrange(grobs=plotlist %>% lapply(\(x) x + theme_void() + theme(legend.position="none")),
             nrow=2)
```

## We only observe 1 sample

What if we could "re-run" the data-generation? Different data, same process

```{r echo=FALSE}
set.seed(2024)
w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)
labeldf <- data.frame(Var1=c(0.15,0.185), Var2=c(0.44,0.42), y=0, lab=c("Y(s[1])", "Y(s[2])"))

ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal() + 
  geom_point(data=labeldf, color="white") +
  geom_text(data=labeldf %>% 
      mutate(
        Var1=c(0.12, 0.22),
        Var2=c(0.4, 0.45)), aes(label=lab), color="white", parse=TRUE)
```

## We only observe 1 sample

What if we could "re-run" the data-generation? Different data, same process

```{r echo=FALSE}
set.seed(1)
w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)
labeldf <- data.frame(Var1=c(0.15,0.185), Var2=c(0.44,0.42), y=0, lab=c("Y(s[1])", "Y(s[2])"))

ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal() + 
  geom_point(data=labeldf, color="white") +
  geom_text(data=labeldf %>% 
      mutate(
        Var1=c(0.12, 0.22),
        Var2=c(0.4, 0.45)), aes(label=lab), color="white", parse=TRUE)
```


## We only observe 1 sample

What if we could "re-run" the data-generation? Different data, same process

```{r echo=FALSE}
set.seed(2)
w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)
labeldf <- data.frame(Var1=c(0.15,0.185), Var2=c(0.44,0.42), y=0, lab=c("Y(s[1])", "Y(s[2])"))

ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal() + 
  geom_point(data=labeldf, color="white") +
  geom_text(data=labeldf %>% 
      mutate(
        Var1=c(0.12, 0.22),
        Var2=c(0.4, 0.45)), aes(label=lab), color="white", parse=TRUE)
```


## We only observe 1 sample

What if we could "re-run" the data-generation? Different data, same process

```{r echo=FALSE}
set.seed(3)
w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)
labeldf <- data.frame(Var1=c(0.15,0.185), Var2=c(0.44,0.42), y=0, lab=c("Y(s[1])", "Y(s[2])"))

ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal() + 
  geom_point(data=labeldf, color="white") +
  geom_text(data=labeldf %>% 
      mutate(
        Var1=c(0.12, 0.22),
        Var2=c(0.4, 0.45)), aes(label=lab), color="white", parse=TRUE)
```



## We only observe 1 sample

What if we could "re-run" the data-generation? Different data, same process

```{r echo=FALSE}
set.seed(4)
w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)
labeldf <- data.frame(Var1=c(0.15,0.185), Var2=c(0.44,0.42), y=0, lab=c("Y(s[1])", "Y(s[2])"))

ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal() + 
  geom_point(data=labeldf, color="white") +
  geom_text(data=labeldf %>% 
      mutate(
        Var1=c(0.12, 0.22),
        Var2=c(0.4, 0.45)), aes(label=lab), color="white", parse=TRUE)
```


## We only observe 1 sample

What if we could "re-run" the data-generation? Different data, same process

```{r echo=FALSE}
set.seed(5)
w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)
labeldf <- data.frame(Var1=c(0.15,0.185), Var2=c(0.44,0.42), y=0, lab=c("Y(s[1])", "Y(s[2])"))

ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal() + 
  geom_point(data=labeldf, color="white") +
  geom_text(data=labeldf %>% 
      mutate(
        Var1=c(0.12, 0.22),
        Var2=c(0.4, 0.45)), aes(label=lab), color="white", parse=TRUE)
```


## We only observe 1 sample

What if we could "re-run" the data-generation? Different data, same process

```{r echo=FALSE}
set.seed(6)
w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)
labeldf <- data.frame(Var1=c(0.15,0.185), Var2=c(0.44,0.42), y=0, lab=c("Y(s[1])", "Y(s[2])"))

ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal() + 
  geom_point(data=labeldf, color="white") +
  geom_text(data=labeldf %>% 
      mutate(
        Var1=c(0.12, 0.22),
        Var2=c(0.4, 0.45)), aes(label=lab), color="white", parse=TRUE)
```


## Exploratory analysis and simulation

-   How "similar" are the random variables $Y(s_1)$ and $Y(s_2)$?
-   In other words, how large do we expect this is:
    $| p(Y(s_2) | Y(s_1)) - p(Y(s_2)) |$?
-   Easier? [How large do we expect this is:
    $Cov( Y(s_1), Y(s_2) )$?]{.orange}

```{r echo=FALSE}
labeldf <- data.frame(Var1=c(0.2,0.7), Var2=c(0.3,0.6), y=0, lab=c("Y(s[1])", "Y(s[2])"))

grid.arrange(grobs=plotlist %>% lapply(\(x) x + 
                                   theme_void() + 
                                   theme(legend.position="none") + 
                                   geom_point(data=labeldf, color="white") +
                                   geom_text(data=labeldf, aes(label=lab), color="white", parse=TRUE, nudge_x=0.02, nudge_y=0.06)),
             nrow=2)

```

## Exploratory analysis and simulation

-   How "similar" are the random variables $Y(s_3)$ and $Y(s_4)$?
-   In other words, how large do we expect this is:
    $| p(Y(s_4) | Y(s_3)) - p(Y(s_4)) |$?
-   Easier? How large do we expect this is: $Cov( Y(s_3), Y(s_4) )$?

```{r echo=FALSE}
labeldf <- data.frame(Var1=c(0.3,0.35), Var2=c(0.45,0.5), y=0, lab=c("Y(s[3])", "Y(s[4])"))

grid.arrange(grobs=plotlist %>% lapply(\(x) x + 
                                   theme_void() + 
                                   theme(legend.position="none") + 
                                   geom_point(data=labeldf, color="white") +
                                   geom_text(data=labeldf %>% 
                                               mutate(
                                                 Var1=c(0.27, 0.43),
                                                 Var2=c(0.4, 0.55)), aes(label=lab), color="white", parse=TRUE)),
             nrow=2)

```

## Exploratory analysis and simulation

-   Distance between random variables (largely) drives how related they
    are
-   If the variables are close to each other in space, then
    $Cov(Y(s), Y(s'))$ should be larger
-   The question of *what model* to use for $Cov(Y(s), Y(s'))$ is
    crucial
-   If we let $Cov(Y(s), Y(s')) = g(s, s')$, how do we choose
    $g(\cdot)$?
-   For example, we could model $g(\cdot)$ as a positive decreasing
    function of $\| s-s' \|$
-   We can think of many other alternative models. More on this later...

## Exploratory analysis: empirical semivariogram

-   How do we measure spatial association in data *before we attempt any
    sort of modeling*?
-   Call $Y(s_i)$ and $Y(s_j)$ the r.v.'s at $s_i$ and $s_j$
-   Suppose the domain is $D = [0,1]^2$
-   Partition the line $(0, m) \subset \Re$ into $K$ disjoint intervals
-   $I_1 = (0=m_0, m_1]$, $I_2 =(m_2, m_3]$, $\dots$,
    $I_K = (m_{K-1}, m_K=m]$.
-   Define $t_k = \frac{m_k - m_{k-1}}{2}$ (the midpoint of the $k$th
    interval)
-   Define $N(t_k) = \{ (s_i, s_j) : \| s_i - s_j \| \in I_k \}$
-   $N(t_k)$ is a set of all pairs of locations whose pairwise distance
    is "approximately" $t_k$

## Exploratory analysis: empirical semivariogram

-   Define $N(t_k) = \{ (s_i, s_j) : \| s_i - s_j \| \in I_k \}$
-   Define the empirical semivariogram as:

$$ \gamma(t_k) = \frac{1}{2 |N(t_k)|} \sum_{(s_i, s_k) \in I_k} (Y(s_i) - Y(s_j))^2 $$

-   [What should $\gamma(t)$ look like for increasing $t$ if there is
    spatial dependence?]{.orange}
-   If spatial dependence, $Y(s_i)$ and $Y(s_j)$ should be close to each
    other if $s_i$ and $s_j$ are not far apart

## Exploratory analysis: empirical semivariogram

-   Define $N(t_k) = \{ (s_i, s_j) : \| s_i - s_j \| \in I_k \}$
-   Define the empirical semivariogram as:

$$ \gamma(t_k) = \frac{1}{2 |N(t_k)|} \sum_{(s_i, s_k) \in I_k} (Y(s_i) - Y(s_j))^2 $$

-   What should $\gamma(t)$ look like for increasing $t$ if there is
    spatial dependence?
-   If spatial dependence, $Y(s_i)$ and $Y(s_j)$ should be close to each
    other if $s_i$ and $s_j$ are not far apart
-   We then expect $\gamma(t)$ to be *increasing* with $t$: as the
    locations are farther from each other, the observations become
    increasingly different

## Exploratory analysis: empirical semivariogram

-   `install.packages("geoR")` (you may need to install XQuartz in your
    system)
-   Plot the empirical semivariogram via `geoR::variog`
-   Let's use Example 8 from earlier

```{r echo=FALSE}
set.seed(2024)
coords <- expand.grid(x <- seq(0,1,length.out=100), x)
simdf <- meshed::rmeshedgp(coords, c(3,1.5,5), block_size=200) %>% dplyr::rename(y=w)

data_plot <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()
```

```{r}
#| fig-align: center

sv <- simdf %>% with(geoR::variog(data = y, coords = cbind(Var1, Var2), messages=FALSE))

sv_df <- data.frame(dists = sv$u, variogram = sv$v, npairs = sv$n, sd = sv$sd)
sv_plot <- ggplot(sv_df, aes(x=dists, y=variogram)) + geom_point(size=2, shape=8) +
  theme_minimal()

grid.arrange(data_plot, sv_plot, nrow=1)
```

## Exploratory analysis: empirical semivariogram

-   Let's use Example 10 from earlier

```{r echo=FALSE}
set.seed(2024)
coords <- expand.grid(x <- seq(0,1,length.out=200), x)

warp1 <- meshed::rmeshedgp(coords, c(10,1.5,.02), block_size=100)
warp2 <- meshed::rmeshedgp(coords, c(2,1.5,.02), block_size=100)

newcoords <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(Var1 = Var1 + warp1$w, Var2 = Var2 + warp2$w)

w <- GpGp::fast_Gp_sim(c(5,1/30,2,0), "matern_isotropic", newcoords %>% as.matrix(), m = 20)
simdf <- warp1 %>% dplyr::select(Var1, Var2) %>% mutate(y=w)

data_plot <- ggplot(simdf, aes(Var1, Var2, fill=y)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal()

sv <- simdf %>% with(geoR::variog(data = y, coords = cbind(Var1, Var2), messages=FALSE))
sv_df <- data.frame(dists = sv$u, variogram = sv$v, npairs = sv$n, sd = sv$sd)
sv_plot <- ggplot(sv_df, aes(x=dists, y=variogram)) + geom_point(size=2, shape=8) +
  theme_minimal()

grid.arrange(data_plot, sv_plot, nrow=1)
```

## Exploratory analysis: example

-   We want to investigate a variable representing the *time it takes
    for vegetation to reach peak greenness*
-   Changes in the greenup time may be explained by climate change

```{r echo=T}
df <- read.csv("data/michigan_greenness.csv")

(data_plot <- ggplot(df, aes(Longitude, Latitude, fill=Greenness)) +
  geom_raster() +
  scale_fill_scico(palette="batlowK") + 
  theme_minimal())
```

## Exploratory analysis: example

```{r}
sv <- df %>% with(geoR::variog(data = Greenness, coords = cbind(Longitude, Latitude), messages=FALSE))
sv_df <- data.frame(dists = sv$u, variogram = sv$v, npairs = sv$n, sd = sv$sd)
sv_plot <- ggplot(sv_df, aes(x=dists, y=variogram)) + geom_point(size=2, shape=8) +
  theme_minimal()

grid.arrange(data_plot, sv_plot, nrow=1)
```
