---
title: "Environmental Data Analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Environmental Data Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r not-cran-check, include = FALSE}
# Skip code evaluation on CRAN (external API calls will fail)
not_cran_check <- identical(Sys.getenv("NOT_CRAN"), "true") &&
  !nzchar(Sys.getenv("_R_CHECK_PACKAGE_NAME_"))
run_vignette <- isTRUE(getOption("pkgdown.in_progress")) || not_cran_check
knitr::opts_chunk$set(eval = run_vignette, purl = run_vignette)
# Inline `r` code runs even when chunks are skipped (eval = FALSE), so
# values computed in chunks go through this helper. `expr` is only
# evaluated when the chunks ran, and falls back if it cannot be computed.
inline_value <- function(expr, fallback = "n/a") {
  if (!run_vignette) {
    return(fallback)
  }
  tryCatch(expr, error = function(e) fallback)
}
```

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,        # standardized width that fits vignette text column
  fig.height = 5,       # default height, override in individual chunks as needed
  fig.dpi = 150,        # higher resolution for crisp rendering
  out.width = "100%",   # constrain to page width, prevents overflow
  warning = FALSE,
  message = FALSE
)
```

```{r setup}
library(mongolstats)
library(mongolmaps)
library(dplyr)
library(ggplot2)
library(tidyr)
library(stringr)

old_options <- nso_options(mongolstats.lang = "en")

# Global theme with proper margins to prevent text cutoff
old_theme <- theme_set(
  theme_minimal(base_size = 11) +
    theme(
      plot.margin = margin(10, 10, 10, 10),
      plot.title = element_text(size = 13, face = "bold"),
      plot.subtitle = element_text(size = 10, color = "grey40"),
      legend.text = element_text(size = 9),
      legend.title = element_text(size = 10)
    )
)
```

## Overview

Air pollution is a leading environmental health risk in Mongolia, above all in Ulaanbaatar in winter. NSO publishes measurements from air quality monitoring stations, river and lake water quality stations, and weather stations that count dust days. Station data describe the air or water at one place; they are a starting point for exposure assessment, not a measure of what people breathe or drink.

This guide shows how to:

- Track fine particulate matter (PM2.5) and sulphur dioxide (SO₂) at Ulaanbaatar stations
- Compare annual concentrations with Mongolia's air quality standards
- Look at river and lake water quality
- Map dust days across the country

## Air Quality Monitoring

> **About the air quality tables**
>
> - **Annual table** (`DT_NSO_2024_135V01`): annual mean concentration by station and pollutant, selected by `Year`.
> - **Monthly station tables** (`DT_NSO_2400_015V2` for NO₂, `DT_NSO_2400_015V5` for PM2.5): for each station and month, the monthly average, the monthly maximum, the number of measurement days ("MN by day"), and the percentage of days above the air quality standard. Month codes count back from the latest month, so select months by their `"YYYY-MM"` labels, listed with `nso_table_periods()`.
>
> The NSO catalogue reuses the ids `DT_NSO_2400_015V1` (SO₂), `V3` (CO) and `V4` (PM10) for unrelated tables in another folder, so those monthly series cannot be fetched by id alone (mongolstats warns about the ambiguity). This vignette uses a snapshot bundled with mongolstats for monthly SO₂.
>
> NSO reports concentrations in mg/m³. We convert them to µg/m³, the unit used by the WHO guidelines (1 mg/m³ = 1,000 µg/m³). Annual means are rounded to 0.01 mg/m³ (10 µg/m³), which is coarse next to standards of 20–40 µg/m³.

### Most Polluted Stations

We rank Ulaanbaatar stations by their 2024 annual mean PM2.5. We use 2024 rather than the latest year because several 2025 annual values are exact zeros at stations that reported high levels the year before, which looks like missing data.

