---
title: "Spatial Interpolation Using Ordinary Kriging"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Spatial Interpolation Using Ordinary Kriging}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Introduction

Spatial interpolation is commonly used to estimate meteorological variables
at locations where no direct observations are available. The
`kriging_inmet()` function implements **ordinary kriging** using weather
station observations and spatial prediction locations provided as `sf`
objects.

This vignette demonstrates how to interpolate accumulated rainfall from
INMET weather stations.

```{r setup packages, warning = FALSE, message = FALSE}
library(climateBR)
library(dplyr)
library(sf)
library(ggplot2)
library(patchwork)

```

# Input Data

The example uses rainfall observations from the 2024 Rio Grande do Sul
flood event together with municipality geometries.

```{r import data, warning = FALSE, message = FALSE}
data("floods_rs")
data("mun_stations_distance")

tmp <- tempfile(fileext = ".rda")
on.exit(unlink(tmp), add = TRUE)

download.file(
  "https://github.com/kaiorb52/dados_municipais/raw/main/mun_24.rda",
  destfile = tmp,
  mode = "wb"
)

load(tmp)

```

The station observations are converted to an `sf` object and projected to
a planar coordinate reference system suitable for distance-based spatial
analysis.

```{r setup, warning = FALSE, message = FALSE}
rain_sf <- floods_rs |>
  st_as_sf(coords = c("long", "lat"), crs = 4326) |>
  st_transform(crs = 29193)

mun_centroid <- mun_24 |>
  select(code_muni, geom) |>
  mutate(
    ponto = st_point_on_surface(geom)
  ) |>
  st_drop_geometry()

mun_grid_sf <- mun_centroid |>
  st_as_sf() |>
  st_transform(crs = 29193)

```

# Running the Kriging Model

The `kriging_inmet()` function computes an empirical variogram, fits a
spherical variogram model, and performs ordinary kriging predictions for
the target locations.

```{r, warning = FALSE, message = FALSE}
krig_result <- kriging_inmet(
  stations_df = rain_sf,
  mun_geo     = mun_grid_sf
)
```

The predicted values (`var1.pred`) and prediction variances (`var1.var`)
are then attached to the municipality dataset.

```{r, warning = FALSE, message = FALSE}
mun_pred <- mun_centroid |>
  select(code_muni) |>
  bind_cols(
    krig_result |>
      st_drop_geometry() |>
      select(
        precip_krig     = var1.pred,
        precip_krig_var = var1.var
      )
  )
```

# Visualizing Kriging Predictions

The following map displays municipality-level rainfall estimates obtained
through ordinary kriging.

```{r, warning = FALSE, message = FALSE}
p1 <- mun_24 |>
  left_join(mun_pred, by = "code_muni") |>
  ggplot() +
  geom_sf(aes(fill = precip_krig), color = NA) +
  labs(title = "Ordinary Kriging") +
  scale_fill_distiller(
    palette = "RdYlGn",
    direction = 1,
    na.value = "grey80"
  ) +
  theme_void() +
  theme(
    legend.position = c(0.9, 0.1)
  )
```

# Comparison with the Nearest-Station Approach

For comparison, the map below assigns each municipality the rainfall value
from its closest weather station.

```{r, warning = FALSE, message = FALSE}
mun_floods <- floods_rs |>
  left_join(
    mun_stations_distance |>
      filter(ano == 2024, i == 1),
    by = c("id_who" = "codigo_wmo")
  )

p2 <- mun_24 |>
  left_join(
    mun_floods |>
      select(code_muni, total_rainfall),
    by = "code_muni"
  ) |>
  ggplot() +
  geom_sf(aes(fill = total_rainfall), color = NA) +
  labs(title = "Nearest Station") +
  scale_fill_distiller(
    palette = "RdYlGn",
    direction = 1,
    na.value = "grey80"
  ) +
  theme_void() +
  theme(
    legend.position = c(0.9, 0.1)
  )
```

The figure below compares rainfall estimates generated using ordinary
kriging against values obtained from the nearest-station approach.

```{r, fig.width=8, fig.height=6, echo=FALSE, warning = FALSE, message = FALSE}

p1 + p2

```

Ordinary kriging generally produces smoother spatial surfaces and
incorporates information from multiple nearby stations, whereas the
nearest-station method assigns the same value to all municipalities linked
to a given station and may introduce abrupt spatial discontinuities.

Each approach has important advantages and limitations. Ordinary kriging 
accounts for spatial dependence and usually generates more realistic spatial 
patterns than others methods. However, because it is a statistical 
interpolation technique, it may produce physically unrealistic estimates, 
such as negative rainfall values, and it often smooths the spatial field, 
reducing the magnitude of extreme precipitation events. As a result, 
very high observed rainfall totals may be underestimated in the interpolated 
surface. In contrast, the nearest-station approach preserves the original 
observed values, including extremes, but can create large areas with identical
rainfall estimates because multiple municipalities may be assigned to the same
station. This can lead to artificial boundaries and abrupt changes between 
neighboring municipalities that do not reflect the continuous nature of 
precipitation processes.

```{r, fig.width=9, fig.height=5, echo=FALSE, warning = FALSE, message = FALSE}
(
mun_pred |> 
  ggplot(aes(x = precip_krig)) +
  geom_histogram(color = "black", fill = "steelblue4", alpha = 0.85) +
  coord_cartesian(xlim = c(0, 650), ylim = c(0, 3500)) +
  labs(title = "Ordinary kriging") +
  theme_minimal() +
mun_floods |> 
  #filter(total_rainfall >= 1) |> 
  ggplot(aes(x = total_rainfall)) +
  geom_histogram(color = "black", fill = "steelblue4", alpha = 0.85) +
  coord_cartesian(xlim = c(0, 650), ylim = c(0, 3500)) +
  labs(title = "Nearest Station") +
  theme_minimal()
)

```