In holstius/Rcaline: Line source dispersion modeling using CALINE

Introduction

Here we apply Rcaline to San Francisco data provided by the Bay Area Air Quality Management District (BAAQMD). This vignette illustrates the basic steps of constructing, running, and visualizing a model with Rcaline. We discuss how to load, visualize, and export model data with the use of several third-party packages.

Model construction

library(knitr)
options(width = 80, digits = 3, continue = " ")
opts_chunk$set(echo = FALSE, message = FALSE, cache = TRUE)
knit_hooks$set(small.mar = function(before, options, envir) {
    if (before) par(mar = c(.1, .1, .1, .1))  # smaller margin on top and right
})

Traffic data

Example traffic data has been supplied for you in the form of a SpatialLinesDataFrame called SF_highways.

library(Rcaline)

data(SanFrancisco, package = "Rcaline")
summary(geometry(SF_highways))
str(SF_highways@data)

It was originally imported from a shapefile. Should you have traffic data of your own in shapefile format, you can use readOGR() from the rgdal package to import it. The next section will describe how to set it up for use in your model.

See help(readOGR, package = "rgdal") and help(SpatialLinesDataFrame, package = "sp") for more on these data structures and their applications.

Meteorological data

Another example input comes in the form of an ASCII file, met_5801.isc, that contains $24 \times 365 = 8760$ records containing a year's worth of data on:

• wind speed (m/s)
• wind bearing (degrees)
• atmospheric stability (Pasquill A-F)
• mixing height (meters)

... at 1-hour scale.

met_file <- system.file("extdata", "BayArea", "Meteorology", "met_5801.isc", package = "Rcaline")
met_data <- ISCFile(met_file)

If you have meteorology data of your own in the same "ISC-ready" file format, you can import it with ISCFile(...) and use it instead. You can also run Rcaline with "manually" specified meteorology. Just modify the code above, substituting your own data.frame for met_data. (Make sure that your data.frame has appropriately named columns. See help(Meteorology) for more.)

SF_conditions <- Meteorology(met_data, use = "urban")
summary(SF_conditions)

You will be warned if there are any wind speeds less than 1.0 m/s. CALINE3 is premised on the assumption that advection dominates diffusion---in other words, that conditions are not calm. Rcaline will still compute results for non-calm conditions. Care should be taken to ensure that the proportion of calm conditions is very small.

Emission data

SF_highways has an attribute TRVol2009 that describes the Annual Average Daily Traffic (AADT) volume along polylines, each of which may be composed of one or more straight segments or links.

CALINE3 models dispersion at the level of individual links, so Rcaline splits each polyline into individual segments, assigning to each segment the attributes of the polyline to which it belonged.

SF_links <- FreeFlowLinks(
  SF_highways,
    vehiclesPerHour = TRVol2009 / 24, # convert from 24 h to 1 h scale
    emissionFactor = 1.0,             # in grams per vehicle-mile
    width = 30.0)                     # in meters

Model inputs can be specified in the form of expressions that depend on the data (as with vehiclesPerHour, above), or as fixed numeric constants. In this example---for illustrative purposes only---we specify a unform road width ($30$ meters) and a uniform emission factor ($1.0 \frac{g}{veh \cdot mi}$).

In California, pollutant-specific emission factors can be obtained from EMFAC.

Receptor locations

Receptors are the locations at which the dispersed pollutant concentration will be calculated. Since CALINE3 is a steady-state Gaussian plume model, rather than a numerical grid-based model, the locations of receptors can be freely specified.

You might want to compute predicted concentrations at a set of specific locations, like geocoded street addresses. Load them as follows:

## Not run:
# fn <- file.choose()
# layer_name <- sub("\\.shp$", "", fn, ignore.case = TRUE)
# features <- spTransform(readOGR(dirname(fn), layer_name), CRS(proj4string(SF_county)))
# stopifnot(inherits(features, "SpatialPoints"))
# receptors <- Receptors(features)

Here we show how to generate receptors on a regular Cartesian grid within 1 km of highway centerlines:

SF_receptor_grid <- ReceptorGrid(
  SF_links, 
  resolution = 250, 
  maxDistance = 1e3)