```{r air-quality-data}
air_annual <- tryCatch(
  nso_data("DT_NSO_2024_135V01", selections = list(Year = "2024"), labels = "en"),
  error = function(e) NULL
)

# Fall back to the bundled snapshot if the API is unavailable
if (is.null(air_annual) || nrow(air_annual) == 0) {
  air_annual <- utils::read.csv(
    system.file("extdata", "air_annual_cached.csv", package = "mongolstats"),
    check.names = FALSE,
    stringsAsFactors = FALSE
  )
}

air_annual <- air_annual |>
  filter(!is.na(value)) |>
  mutate(
    Station = str_trim(`Station location_en`),
    Pollutant = str_trim(Indicator_en),
    ug_m3 = value * 1000 # mg/m3 to ug/m3
  )

# Ulaanbaatar stations
ub_stations <- c(
  "Misheel-Expo center", "West crossroad", "1st micro district",
  "13th micro district", "32nd Toirog", "Ofitseruudiin ordon",
  "Kharkhorin market", "Urgakh naran micro district", "Dambdarjaa",
  "Khailaast", "Tolgoit", "Zuragt", "Amgalan", "Nisekh",
  "Tavan buudal", "Bayankhoshuu", "Sharkhad", "100 ail"
)

# Up to 10 Ulaanbaatar stations with the highest annual PM2.5
top_stations <- air_annual |>
  filter(Station %in% ub_stations, str_detect(Pollutant, "PM2.5")) |>
  slice_max(ug_m3, n = 10, with_ties = FALSE) |>
  pull(Station)

top_stations
```

Only `r inline_value(length(top_stations))` Ulaanbaatar stations reported an annual PM2.5 mean for 2024, so all of them are included.

### PM2.5: Monthly Trends

Fine particulate matter (PM2.5) is the air pollutant with the strongest evidence of harm: long-term exposure raises the risk of heart and lung disease, lung cancer, and death. We keep only months with measurements on at least 75% of days, a common completeness rule for a valid monthly mean.

```{r pm25-data}
# Labels of the months from `from` onward
recent_months <- function(tbl, from = "2021-01-01") {
  periods <- nso_table_periods(tbl)
  periods$label[!is.na(periods$date) & periods$date >= as.Date(from)]
}

pm25_monthly <- tryCatch(
  nso_data(
    "DT_NSO_2400_015V5",
    selections = list(Month = recent_months("DT_NSO_2400_015V5")),
    labels = "en"
  ),
  error = function(e) NULL
)

# Bundled snapshot: fallback for PM2.5 and the only source of monthly SO2
air_cached <- utils::read.csv(
  # gzip-compressed to keep the installed package small;
  # read.csv() decompresses transparently
  system.file("extdata", "air_monthly_cached.csv.gz", package = "mongolstats"),
  check.names = FALSE,
  stringsAsFactors = FALSE
)

if (is.null(pm25_monthly) || nrow(pm25_monthly) == 0) {
  pm25_monthly <- filter(air_cached, str_detect(Indicator_en, "PM2.5"))
}

# One row per station and month, with one column per measure
reshape_monthly <- function(data) {
  data |>
    mutate(
      Station = str_trim(Location_en),
      measure = str_extract(Indicator_en, "(Average|Maximum|MN by day|Exceeding percent of AQS)$"),
      date = as.Date(paste0(Month_en, "-01"))
    ) |>
    filter(Station %in% top_stations, !is.na(measure)) |>
    select(Station, date, measure, value) |>
    pivot_wider(names_from = measure, values_from = value) |>
    rename(
      average = Average, maximum = Maximum,
      days_measured = `MN by day`, pct_days_above = `Exceeding percent of AQS`
    ) |>
    mutate(average = average * 1000, maximum = maximum * 1000) # mg/m3 to ug/m3
}

pm25 <- reshape_monthly(pm25_monthly) |>
  filter(date >= as.Date("2021-01-01")) |>
  mutate(complete = days_measured >= 0.75 * lubridate::days_in_month(date))

pm25_complete <- filter(pm25, complete, !is.na(average))

# For line charts: every station-month, with excluded months as gaps (NA)
# so lines break instead of joining across missing data
pm25_lines <- pm25 |>
  complete(Station, date = seq(min(date), max(date), by = "month")) |>
  mutate(average = if_else(complete %in% TRUE, average, NA_real_))
```

