Understanding speed-dependent inventories with emis

This article explains, step by step, how VEIN builds an hourly, per-street inventory with speed-dependent emission factors, and how emis_speed() makes it fast enough for country-scale networks. Every chunk runs; you can copy and paste them.

1. The idea

A local emission factor (for example the CETESB factors used in Brazil) is a single number per vehicle age, measured over a driving cycle. Real emissions depend on the average speed of the link. The EMEP/EEA guidebook provides speed curves for the same technologies. VEIN combines both: the EMEP/EEA curve gives the shape as a function of speed, rescaled so that at the speed of the local driving cycle it reproduces the local emission factor.

Doing this for every street and every hour used to be slow because the curve was interpreted in R once per street, hour and age. emis_speed() compiles the curves and evaluates them in a C/OpenMP kernel.

2. Local emission factors by age

ef_cetesb() returns one emission factor per age. full = TRUE also gives the Euro equivalence standard of each age, which we need to pick the right curve.

library(vein)

CO <- ef_cetesb(p = "CO", veh = "PC_G", year = 2018, agemax = 40, full = TRUE)
head(CO[, c("Age", "Year", "Proconve_PC", "EqEuro_PC", "CO")], 8)
#>      Age Year Proconve_PC EqEuro_PC           CO
#> 1051   3 2018          L6         V 0.173 [g/km]
#> 1052   4 2017          L6         V 0.141 [g/km]
#> 1053   5 2016          L6         V 0.114 [g/km]
#> 1054   6 2015          L6         V 0.155 [g/km]
#> 1055   7 2014          L5        IV 0.211 [g/km]
#> 1056   8 2013          L5        IV 0.241 [g/km]
#> 1057   9 2012          L5        IV 0.274 [g/km]
#> 1058  10 2011          L5       III 0.274 [g/km]

More than one call to ef_cetesb() with the same year/scale is cheap: the large, pollutant-independent CETESB table is cached internally.

3. Scaling a curve to the driving cycle (SDC)

ef_ldv_scaled() (light vehicles) and ef_hdv_scaled() (heavy vehicles and buses) return one function per age. Each function is EF_a(V) = k_a * curve_a(V) where k_a = EF_local[a] / curve_a(SDC), so that at the driving-cycle speed SDC we recover the local factor exactly.

lef <- ef_ldv_scaled(
  dfcol = CO$CO, SDC = 34.12,
  v = "PC", t = "4S", cc = "<=1400", f = "G", eu = CO$EqEuro_PC, p = "CO"
)

length(lef) # one function per age
#> [1] 40
# at the driving-cycle speed the scaled function equals the local factor
c(local = as.numeric(CO$CO[1]), scaled = as.numeric(lef[[1]](34.12)))
#>  local scaled 
#>  0.173  0.173

Once scaled, the factor depends on speed. Ages with a newer Euro standard have a different curve, so the shape also changes with age:

V <- seq(0, 130, by = 10)
round(sapply(c(1, 10, 20, 40), function(a) as.numeric(lef[[a]](V))), 3)
#>        [,1]  [,2]  [,3]   [,4]
#>  [1,] 0.190 0.224 1.906 71.501
#>  [2,] 0.190 0.224 1.906 71.501
#>  [3,] 0.162 0.195 1.072 46.202
#>  [4,] 0.167 0.196 0.802 35.787
#>  [5,] 0.182 0.206 0.678 29.855
#>  [6,] 0.195 0.222 0.616 25.939
#>  [7,] 0.202 0.244 0.591 23.125
#>  [8,] 0.207 0.274 0.591 20.984
#>  [9,] 0.222 0.314 0.615 19.291
#> [10,] 0.262 0.370 0.670 17.912
#> [11,] 0.346 0.454 0.769 16.846
#> [12,] 0.500 0.589 0.954 18.062
#> [13,] 0.746 0.844 1.351 19.277
#> [14,] 1.113 1.500 2.645 20.493

4. What `emis_speed()` consumes: compiled programs

Evaluating R closures for millions of street-hours is slow. Instead we compile each curve to a tiny stack machine and group the ages that share the same equation. programs_only = TRUE returns just that bundle.

prog <- ef_ldv_scaled(
  dfcol = CO$CO, SDC = 34.12,
  v = "PC", t = "4S", cc = "<=1400", f = "G", eu = CO$EqEuro_PC, p = "CO",
  programs_only = TRUE
)
class(prog)
#> [1] "speed_programs"
c(ages = prog$n, groups = prog$G) # 40 ages collapse into a few equations
#>   ages groups 
#>     40      6
str(prog, max.level = 1)
#> List of 13
#>  $ n     : int 40
#>  $ G     : int 6
#>  $ gid   : int [1:40] 0 0 0 0 1 1 1 2 2 2 ...
#>  $ kk_age: num [1:40] 0.71 0.578 0.468 0.636 1.153 ...
#>  $ code  : int [1:152] 10 17 1 0 24 22 11 17 1 1 ...
#>  $ clen  : int [1:6] 33 26 26 26 26 15
#>  $ consts: num [1:17] 5 4 3 2 2 1 2 2 1 2 ...
#>  $ soff  : int [1:6] 0 4 7 10 13 16
#>  $ cofs  : num [1:36] -1.35e-10 7.86e-08 -1.22e-05 7.75e-04 -1.97e-02 3.98e-01 1.36e-01 -1.41e-02 -8.91e-04 4.99e-05 ...
#>  $ x     : num [1:6] 0 0 0 0 0 0
#>  $ minv  : num [1:6] 10 10 10 10 10 10
#>  $ maxv  : num [1:6] 130 130 130 130 130 130
#>  $ progs :List of 40
#>  - attr(*, "class")= chr "speed_programs"

