Analiștii NOAA urmăresc în fiecare zi vaporile de fum de incendiu din imaginile din satelit și le publică ca fișiere de formă GIS deschise. Fișierele trăiesc la o adresă URL previzibilă, bazată pe dată, astfel încât să puteți trage harta de fum a zilei curente direct în R și să o reprezentați peste SUA în aproximativ 30 de rânduri. Noi folosim sf pentru a citi datele spațiale și ggplot2 să-l deseneze.
Sursa de date
Sistemul de cartografiere a pericolelor (HMS), condus de NOAA/NESDIS, este o analiză zilnică în care operatorii urmăresc fumul peste America de Nord din imaginile satelitare GOES. Fiecare penaj este un poligon cu patru atribute: satelitul, ora de început și de sfârșit și a densitate clasa – Light, Mediumsau Heavy. Shapefile-ul zilnic este postat aici:
.../Smoke_Polygons/Shapefile/YYYY/MM/hms_smokeYYYYMMDD.zip
Deoarece numele fișierului este doar data, putem construi adresa URL în mod programatic și obținem întotdeauna cea mai recentă hartă. Rețineți că HMS finalizează analiza unei anumite zile în dimineața următoare (ora de Est), așa că fișierul de astăzi poate să nu existe până atunci – ajutorul robust de la sfârșit se ocupă de asta.
Postări similare pe DataScience+:
Pachete
library(sf) library(dplyr) library(ggplot2) library(maps)
Obțineți datele
Formatăm data în URL, descarcăm și dezarhivăm shapefile într-un folder temporar și îl citim cu st_read(). HMS numește fișierele după data calendaristică, deci Sys.Date() ne oferă harta de astăzi.
day <- Sys.Date()
ymd <- format(day, "%Y%m%d")
url <- sprintf(
"https://satepsanone.nesdis.noaa.gov/pub/FIRE/web/HMS/Smoke_Polygons/Shapefile/%s/%s/hms_smoke%s.zip",
format(day, "%Y"), format(day, "%m"), ymd)
zip <- file.path(tempdir(), basename(url))
dir <- file.path(tempdir(), paste0("hms_", ymd))
download.file(url, zip, mode = "wb", quiet = TRUE)
unzip(zip, exdir = dir)
smoke <- st_read(dir, quiet = TRUE)
Inspectează ce avem
Priviți întotdeauna datele înainte de a le trasa. st_read() returnează o sf obiect — un cadru de date cu a geometry coloană — așa că instrumentele obișnuite funcționează.
smoke # note the CRS line: WGS 84 / EPSG:4326 ## Simple feature collection with 98 features and 4 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: -144.505 ymin: 12.23946 xmax: -11.87108 ymax: 85.46261 ## Geodetic CRS: WGS 84 ## First 10 features: ## Satellite Start End Density geometry ## 1 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-71.21715 19.0450... ## 2 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-71.99762 18.4664... ## 3 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-71.50895 18.6244... ## 4 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-76.61889 20.6635... ## 5 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-77.00836 21.2569... ## 6 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-80.49727 26.2337... ## 7 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-80.88491 26.4339... ## 8 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-82.60996 32.8358... ## 9 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-123.1322 42.0559... ## 10 GOES-WEST 2026198 1200 2026198 1500 Light POLYGON ((-103.8145 37.1463... table(smoke$Density) # how many plumes of each density today ## ## Heavy Light Medium ## 35 29 34
Două lucruri de observat. În primul rând, CRS este deja EPSG:4326 (longitudine/latitudine simplă), deci nu este necesară nicio reproiectare. Doilea, Density ia trei valori – Light, Medium, Heavy — care este variabila după care vom colora. The Start/End câmpurile folosesc a YYYYDDD HHMM Format UTC (ziua anului), dar nu vom avea nevoie de ele aici.
Ordonați clasele de densitate
Noi constrângem Density la un factor ordonat şi arrange() pe ea atât de grele sunt desenate dura (pe deasupra celor mai ușoare) și sunați st_make_valid() deoarece poligoane desenate manual se auto-intersectează ocazional și altfel ar rupe recolta.
smoke <- smoke %>%
filter(Density %in% c("Light", "Medium", "Heavy")) %>%
mutate(Density = factor(Density, levels = c("Light", "Medium", "Heavy"))) %>%
st_make_valid() %>%
arrange(Density)
Harta-l
Contururile de stat provin din maps pachet (CRAN pur, fără cheie API). Decupăm fumul într-o casetă de delimitare inferioară cu 48, astfel încât complotul să nu fie dominat de penaj deasupra Canadei și oceanelor, apoi stratificăm stările dedesubt și fumul deasupra, colorat după densitate.
states <- st_as_sf(map("state", plot = FALSE, fill = TRUE))
bbox <- st_bbox(c(xmin = -125, xmax = -66, ymin = 24, ymax = 50), crs = 4326)
smoke_us <- st_crop(smoke, bbox)
ggplot() +
geom_sf(data = states, fill = "grey97", color = "grey80", linewidth = 0.2) +
geom_sf(data = smoke_us, aes(fill = Density), color = NA, alpha = 0.6) +
scale_fill_manual(values = c(Light = "#FFD24D", Medium = "#FB8C00", Heavy = "#C62828")) +
coord_sf(xlim = c(-125, -66), ylim = c(24, 50), expand = FALSE) +
labs(
title = "Wildfire smoke over the U.S.",
subtitle = format(day, "NOAA HMS smoke plumes, %B %d, %Y"),
fill = "Smoke density",
caption = "Source: NOAA/NESDIS Hazard Mapping System (HMS)"
) +
theme_minimal(base_size = 13) +
theme(axis.text = element_blank(), panel.grid = element_blank(),
plot.title = element_text(face = "bold"))
Citind harta
Analiza de astăzi are 98 de pene de fum peste America de Nord, 35 dintre ele clasificate drept grele, iar 57 intersectând cele 48 inferioare. Harta transformă acel tabel în geografie: acolo unde poligoanele se stivuiesc și se întunecă, fumul este mai gros, iar penele urmăresc drumul pe care fumul a parcurs-o de la sursa incendiilor – adesea la sute de mile în aval. Reluați codul într-o zi diferită și atât numerele, cât și imaginea se schimbă; acesta este scopul construirii URL-ului de la Sys.Date().
Deoarece numărările și harta sunt generate din același obiect, paragraful de mai sus descrie întotdeauna harta pe care o vedeți – nimic nu este codificat.
O descărcare robustă
Rulați acest lucru dimineața devreme și este posibil ca fișierul de astăzi să nu fie încă postat. Acest ajutor se întoarce zi de zi până când găsește unul care există, astfel încât scenariul nu moare niciodată la o dată lipsă:
get_hms_smoke <- function(day = Sys.Date(), max_back = 3) {
for (d in seq(0, max_back)) {
date <- day - d
ymd <- format(date, "%Y%m%d")
url <- sprintf(
"https://satepsanone.nesdis.noaa.gov/pub/FIRE/web/HMS/Smoke_Polygons/Shapefile/%s/%s/hms_smoke%s.zip",
format(date, "%Y"), format(date, "%m"), ymd)
zip <- file.path(tempdir(), basename(url))
ok <- tryCatch({ download.file(url, zip, mode = "wb", quiet = TRUE); TRUE },
error = function(e) FALSE)
if (ok && file.exists(zip) && file.size(zip) > 0) {
dir <- file.path(tempdir(), paste0("hms_", ymd))
unzip(zip, exdir = dir)
message("Using HMS smoke for ", date)
return(st_read(dir, quiet = TRUE))
}
}
stop("No HMS smoke file found in the last ", max_back, " days.")
}
Avertisment: fum în sus vs. fum pe care îl respiri
HMS este un satelit produs — cartografiază fumul văzut de sus, la orice altitudine. Un penaj gros de pe hartă poate sta la cinci kilometri în sus și poate lăsa aerul limpede la nivelul solului, așa că HMS este o măsură de fum transportnu de ceea ce respiră oamenii. Pentru calitatea aerului de suprafață, asociați-l cu PM2,5 la sol de la EPA/USFS AirNow Fire and Smoke Map, care are și un API – o analiză naturală de urmărire.
Asta e tot. Acum aveți un script care extrage un produs satelit analizat manual și îl transformă într-o hartă națională și se rulează din nou pentru orice zi doriți.
Acest articol a fost publicat pentru prima dată pe DataScience+, o comunitate de autori de tutoriale R și Python. Aveți o tehnică de știință a datelor care merită împărtășită? Scrieți pentru noi – nu este necesară prezentarea.