Of `r inline_value(nrow(pm25))` station-months, `r inline_value(nrow(pm25_complete))` meet the completeness rule.

```{r pm25-trends, fig.height=8, fig.alt="Small multiples of monthly mean PM2.5 at the most polluted Ulaanbaatar stations since 2021, with reference lines at Mongolia's annual standard of 25 micrograms per cubic metre and the WHO guideline of 5. Every station shows large winter peaks."}
pm25_lines |>
  ggplot(aes(x = date, y = average)) +
  geom_hline(yintercept = 25, linetype = "dashed", color = "darkred") + # Mongolia annual standard
  geom_hline(yintercept = 5, linetype = "dotted", color = "grey30") + # WHO 2021 annual guideline
  geom_line(color = "firebrick", linewidth = 0.6) +
  geom_point(color = "firebrick", size = 0.8) +
  facet_wrap(~Station, ncol = 2) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(limits = c(0, NA)) +
  labs(
    title = "Monthly Mean PM2.5 at Ulaanbaatar's Most Polluted Stations",
    subtitle = "Dashed = Mongolia annual standard (25 µg/m³); dotted = WHO annual guideline (5 µg/m³)",
    x = NULL,
    y = "PM2.5 (µg/m³)",
    caption = "Months with measurements on at least 75% of days. Source: NSO Mongolia, DT_NSO_2400_015V5"
  ) +
  theme(strip.text = element_text(face = "bold"), panel.grid.minor = element_blank())
```

The reference lines are for **annual** means, so they give scale rather than a test: compliance is judged on the annual mean (next section), and a monthly mean above the line marks a polluted month, not a legal exceedance. Mongolia's annual PM2.5 standard (MNS 4585:2016) is 25 µg/m³; the WHO's 2021 guideline is 5 µg/m³. The monthly maximum that NSO also publishes is left out because its definition (highest daily mean or highest single reading) is not documented.

### PM2.5: Days Above the Standard

The monthly tables also report the percentage of measured days on which PM2.5 exceeded Mongolia's air quality standard. Unlike a monthly mean, this compares like with like: daily values against a daily standard.

```{r pm25-exceedance, fig.height=4.5, fig.alt="Heat map of the percentage of days above the PM2.5 standard by station and month. Winter months are dark red at every station, and summer months are near zero."}
pm25_complete |>
  ggplot(aes(x = date, y = Station, fill = pct_days_above)) +
  geom_tile(color = "white", linewidth = 0.2) +
  scale_fill_gradient(
    low = "#fbe3e3", high = "#8b1a1a", na.value = "grey90",
    limits = c(0, 100), name = "% of days\nabove standard"
  ) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y", expand = c(0, 0)) +
  labs(
    title = "Share of Days with PM2.5 Above Mongolia's Standard",
    subtitle = "Each tile is one station-month; grey = not reported; blank = fewer than 75% of days measured",
    x = NULL, y = NULL,
    caption = "Source: NSO Mongolia, DT_NSO_2400_015V5"
  ) +
  theme(panel.grid = element_blank())
```

```{r pm25-facts, include = FALSE}
season_of <- function(date) {
  m <- as.integer(format(date, "%m"))
  case_when(m %in% c(12, 1, 2) ~ "Dec-Feb", m %in% 6:8 ~ "Jun-Aug", .default = "Other")
}
pm25_season <- if (run_vignette) {
  pm25_complete |>
    mutate(season = season_of(date)) |>
    summarise(
      median_pct_above = median(pct_days_above, na.rm = TRUE),
      mean_ug = mean(average),
      .by = season
    )
}
season_value <- function(s, col) pm25_season[[col]][pm25_season$season == s]
```

