Skip to content

Instantly share code, notes, and snippets.

@steveharoz
Last active October 26, 2023 21:06
Show Gist options
  • Save steveharoz/ddd6ec23628e41f06998732b93efc807 to your computer and use it in GitHub Desktop.
Save steveharoz/ddd6ec23628e41f06998732b93efc807 to your computer and use it in GitHub Desktop.
NOAA Storm and Hurricane data for this year
library(tidyverse)
library(rvest)
library(sf)
LOCATION = "https://ftp.nhc.noaa.gov/atcf/btk"
page = read_html(LOCATION)
files = page %>%
html_elements("a[href*='.dat']") %>%
html_text(trim = TRUE)
files = file.path(LOCATION, files)
# File format: https://www.nrlmry.navy.mil/atcf_web/docs/database/new/abdeck.txt
column_types = list(
basin = col_character(), # e.g. WP, IO, SH, CP, EP, AL, LS
cyclone_number = col_integer(), # 1-99
date = col_character(),
technum_min = col_integer(), # TECHNUM/MIN- objective technique sorting number, minutes for best track: 00 - 99
tech = col_character(),
forcast_period = col_integer(), # forecast period: -24 through 240 hours, 0 for best-track, negative taus used for CARQ and WRNG records.
lat = col_character(), # Latitude for the DTG: 0 - 900 tenths of degrees,
long = col_character(), # Longitude for the DTG: 0 - 1800 tenths of degrees,
time = col_character(),
wind_speed_max = col_integer(), # Maximum sustained wind speed in knots: 0 - 300 kts.
min_sea_level_pressure = col_character(), # Minimum sea level pressure, 850 - 1050 mb.
type = col_character(),
# DB - disturbance,
# TD - tropical depression,
# TS - tropical storm,
# TY - typhoon,
# ST - super typhoon,
# TC - tropical cyclone,
# HU - hurricane,
# SD - subtropical depression,
# SS - subtropical storm,
# EX - extratropical systems,
# PT - post tropical,
# IN - inland,
# DS - dissipating,
# LO - low,
# WV - tropical wave,
# ET - extrapolated,
# MD - monsoon depression,
# XX - unknown.
# check type???
rad = col_character(), # Wind intensity for the radii defined in this record: 34, 50 or 64 kt.
wind_code = col_character(), # radius code: AAA - full circle, NEQ, SEQ, SWQ, NWQ - quadrant
rad1 = col_integer(),
rad2 = col_integer(),
rad3 = col_integer(),
rad4 = col_integer(),
# If full circle, radius of specified wind intensity, or radius of
# first quadrant wind intensity as specified by WINDCODE. 0 - 999 n mi
pressure_outer = col_integer(), # pressure in millibars of the last closed isobar, 900 - 1050 mb.
radius_outer = col_integer(), # radius of the last closed isobar, 0 - 999 n mi.
radius_of_max_wind = col_integer(), # radius of max winds, 0 - 999 n mi.
gusts = col_integer(), # 0 - 999 kt.
eye_diameter = col_character(), # 0 - 120 n mi.
subregion = col_character(),
# subregion code: W,A,B,S,P,C,E,L,Q.
# A - Arabian Sea
# B - Bay of Bengal
# C - Central Pacific
# E - Eastern Pacific
# L - Atlantic
# P - South Pacific (135E - 120W)
# Q - South Atlantic
# S - South IO (20E - 135E)
# W - Western Pacific
max_seas = col_integer(), # 0 - 999 ft.
initials = col_character(), # Forecaster's initials used for tau 0 WRNG or OFCL
direction = col_integer(), # 0 - 359 degrees.
speed = col_character(), # 0 - 999 kts.
name = col_character(),
# literal storm name, number, NONAME or INVEST, or TCcyx where:
# cy = Annual cyclone number 01 - 99
# x = Subregion code: W,A,B,S,P,C,E,L,Q.
depth = col_character(), # system depth,
# D - deep,
# M - medium,
# S - shallow,
# X - unknown
seas = col_integer(), # Wave height for radii defined in SEAS1 - SEAS4, 0 - 99 ft.
seas_code = col_character(),
# AA - full circle
# NEQ, SEQ, SWQ, NWQ - quadrant
seas1 = col_integer(), # 1st quadrant seas radius as defined by SEASCODE, 0 - 999 n mi.
seas2 = col_integer(), # 2nd quadrant seas radius as defined by SEASCODE, 0 - 999 n mi.
seas3 = col_integer(), # 3rd quadrant seas radius as defined by SEASCODE, 0 - 999 n mi.
seas4 = col_character(), # 4th quadrant seas radius as defined by SEASCODE, 0 - 999 n mi.
user_defined1 = col_character(),
user_data1 = col_character(),
user_defined2 = col_character(),
user_data2 = col_character(),
user_defined3 = col_character(),
user_data3 = col_character(),
user_defined4 = col_character(),
user_data4 = col_character(),
user_defined5 = col_character(),
user_data5 = col_character(),
column39 = col_character(),
column40 = col_character()
)
column_names = names(column_types)
### read files
data = map(files, .progress = TRUE, function(filename) {
data_raw = read_lines(filename) %>% paste(collapse = "\n")
read_csv(data_raw, col_names = column_names, col_types = column_types, na = c("", "-99", "-999"))
})
data = bind_rows(data)
### clean data
data = data %>%
separate_wider_regex(lat, c(latitude="\\d*", lat_ns=".*")) %>%
separate_wider_regex(long, c(longitude="\\d*", long_ew=".*")) %>%
mutate(latitude = as.numeric(latitude) * ifelse(lat_ns == "N", .1, -.1)) %>%
mutate(longitude = as.numeric(longitude) * ifelse(long_ew == "E", .1, -.1)) %>%
select(-lat_ns, -long_ew)
### plot
# Load world map
world = st_as_sf(maps::map("world", plot = FALSE, fill = TRUE))
# check world map
ggplot(world) +
geom_sf(color = "white", fill = "gray20") +
theme_void()
ggplot(data %>% filter(basin == "AL")) +
geom_sf(data = world, color = "#222222", fill = "#555555") +
geom_path(aes(x = longitude, y = latitude, color = paste(basin, cyclone_number)), linewidth = 1) +
coord_sf(xlim = c(-110, -10), ylim = c(0, 55)) +
guides(color = "none") +
labs(x = NULL, y = NULL, title = "Atlantic Storm in 2023") +
theme_dark(13) +
theme(panel.background = element_rect("#222222"), panel.grid.major = element_line(alpha("black", 0.5)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment