---
title: "Finding the Closest INMET Station for Each Municipality"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Finding the Closest INMET Station for Each Municipality}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Introduction

Many climate analyses need to assign each municipality to its closest weather station. `nearest_stations()` automates this by finding the nearest INMET station(s) for every municipality.

```{r setup packages, warning = FALSE, message = FALSE}

library(climateBR)
library(dplyr)
library(ggplot2)
library(sf)

```

# Preparing the data

The package includes datasets of Brazilian municipalities (`municipality`) and INMET weather stations (`inmet_stations`). This example uses the state of Pernambuco (PE).

```{r}
data("municipality")
data("inmet_stations")

municipality_pe <- municipality |>
  filter(
    state_muni == "PE",
    code_ibge7 != "2605459" # Fernando de Noronha
  )

inmet_stations_pe <- inmet_stations |>
  filter(state_station == "PE")
```

# Finding the nearest station

Use `nearest_stations()` with `n = 1` to find the closest station for each municipality.

```{r}
mun_stations_pe <- nearest_stations(
  municipality = municipality_pe,
  inmet_stations = inmet_stations_pe,
  n = 1
)

head(mun_stations_pe)
```

# Municipality–station connections

The top map colors each municipality by its nearest station. The bottom map draws a line from each municipality centroid (×) to its nearest station (◆).

```{r, echo=FALSE, warning = FALSE, message = FALSE}
tmp <- tempfile(fileext = ".rda")
on.exit(unlink(tmp), add = TRUE)

# Alternative to: geobr::read_municipality(year = 2024)
download.file(
  "https://github.com/kaiorb52/dados_municipais/raw/main/mun_24.rda",
  destfile = tmp,
  mode = "wb"
)

load(tmp)

mun_24_pe <- mun_24 |>
  filter(abbrev_state == "PE") |>
  select(code_muni, geom)

mun_stations_pe2 <- mun_24_pe |>
  left_join(
    mun_stations_pe,
    by = c("code_muni" = "code_ibge7")
  )
```

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

p1 <- ggplot() +
  geom_sf(
    data = mun_stations_pe2,
    aes(fill = code_wmo),
    color = "grey20"
  ) +
  coord_sf(
    xlim = c(-41.25, -35),
    ylim = c(-10, -6.5)
  ) +
  geom_text(
    data = inmet_stations_pe,
    aes(lon, lat, label = code_wmo),
    shape = 18,
    color = "black",
    size = 3.25,
    fontface = "bold"
  ) +
  theme_void() +
  guides(fill = "none") + 
  scale_fill_manual(
    values = c(
      "steelblue",
      "firebrick",
      "seagreen",
      "orange",
      "mediumpurple",
      "grey",
      "white",
      "darkgreen",
      "antiquewhite",
      "aquamarine",
      "brown",
      "coral",
      "royalblue"
    )
  )

distance_lines <- mun_stations_pe |>
  left_join(
    inmet_stations_pe |>
      select(
        code_wmo,
        station_lat = lat,
        station_lon = lon
      ),
    by = "code_wmo"
  ) |>
  left_join(
    municipality_pe |>
      select(
        code_ibge7,
        mun_lat = lat,
        mun_lon = lon
      ),
    by = "code_ibge7"
  ) |>
  rowwise() |>
  mutate(
    geometry = st_sfc(
      st_linestring(
        matrix(
          c(
            mun_lon, mun_lat,
            station_lon, station_lat
          ),
          ncol = 2,
          byrow = TRUE
        )
      ),
      crs = 4326
    )
  ) |>
  st_as_sf()

p2 <- ggplot() +
  geom_sf(data = mun_24_pe, fill = "grey90", color="grey25") +
  geom_sf(data = distance_lines, linewidth = 0.12) +
  geom_point(
    data = municipality_pe,
    aes(lon, lat),
    shape = 4,
    color = "black",
    size = 0.4
  ) +
  geom_point(
    data = inmet_stations_pe,
    aes(lon, lat),
    shape = 18,
    color = "firebrick2",
    size = 3
  ) +
  coord_sf(
    xlim = c(-41.25, -35),
    ylim = c(-10, -6.5)
  ) +
  theme_void()

library(patchwork)
(p1 / p2)

```