> **Key findings:** In December–February, the median station-month had **`r inline_value(round(season_value("Dec-Feb", "median_pct_above")))`%** of measured days above the PM2.5 standard, and a mean of `r inline_value(round(season_value("Dec-Feb", "mean_ug")))` µg/m³. In June–August, the median was `r inline_value(round(season_value("Jun-Aug", "median_pct_above")))`% of days, with a mean of `r inline_value(round(season_value("Jun-Aug", "mean_ug")))` µg/m³. Exposure is concentrated in the heating season, when household coal and wood stoves, heating plants and temperature inversions combine.

### Sulphur Dioxide (SO₂)

SO₂ comes mainly from burning coal, so it tracks fuel use for heating. Monthly SO₂ comes from the bundled snapshot (see above). Its "measurement days" column counts readings rather than days, so we do not apply the completeness rule here.

```{r so2-trends, fig.height=8, fig.alt="Small multiples of monthly mean SO2 at the most polluted Ulaanbaatar stations from the bundled snapshot, with a reference line at Mongolia's annual standard of 20 micrograms per cubic metre. Concentrations rise sharply each winter."}
so2 <- air_cached |>
  filter(str_detect(Indicator_en, "Sulphur dioxide")) |>
  reshape_monthly() |>
  filter(!is.na(average))

so2 |>
  ggplot(aes(x = date, y = average)) +
  geom_hline(yintercept = 20, linetype = "dashed", color = "darkorange3") + # Mongolia annual standard
  geom_line(color = "darkorange", linewidth = 0.6) +
  geom_point(color = "darkorange", size = 0.8) +
  facet_wrap(~Station, ncol = 2) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(limits = c(0, NA)) +
  labs(
    title = "Monthly Mean SO₂ at Ulaanbaatar's Most Polluted Stations",
    subtitle = "Bundled snapshot; dashed = Mongolia annual standard (20 µg/m³)",
    x = NULL,
    y = "SO₂ (µg/m³)",
    caption = "Source: NSO Mongolia (snapshot of DT_NSO_2400_015V1)"
  ) +
  theme(strip.text = element_text(face = "bold"), panel.grid.minor = element_blank())
```

```{r so2-facts, include = FALSE}
so2_ratio <- if (run_vignette) {
  so2 |>
    mutate(season = season_of(date)) |>
    summarise(
      winter = mean(average[season == "Dec-Feb"]),
      summer = mean(average[season == "Jun-Aug"]),
      .by = Station
    ) |>
    mutate(ratio = winter / summer)
}
```

> **Seasonal pattern:** Mean SO₂ in December–February is `r inline_value(round(min(so2_ratio$ratio)))` to `r inline_value(round(max(so2_ratio$ratio)))` times the June–August mean, depending on the station (median `r inline_value(round(median(so2_ratio$ratio)))` times). This swing, larger than for PM2.5, points to coal burning for heating as the main winter source.

### Annual Means Against the Standards

The exceedance factor is the annual mean divided by Mongolia's annual standard (MNS 4585:2016). A factor above 1 means the station did not meet the standard.

```{r compliance-check, fig.alt="Bar chart of exceedance factors (annual mean divided by Mongolia's annual standard) for PM2.5, SO2 and NO2 at the most polluted Ulaanbaatar stations in 2024"}
# Mongolia's annual standards (MNS 4585:2016), ug/m3
annual_standard <- c("PM2.5" = 25, "SO2" = 20, "NO2" = 40)

compliance <- air_annual |>
  filter(Station %in% top_stations) |>
  mutate(
    Pollutant_Type = case_when(
      str_detect(Pollutant, "PM2.5") ~ "PM2.5",
      str_detect(Pollutant, "Sulphur dioxide") ~ "SO2",
      str_detect(Pollutant, "Nitrogen dioxide") ~ "NO2",
      .default = NA_character_
    )
  ) |>
  filter(!is.na(Pollutant_Type)) |>
  mutate(
    Standard = annual_standard[Pollutant_Type],
    Exceedance_Factor = ug_m3 / Standard,
    Status = if_else(Exceedance_Factor > 1, "Above standard", "Meets standard")
  )

compliance |>
  ggplot(aes(x = Exceedance_Factor, y = reorder(Station, Exceedance_Factor), fill = Status)) +
  geom_col() +
  geom_vline(xintercept = 1, linetype = "dashed", color = "black") +
  geom_text(aes(label = round(Exceedance_Factor, 1)), hjust = -0.15, size = 2.8) +
  facet_wrap(~Pollutant_Type, scales = "free_x") +
  scale_fill_manual(values = c("Meets standard" = "#2a78d6", "Above standard" = "#c0392b"), name = NULL) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.2))) +
  labs(
    title = "Annual Mean Divided by Mongolia's Annual Standard, 2024",
    subtitle = "Values above 1 (dashed line) do not meet the standard",
    x = "Exceedance factor", y = NULL,
    caption = "Annual means are rounded to 10 µg/m³. Source: NSO Mongolia, DT_NSO_2024_135V01"
  ) +
  theme(legend.position = "bottom", panel.grid.major.y = element_blank())
```

```{r compliance-facts, include = FALSE}
compliance_summary <- if (run_vignette) {
  compliance |>
    summarise(n = n(), above = sum(Status == "Above standard"), max_factor = max(Exceedance_Factor), .by = Pollutant_Type)
}
cs <- function(p, col) compliance_summary[[col]][compliance_summary$Pollutant_Type == p]
```

> **Results:** In 2024, `r inline_value(cs("PM2.5", "above"))` of `r inline_value(cs("PM2.5", "n"))` stations exceeded the PM2.5 standard (up to `r inline_value(round(cs("PM2.5", "max_factor"), 1))` times), `r inline_value(cs("SO2", "above"))` of `r inline_value(cs("SO2", "n"))` exceeded the SO₂ standard (up to `r inline_value(round(cs("SO2", "max_factor"), 1))` times), and `r inline_value(cs("NO2", "above"))` of `r inline_value(cs("NO2", "n"))` exceeded the NO₂ standard. Because annual means are rounded to 10 µg/m³, a factor near 1 cannot separate compliance from exceedance: a reported 20 µg/m³ of PM2.5 could be anywhere from 15 to 25. All stations remain far above the WHO's 2021 annual PM2.5 guideline of 5 µg/m³.

## Water Quality Monitoring

NSO publishes river and lake water quality by station. We look first at electrical conductivity, a measure of dissolved minerals and salts.

```{r water-quality, fig.alt="Bar chart of the 15 water quality stations with the highest electrical conductivity"}
water <- nso_data("DT_NSO_2300_005V12", selections = list(), labels = "en") |>
  filter(!is.na(value)) |>
  mutate(
    Station = str_trim(`Station location_en`),
    Indicator = str_trim(Indicator_en)
  )
water_year <- max(water$Year_en)

conductivity <- water |>
  filter(str_detect(Indicator, "Electrical conductivity")) |>
  slice_max(order_by = value, n = 15)

conductivity |>
  ggplot(aes(x = value, y = reorder(Station, value))) +
  geom_col(fill = "#2a78d6") +
  geom_text(aes(label = scales::label_comma()(value)), hjust = -0.1, size = 3) +
  scale_x_continuous(labels = scales::label_comma(), expand = expansion(mult = c(0, 0.15))) +
  labs(
    title = paste0("Top 15 Water Stations by Electrical Conductivity, ", water_year),
    subtitle = "Higher values mean more dissolved minerals and salts",
    x = "Conductivity (µS/cm)", y = NULL,
    caption = "Source: NSO Mongolia, DT_NSO_2300_005V12"
  ) +
  theme(panel.grid.major.y = element_blank(), panel.grid.minor = element_blank())
```