The important fields:

code, consts: the RPN bytecode of each distinct equation.
cofs, minv, maxv: coefficients a..f and the speed clamp per group.
gid: for each age, which group (equation) it uses.
kk_age: the per-age scaling constant k_a.

5. Speeds per street and hour

Emissions need a speed for each link and each hour of the week. temp_fact() expands the traffic flow with a temporal profile and netspeed() turns flow into speed with a BPR function.

data(net)
data(pc_profile)

total <- net$ldv + net$hdv
flow  <- temp_fact(total, pc_profile)                                  # 168 h
speed <- netspeed(flow, net$ps, net$ffs, net$capacity, net$lkm, alpha = 1)
dim(speed) # streets x (24 h * 7 days)
#> [1] 1505  168

6. Running the inventory

veh is a matrix of vehicle flow (veh/h) by street (rows) and age (columns). emis_speed() returns the emissions by street x hour and, with by_age = TRUE, the totals by age x hour.

A <- 40
veh <- matrix(as.numeric(net$ldv) / A, nrow = nrow(net), ncol = A)

E <- emis_speed(
  veh = veh, lkm = net$lkm, ef = lef,
  speed = speed, profile = pc_profile,
  by_age = TRUE, nt = 2
)

dim(E$streets) # streets x hours
#> [1] 1505  168
dim(E$veh)     # ages x hours
#> [1]  40 168
sum(E$streets) # grams per hour over the network
#> [1] 795654179

The units of E$streets are g/h because veh is veh/h, lkm is km and the emission factor is g/km.

7. It matches the reference implementation

emis_speed() is a drop-in replacement for the speed branch of emis(); the results agree to machine precision.

Eref <- emis(
  veh = veh, lkm = net$lkm, ef = lef,
  speed = speed, profile = pc_profile,
  simplify = TRUE, agemax = A
)

max(abs(apply(Eref, c(1, 3), sum) - E$streets))
#> [1] 1.047738e-09

8. Under the hood

Compile. Each Y equation in sysdata$ldv / sysdata$hdv_* is parsed once into RPN bytecode (compile_ef_expr()); the result is cached.
Group. Ages with identical equation and coefficients (usually several years within one Euro standard) are merged, so a curve is evaluated once per street-hour instead of once per age. For this example the 40 ages need only 6 equations.
Evaluate. src/e_speed.c runs the bytecode in a tight OpenMP loop and accumulates the totals directly, never materialising the full street x age x hour array.

A quick timing comparison:

S <- 1000
vh <- matrix(50, S, A)
sp <- Speed(matrix(40, S, H <- 168))
pf <- matrix(1, 24, 7)
lkm1 <- units::set_units(rep(1, S), "km")

t_ref <- system.time(
  E1 <- emis(veh = vh, lkm = lkm1, ef = lef, speed = sp, profile = pf,
             simplify = TRUE, agemax = A)
)["elapsed"]
t_new <- system.time(
  E2 <- emis_speed(veh = vh, lkm = lkm1, ef = lef, speed = sp, profile = pf,
                   agemax = A, nt = 2)
)["elapsed"]
round(c(emis = t_ref, emis_speed = t_new), 3)
#>       emis.elapsed emis_speed.elapsed 
#>              1.100              0.058

9. Country-scale workflow: the `brazil_bu_speed` project

The complete workflow is packaged as a project you can download and run:

get_project(directory = tempdir(), case = "brazil_bu_speed")

Key files:

config/ef_speed_mapping.R maps every CETESB category (e.g. TRUCKS_M_D) to an EMEP/EEA category (v, t, cc/g, f) and the Euro column of ef_cetesb(). It also defines ef_cetesb_speed(), the scaled-EF builder.
scripts/speed.R computes the 168 hourly speeds with netspeed().
scripts/exhaust_speed.R loops over categories and pollutants, runs emis_speed() and falls back to a constant emission factor when no EMEP curve applies.
scripts/post_speed.R grids the hourly street emissions and writes totals.

10. Tips and pitfalls

Not every curve is usable. For some pollutant/Euro combinations the EMEP curve is zero at SDC (for example pre-Euro PM), so the scaling constant is not finite. emis_speed() stops with a clear message; the project detects this and uses the constant CETESB factor instead.
Threads. nt controls OpenMP threads. The default is half of check_nt(), which works well when R and the kernel must share the machine.
Memory. For huge networks the per-street result is only streets x hours (e.g. 200 000 x 168 is about 270 MB), because ages are summed in C. Ask for by_age = TRUE to also get the small ages x hours totals.
Speed layout. The columns of speed must match the flattened profile (hours of day 1, then day 2, …). temp_fact() produces this order.

Understanding speed-dependent inventories with emis_speed