Wind direction and speed • aire.zmvm

We can use the idw360 function to extrapolate wind direction and speed within a grid. First, we load the packages:

## Auto-install required R packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(aire.zmvm, dplyr, ggplot2, lubridate, stringr, sp, mapproj,
               gstat, ggmap)

and then download the wind speed and direction data:

# download wind speed and direction and filter March 28, 2018 7:00 PM
wsp <- get_station_month_data(criterion = "HORARIOS",
                              pollutant = "WSP",
                              2018, 3)%>%
  filter(date == "2018-03-28" & hour == 19) %>%
  na.omit()
wdr <- get_station_month_data(criterion = "HORARIOS", 
                        pollutant = "WDR", 
                        year = 2018, 3) %>%
  filter(date == "2018-03-28" & hour == 19)  %>%
  na.omit()

This is what the data looks like:

knitr::kable(head(wsp))

	date	hour	station_code	pollutant	unit	value
2	2018-03-28	19	AJM	WSP	m/s	10.2
5	2018-03-28	19	BJU	WSP	m/s	4.6
9	2018-03-28	19	MON	WSP	m/s	6.5
10	2018-03-28	19	CHO	WSP	m/s	3.5
12	2018-03-28	19	CUA	WSP	m/s	4.3
14	2018-03-28	19	FAC	WSP	m/s	3.8

knitr::kable(head(wdr))

	date	hour	station_code	pollutant	unit	value
2	2018-03-28	19	AJM	WDR	°	198
5	2018-03-28	19	BJU	WDR	°	169
9	2018-03-28	19	MON	WDR	°	165
10	2018-03-28	19	CHO	WDR	°	150
12	2018-03-28	19	CUA	WDR	°	246
14	2018-03-28	19	FAC	WDR	°	205

We can create a point grid to extrapolate the wind values:

## Location of sensors. First column x/longitude, second y/latitude
locations <- stations[, c("station_code", "lon", "lat")]
locations <- merge(wdr, locations)
coordinates(locations) <- ~lon+lat
proj4string(locations) <- CRS("+proj=longlat +ellps=WGS84 +no_defs +towgs84=0,0,0")

## The grid for which to extrapolate values
pixels = 30
grid <- expand.grid(lon = seq((min(coordinates(locations)[, 1]) - .1),
                              (max(coordinates(locations)[, 1]) + .1),
                              length.out = pixels),
                    lat = seq((min(coordinates(locations)[, 2]) - .1),
                              (max(coordinates(locations)[, 2]) + .1),
                              length.out = pixels))
grid<- SpatialPoints(grid)
proj4string(grid) <- CRS("+proj=longlat +ellps=WGS84 +no_defs +towgs84=0,0,0")

plot(grid, main = "Grid covering Mexico City and Locations of\nWind Measuring Stations (red)",
     col = "#666666", lwd = .2)
# Plot the locations of the stations
geog <- stations[stations$station_code %in% unique(wsp$station_code), c("lat", "lon")]
coordinates(geog) <- ~lon+lat
points(geog, pch = 22, col = "#333333", bg = "tomato1", cex = 1.2)

Extrapolate using inverse distance weighting:

# extrapolate wind direction
res <- idw360(wdr$value, locations, grid)

# extrapolate wind speed
wsp <- left_join(wsp, stations, by = "station_code")
wsp <- wsp[!is.na(wsp$value),]
wsp <- wsp[,c("lat", "lon", "value")]
coordinates(wsp) <- ~lon+lat
proj4string(wsp) <- CRS("+proj=longlat +ellps=WGS84 +no_defs +towgs84=0,0,0")
wsp <- idw(value ~ 1, wsp, grid, idp = 2, debug.level = 0)
wsp = as.data.frame(wsp)
names(wsp) <- c("lon", "lat", "var1.wsp", "var1.var")

# Combine wind direction and speed with the grid
df <- cbind(res, as.data.frame(grid), wsp = wsp$var1.wsp)

and finally we plot the data (Note that you may have to install a beta version of ggmap by running: devtools::install_github("dkahle/ggmap", ref = "tidyup"))

df$group <- 1
register_stadiamaps(Sys.getenv("STADIA_KEY"), write = FALSE)
qmplot(lon, lat, data = df, geom = "blank", 
       maptype = "alidade_smooth", source = "stadia", zoom = 10, mapcolor = "bw") +
  geom_point(size = .1) +
  # The wind direction compass starts where the 90 degree mark is located
  geom_spoke(aes(angle = ((270 -  pred) %% 360) * pi / 180,
                 radius = wsp * .003),
             arrow = arrow(length = unit(0.09, "cm"))) +
  ggtitle("Wind direction and speed: March 28, 2018 at 7:00 PM")