The highest values, at `r inline_value(knitr::combine_words(str_extract(conductivity$Station[1:3], "^[^,]+")))`, are lakes. Lakes in closed basins with no outflow, such as Uvs and Khyargas in western Mongolia, are naturally salty, because water leaves by evaporation and the salts stay behind. Conductivity has no health-based limit and does not by itself indicate pollution. To look for pollution, use indicators of organic waste and sewage, such as biochemical oxygen demand (BOD₅) and ammonium nitrogen (NH₄-N).

### The Tuul River Through Ulaanbaatar

```{r water-tuul}
tuul <- water |>
  filter(str_starts(Station, "Tuul")) |>
  mutate(Indicator = case_when(
    str_detect(Indicator, "BOD5") ~ "BOD5 (mg/L)",
    str_detect(Indicator, "NH4-N") ~ "NH4-N (mg/L)",
    str_detect(Indicator, "Water quality index") ~ "Water quality index",
    .default = NA_character_
  )) |>
  filter(!is.na(Indicator)) |>
  select(Station, Indicator, value) |>
  pivot_wider(names_from = Indicator, values_from = value) |>
  arrange(`NH4-N (mg/L)`)

knitr::kable(tuul, caption = paste0("Tuul River stations, ", water_year))
```

```{r tuul-facts, include = FALSE}
tuul_value <- function(station, col) tuul[[col]][tuul$Station == station]
```

At **Zaisan**, in the south of Ulaanbaatar, the Tuul has a BOD₅ of `r inline_value(tuul_value("Tuul, Zaisan", "BOD5 (mg/L)"))` mg/L and ammonium of `r inline_value(tuul_value("Tuul, Zaisan", "NH4-N (mg/L)"))` mg/L. At **Songino**, downstream and west of the city, these rise to `r inline_value(tuul_value("Tuul, Songino", "BOD5 (mg/L)"))` and `r inline_value(tuul_value("Tuul, Songino", "NH4-N (mg/L)"))` mg/L. This gradient is consistent with urban wastewater entering the river as it passes the city, and it matters for anyone using river water downstream for livestock, irrigation or recreation. A single year of measurements per station cannot show trends; check other years before drawing conclusions.

## Dust Days

Dust storms carry coarse particles (PM10) and are linked to respiratory and cardiovascular illness. Dust is most frequent in the Gobi and steppe regions in spring. Dust days vary a lot from year to year, so we average over 10 years rather than ranking a single year. NSO leaves a cell blank when a station had no dust days in a month, and a few months are missing from the table altogether (mostly summer months); we count both as zero.

```{r dust-days, fig.alt="Bar chart of the 10 weather stations with the most dust days per year on average over the last 10 complete years, with lines showing the lowest and highest year"}
dust <- nso_data("DT_NSO_2400_028V1", selections = list(), labels = "en") |>
  mutate(
    Station = str_trim(`Station location_en`),
    year = as.integer(substr(Month_en, 1, 4))
  )

# The last 10 calendar years before the current, partial year
dust_years <- sort(unique(dust$year)) |>
  head(-1) |>
  tail(10)

# Blank cells, and the few months missing from the table altogether, are
# taken to mean no dust days and count as zero
dust_annual <- dust |>
  filter(year %in% dust_years) |>
  summarise(dust_days = sum(value, na.rm = TRUE), .by = c(Station, year))

dust_summary <- dust_annual |>
  summarise(
    mean_days = mean(dust_days),
    min_days = min(dust_days),
    max_days = max(dust_days),
    .by = Station
  ) |>
  arrange(desc(mean_days))

dust_summary |>
  slice_max(mean_days, n = 10) |>
  ggplot(aes(x = mean_days, y = reorder(Station, mean_days))) +
  geom_errorbar(aes(xmin = min_days, xmax = max_days), width = 0.3, color = "grey60", orientation = "y") +
  geom_point(color = "#b8860b", size = 3) +
  labs(
    title = paste0("Stations with the Most Dust Days, ", min(dust_years), "–", max(dust_years)),
    subtitle = "Average dust days per year; lines show the lowest and highest year",
    x = "Dust days per year", y = NULL,
    caption = "Source: NSO Mongolia, DT_NSO_2400_028V1"
  ) +
  theme(panel.grid.major.y = element_blank(), panel.grid.minor = element_blank())
```