However, relying on a Cartesian grid can result in apparent "hot spots" that are just the result of certain grid points falling very close to links.

An alternative approach is to generate receptor locations at regular distances from the road network. This approach also allows receptors to be packed more densely close to the roadway, so that modeling resources can be concentrated on the more interesting locations (i.e., those closer to the roadways).

SF_receptor_rings <- ReceptorRings(
  SF_links, 
  distances = c(100, 250, 500, 1000))

plot(SF_county, col = "light gray", border = NA)
lines(SF_links)
points(SF_receptor_grid, pch = '+', cex = 0.5)

plot(SF_county, col = "light gray", border = NA)
lines(SF_links)
points(SF_receptor_rings, pch = '+', cex = 0.5)

If you were interested in estimating aggregate exposures at a population level, you could also construct an irregular grid by sampling locations from predefined regions, like ZIP codes or census tracts.

Terrain and pollutant properties

For detailed information on terrain and pollutant characteristics, consult the CALINE3 User's Guide. Here we specify some reasonable default values:

urban_terrain <- Terrain(surfaceRoughness = 80.0)

fine_PM <- Pollutant("PM2.5", molecularWeight = NA)
show(fine_PM)

Predicting concentrations

The CALINE3 algorithm is CPU-intensive. Although it is beyond the scope of this example, you can use the foreach package to do computations in parallel, using multiple cores or networked hosts. (A future vignette will illustrate this technique.)

Running the model

We use the predict method to actually run the model. Since the model will actually be run once for every meteorological condition we supply ($n = r nrow(SF_conditions)$), it can be helpful to make a first pass using only a small sample of the meteorology.

sample_rows <- function(x, p) {
  x[sample(1:nrow(x), p * nrow(x)),]
}

SF_model <- Caline3Model(
  SF_links, 
  sample_rows(SF_conditions, p = 0.01), 
  SF_receptor_rings, 
  urban_terrain, 
  fine_PM)

pred_matrix <- predict(SF_model, units = "mg/m^3")

The result of running the model is an $M \times N$ HourlyConcentrations matrix, where $M$ is the number of meteorological conditions and $N$ is the number of receptors. Each cell describes the predicted concentration at that receptor during those conditions.

Summarizing results

The aggregate.HourlyConcentrations helper function will compute several statistics by default, including the mean and maximum.

pred_stats <- aggregate(pred_matrix)
show(colnames(pred_stats))

Casting the resulting AggregatedConcentrations object to a SpatialPointsDataFrame then re-binds the statistics to the receptor locations.

pred_spatial <- as(pred_stats, "SpatialPointsDataFrame")
show(pred_spatial[1:3, c("distance", "min", "mean", "median", "max")])

Afterwards, we can select a statistic of interest and explore the distribution. Here we focus on exploring results graphically, although they could also be tabulated or subjected to statistical tests.

Distance-to-roadway

Within the results is a variable, distance, that describes the distance-to-roadway for each receptor. (Recall that we specified these distances when constructing the receptor grid.) Using this, we can group the receptors into specific classes and explore the distribution of predicted concentrations in each.

suppressPackageStartupMessages(library(ggplot2))
ggplot(pred_spatial@data) + geom_histogram(aes(mean)) + facet_wrap(~ distance)

Mapping

Rcaline provides a helper function ggplot.AggregatedConcentrations that automates the re-binding of predicted concentrations with receptor locations, as well as the plotting of model inputs such as receptors and links. Here is an example of a "bubble" plot generated with ggplot.AggregatedConcentrations:

label_kilo <- function (x) sprintf("%0.0f", x / 1000.0)
ggplot(pred_stats) + 
  geom_point(aes(x, y, size = mean, color = mean, order = mean, alpha = mean)) +
  scale_color_gradient(expression(paste(bold(hat(Z)(s)), " ", mg/m^3)), low = "#BEBEBE00", high = "brown") + 
  scale_size(guide = "none") + scale_alpha(guide = "none") + 
  scale_x_continuous("UTM10 Easting (km)", labels = label_kilo) +
  scale_y_continuous("UTM10 Northing (km)", labels = label_kilo)

holstius/Rcaline documentation built on May 17, 2019, 4:39 p.m.