/
getMeta.R
187 lines (149 loc) · 6.84 KB
/
getMeta.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
##' Get information on meteorological sites
##'
##' This function is primarily used to find a site code that can be
##' used to access data using \code{\link{importNOAA}}. Sites searches
##' of approximately 30,000 sites can be carried out based on the site
##' name and based on the nearest locations based on user-supplied
##' latitude and logitude.
##' @title Find a ISD site code and other meta data
##' @param site A site name search string e.g. \code{site =
##' "heathrow"}. The search strings and be partial and can be upper
##' or lower case e.g. \code{site = "HEATHR"}.
##' @param lat A latitude in decimal degrees to search. Takes the
##' values -90 to 90.
##' @param lon A longitude in decimal degrees to search. Takes values
##' -180 to 180. Negative numbers are west of the Greenwich
##' meridian.
##' @param country The country code. This is a two letter code. For a
##' full listing see
##' \url{ftp://ftp.ncdc.noaa.gov/pub/data/noaa/country-list.txt}.
##' @param state The state code. This is a two letter code.
##' @param n The number of nearest sites to search based on
##' \code{latitude} and \code{longitude}.
##' @param end.year To help filter sites based on how recent the
##' available data are. \code{end.year} can be "current", "any" or a
##' numeric year such as 2016, or a range of years e.g. 1990:2016
##' (which would select any site that had an end date in that range.
##' \strong{By default only sites that have some data for the
##' current year are returned}.
##' @param plot If \code{TRUE} will plot sites on an interactive
##' leaflet map.
##' @param fresh Should the meta data be read from the NOAA server or
##' the \code{worldmet} package?. If \code{FALSE} it is read from
##' the package version, which is fast but might be out of date. If
##' \code{TRUE} the data are read from the NOAA server.
##' @param returnMap Should the leaflet map be returned instead of the
##' meta data? Default is \code{FALSE}.
##' @return A data frame is returned with all available meta data,
##' mostly importantly including a \code{code} that can be supplied
##' to \code{\link{importNOAA}}. If latitude and longitude searches
##' are made an approximate distance, \code{dist} in km is also
##' returned.
##' @export
##' @author David Carslaw
##' @examples
##'
##' \dontrun{
##' ## search for sites with name beijing
##' getMeta(site = "beijing")
##' }
##'
##' \dontrun{
##' ## search for near a specified lat/lon - near Beijing airport
##' ## returns 'n' nearest by default
##' getMeta(lat = 40, lon = 116.9)
##' }
getMeta <- function(site = "heathrow", lat = NA, lon = NA,
country = NA, state = NA, n = 10, end.year = "current",
plot = TRUE, fresh = TRUE, returnMap = FALSE) {
## read the meta data
# check year
if (!any(end.year %in% c("current", "all"))) {
if (!is.numeric(end.year)) {
stop("end.year should be one of 'current', 'all' or a numeric 4-digit year such as 2016.")
}
}
# we base the current year as the max available in the meta data
if ("current" %in% end.year) end.year <- max(as.numeric(format(meta$END, "%Y")), na.rm = TRUE)
if ("all" %in% end.year) end.year <- 1900:2100
## download the file, else use the package version
if (fresh) meta <- getMetaLive()
## search based on name of site
if (!missing(site)) {
## search for station
meta <- meta[grep(site, meta$STATION, ignore.case = TRUE), ]
}
## search based on country codes
if (!missing(country) && !is.na(country)) {
## search for country
id <- which(meta$CTRY %in% toupper(country))
meta <- meta[id, ]
}
## search based on state codes
if (!missing(state)) {
## search for state
id <- which(meta$ST %in% toupper(state))
meta <- meta[id, ]
}
# make sure no missing lat / lon
id <- which(is.na(meta$LON))
if (length(id) > 0) meta <- meta[-id, ]
id <- which(is.na(meta$LAT))
if (length(id) > 0) meta <- meta[-id, ]
# filter by end year
id <- which(format(meta$END, "%Y") %in% end.year)
meta <- meta[id, ]
## approximate distance to site
if (!missing(lat) && !missing(lon)) {
r <- 6371 # radius of the Earth
## Coordinates need to be in radians
meta$longR <- meta$LON * pi / 180
meta$latR <- meta$LAT * pi / 180
LON <- lon * pi / 180
LAT <- lat * pi / 180
meta$dist <- acos(sin(LAT) * sin(meta$latR) + cos(LAT) *
cos(meta$latR) * cos(meta$longR - LON)) * r
## sort and retrun top n nearest
meta <- head(openair:::sortDataFrame(meta, key = "dist"), n)
}
dat <- dplyr::rename(meta, latitude = LAT, longitude = LON)
if (plot) {
if (!"dist" %in% names(dat)) dat$dist <- NA
content <- paste(paste(dat$STATION,
paste("Code:", dat$code),
paste("Start:", dat$BEGIN),
paste("End:", dat$END),
paste("Distance (km)", round(dat$dist, 1)),
sep = "<br/>"))
m <- leaflet(dat) %>% addTiles() %>%
addMarkers(~longitude, ~latitude, popup = content)
if (!is.na(lat) && !is.na(lon))
m <- m %>% addCircles(lng = lon, lat = lat, weight = 20, radius = 200,
stroke = TRUE, color = "red",
popup = paste("Search location",
paste("Lat =", lat),
paste("Lon =", lon), sep = "<br/>"))
print(m)
}
if (returnMap) return(m) else return(dat)
}
getMetaLive <- function(...) {
## downloads the whole thing fresh
url <- "https://www1.ncdc.noaa.gov/pub/data/noaa/isd-history.csv"
meta <- suppressMessages(read.csv(url, skip = 21))
## names in the meta file
names(meta) <- c("USAF", "WBAN","STATION", "CTRY", "ST", "CALL", "LAT",
"LON", "ELEV(M)", "BEGIN", "END")
## full character string of site id
meta$USAF <- formatC(meta$USAF, width = 6, format = "d", flag = "0")
## start/end date of measurements
meta$BEGIN <- as.Date(as.character(meta$BEGIN), format = "%Y%m%d")
meta$END <- as.Date(as.character(meta$END), format = "%Y%m%d")
## code used to query data
meta$code <- paste(meta$USAF, meta$WBAN, sep = "-")
return(meta)
}
# how to update meta data
# meta <- getMeta(fresh = TRUE)
# devtools::use_data(meta, overwrite = TRUE, internal = TRUE)
# devtools::use_data(meta, overwrite = TRUE)