The three dustiest stations, `r inline_value(knitr::combine_words(dust_summary$Station[1:3]))`, average `r inline_value(round(min(dust_summary$mean_days[1:3])))` to `r inline_value(round(max(dust_summary$mean_days[1:3])))` dust days a year. The range bars show how much a single year can mislead: the same station can have a handful of dust days one year and several dozen the next.

### Mapping Dust Days

The weather stations are in towns, so mongolmaps' settlement points can place them on a map. Station names are matched to the capital and aimag centres first, then to soum centres where the name is unique.

```{r dust-map, fig.alt="Map of Mongolia with weather stations shown as points sized and coloured by average dust days per year; the largest values are in the Gobi in the south"}
settlements <- mn_settlements() |>
  mutate(rank = match(type, c("capital", "aimag_centre", "soum_centre")))

# Khentii's aimag centre, Undurkhaan, was renamed Chinggis city in 2013
dust_points <- dust_summary |>
  mutate(name = if_else(Station == "Undurkhaan", "Chingis", Station)) |>
  inner_join(settlements, by = join_by(name), relationship = "one-to-many") |>
  # Keep the highest-ranked match; drop names shared by several soum centres
  filter(rank == min(rank), .by = Station) |>
  filter(n() == 1, .by = Station) |>
  sf::st_as_sf()

unplaced <- setdiff(dust_summary$Station, dust_points$Station)

mn_map(mn_aimags(), title = "Average Dust Days per Year") +
  geom_sf(data = dust_points, aes(size = mean_days, color = mean_days), alpha = 0.85) +
  scale_color_gradient(low = "#f3d9a4", high = "#8a5a00", name = "Dust days\nper year") +
  scale_size_area(max_size = 8, guide = "none") +
  labs(subtitle = paste0(min(dust_years), "–", max(dust_years), ". Data: NSO Mongolia, DT_NSO_2400_028V1"))
```

`r inline_value(length(unplaced))` station(s) could not be placed from settlement names (`r inline_value(knitr::combine_words(unplaced))`). For exposure studies, use the weather stations' own coordinates from the National Agency for Meteorology and Environmental Monitoring.

## Next Steps

- **Ulaanbaatar districts**: See the [UB District Analysis](ub-health-environment.html) for the monitoring network by district
- **Discover More Tables**: See the [Discovery Guide](discovery.html) for search techniques
- **Mapping**: Learn spatial analysis in the [Mapping Guide](mapping.html)
- **Reference**: Browse all functions in the [Reference](https://temuulene.github.io/mongolstats/reference/index.html)

## Key Environmental Tables

```{r env-tables, echo=FALSE}
tibble::tribble(
  ~Indicator, ~Table_ID, ~Geographic_Level,
  "Annual mean air pollutant concentrations", "DT_NSO_2024_135V01", "Station",
  "Monthly NO2 by station", "DT_NSO_2400_015V2", "Station",
  "Monthly PM2.5 by station", "DT_NSO_2400_015V5", "Station",
  "Water quality", "DT_NSO_2300_005V12", "Station",
  "Dust days", "DT_NSO_2400_028V1", "Station",
  "Forest fires", "DT_NSO_2400_005V1", "Aimag"
) |>
  knitr::kable()
```

```{r restore-session, include=FALSE}
options(old_options)
ggplot2::theme_set(old_theme)
```