---
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")

# Alternative to: geobr::read_municipality(year = 2024)
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, which is required for distance-based
spatial analysis.

```{r setup, warning = FALSE, message = FALSE}
rain_sf <- floods_rs |>
  st_as_sf(coords = c("lon", "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}

data("inmet_stations")

mun_stations2 <- mun_stations |> 
  left_join(inmet_stations |> select(code_wmo, creation_year)) |> 
  filter(creation_year <= 2024) |> 
  arrange(distance) |> 
  distinct(code_ibge7, .keep_all = TRUE)

mun_floods <- floods_rs |>
  left_join(
    mun_stations2,
    by = c("code_wmo" = "code_wmo")
  )

p2 <- mun_24 |>
  left_join(
    mun_floods |>
      select(code_ibge7, total_rainfall),
    by = c("code_muni" = "code_ibge7")
  ) |>
  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 advantages and limitations:

-   **Ordinary kriging** accounts for spatial dependence and usually produces
    more realistic spatial patterns. However, as a statistical interpolation
    method, it can produce physically impossible values (such as negative
    rainfall) and tends to smooth the surface, so very high rainfall totals
    may be underestimated.
-   **Nearest station** preserves the observed values, including extremes.
    However, many municipalities may share the same station, which creates
    large areas with identical values and artificial boundaries between
    neighboring municipalities. Real rainfall does not change abruptly at
    these boundaries.

The histograms below compare the distribution of municipality-level
estimates from both approaches.

```{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()
)

```