How to generate inputs for the MUNICH model from the VEIN model

Question

MUNICH is the Model of Urban Network of Intersecting Canyons and Highways (MUNICH) (Kim et al., 2018) and VEIN is the Vehicular Emissions INventory model (Ibarra-Espinosa et al, 2018), an R package.

VEIN estimates vehicular emissions at street level, generating useful inputs for the MUNICH model, however, some pre-processing is necessary. Therefore, the model eixport was created to help in this processing.

As a maintainer I receive many emails asking how to generate inputs for the MUNICH model

• Kim, Y., Wu, Y., Seigneur, C., and Roustan, Y.: Multi-scale modeling of urban air pollution: development and application of a Street-in-Grid model (v1.0) by coupling MUNICH (v1.0) and Polair3D (v1.8.1), Geosci. Model Dev., 11, 611-629, https://doi.org/10.5194/gmd-11-611-2018, 2018.
• Ibarra-Espinosa, S., Ynoue, R., O'Sullivan, S., Pebesma, E., Andrade, M. D. F., and Osses, M.: VEIN v0.2.2: an R package for bottom–up vehicular emissions inventories, Geosci. Model Dev., 11, 2209–2229, https://doi.org/10.5194/gmd-11-2209-2018, 2018.
• Ibarra-Espinosa et al., (2018b). eixport: An R package to export emissions to atmospheric models. Journal of Open Source Software, 3(24), 607, https://doi.org/10.21105/joss.00607

Answer 1

First, download a project with vein and estimate vehicular emissions. Some projects are here: https://atmoschem.github.io/vein/reference/get_project.html

library(vein)
get_project(directory = "awesomecity")

Then open the file main.Rproj with Rstudio and source the file main.R

Then, run the following code

library(data.table)
library(eixport)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(writexl)
library(vein)
library(units)
#> udunits system database from /usr/share/xml/udunits
# epsg 31983, projected UTM
co <- readRDS("post/streets/CO.rds")
st_crs(co)
#> Coordinate Reference System:
#>   User input: EPSG:31983 
#>   wkt:
#> PROJCS["SIRGAS 2000 / UTM zone 23S",
#>     GEOGCS["SIRGAS 2000",
#>         DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
#>             SPHEROID["GRS 1980",6378137,298.257222101,
#>                 AUTHORITY["EPSG","7019"]],
#>             TOWGS84[0,0,0,0,0,0,0],
#>             AUTHORITY["EPSG","6674"]],
#>         PRIMEM["Greenwich",0,
#>             AUTHORITY["EPSG","8901"]],
#>         UNIT["degree",0.0174532925199433,
#>             AUTHORITY["EPSG","9122"]],
#>         AUTHORITY["EPSG","4674"]],
#>     PROJECTION["Transverse_Mercator"],
#>     PARAMETER["latitude_of_origin",0],
#>     PARAMETER["central_meridian",-45],
#>     PARAMETER["scale_factor",0.9996],
#>     PARAMETER["false_easting",500000],
#>     PARAMETER["false_northing",10000000],
#>     UNIT["metre",1,
#>         AUTHORITY["EPSG","9001"]],
#>     AXIS["Easting",EAST],
#>     AXIS["Northing",NORTH],
#>     AUTHORITY["EPSG","31983"]]
head(co, 1) # units g/h
#> Simple feature collection with 1 feature and 25 fields
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 328001.4 ymin: 7391129 xmax: 328036.5 ymax: 7391204
#> projected CRS:  SIRGAS 2000 / UTM zone 23S
#>        id            V1             V2             V3             V4
#> 1 1 [1/h] 275.088 [g/h] 140.0241 [g/h] 84.63286 [g/h] 86.81989 [g/h]
#>               V5             V6             V7         V8             V9
#> 1 195.8425 [g/h] 680.7458 [g/h] 1423.767 [g/h] 1719 [g/h] 1723.265 [g/h]
#>              V10            V11            V12            V13          V14
#> 1 1748.996 [g/h] 1600.077 [g/h] 1448.586 [g/h] 1378.825 [g/h] 1551.5 [g/h]
#>              V15            V16            V17            V18            V19
#> 1 1556.648 [g/h] 1538.186 [g/h] 1698.133 [g/h] 1655.548 [g/h] 1601.138 [g/h]
#>              V20            V21            V22            V23            V24
#> 1 1398.703 [g/h] 941.4398 [g/h] 709.3241 [g/h] 581.2601 [g/h] 407.1127 [g/h]
#>                         geometry
#> 1 LINESTRING (328001.4 739120...

# example of  polygon as buffer at center
st_bbox(co) %>%
  st_as_sfc() %>%
  st_centroid() %>%
  st_buffer(dist = 2000) -> polygon

# if you have another polygon, read it
# polygon <- st_read("/path/to/polygon.gpkg")

cob <- st_intersection(co, polygon)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

# Remove ID
co$id <- NULL

# split on vertex conserving mass
cos <- st_explode(cob)
#> Sum: 16542955.31
dim(cos)
#> [1] 372  27
# it adds length LKM and LKM2
# Should we remove streets shorter than 10 mts?
cos <- cos[as.numeric(cos$LKM2) > 10, ]

# remove LKM and LKM2
cos$LKM <- cos$LKM2 <- NULL

# We need to add units back, g/h
cos <- vein::Emissions(cos, time = "1/h")
head(cos, 1)
#> Simple feature collection with 1 feature and 24 fields
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 331916.3 ymin: 7385767 xmax: 331922.9 ymax: 7385775
#> projected CRS:  SIRGAS 2000 / UTM zone 23S
#>               V1             V2             V3             V4             V5
#> 1 59.38562 [g/h] 30.22821 [g/h] 18.27042 [g/h] 18.74256 [g/h] 42.27821 [g/h]
#>               V6             V7             V8             V9            V10
#> 1 146.9585 [g/h] 307.3609 [g/h] 371.0953 [g/h] 372.0159 [g/h] 377.5708 [g/h]
#>              V11            V12            V13            V14            V15
#> 1 345.4224 [g/h] 312.7188 [g/h] 297.6588 [g/h] 334.9356 [g/h] 336.0469 [g/h]
#>              V16            V17            V18            V19          V20
#> 1 332.0614 [g/h] 366.5907 [g/h] 357.3973 [g/h] 345.6514 [g/h] 301.95 [g/h]
#>              V21            V22            V23            V24
#> 1 203.2367 [g/h] 153.1279 [g/h] 125.4816 [g/h] 87.88691 [g/h]
#>                         geometry
#> 1 LINESTRING (331916.3 738577...
# now we transform for the REQUIRED UNITS FOR MUNICH
# ug/km/h

# we have 24 hours in this case
hours <- paste0("V", 1:24)

for (i in seq_along(hours)) {
  cos[[hours[i]]] <- set_units(cos[[hours[i]]], ug / h)
  cos[[hours[i]]] <- cos[[hours[i]]] / set_units(st_length(cos), km)
}
plot(cos["V1"], axes = T)

etm <- to_munich(sdf = cos)
names(etm)
#> [1] "Emissions" "Street"
head(etm$Emissions, 1)
#>   i idbrin typo        xa        ya        xb        yb                   V1
#> 1 1      1    0 -46.64773 -23.62993 -46.64766 -23.63001 5462199660 [ug/h/km]
#>                     V2                   V3                   V4
#> 1 2780345076 [ug/h/km] 1680485867 [ug/h/km] 1723911918 [ug/h/km]
#>                     V5                    V6                    V7
#> 1 3888685780 [ug/h/km] 13517016236 [ug/h/km] 28270586844 [ug/h/km]
#>                      V8                    V9                   V10
#> 1 34132780053 [ug/h/km] 34217462106 [ug/h/km] 34728386913 [ug/h/km]
#>                     V11                   V12                   V13
#> 1 31771425772 [ug/h/km] 28763397272 [ug/h/km] 27378203813 [ug/h/km]
#>                     V14                   V15                   V16
#> 1 30806871891 [ug/h/km] 30909088139 [ug/h/km] 30542507720 [ug/h/km]
#>                     V17                   V18                   V19
#> 1 33718452534 [ug/h/km] 32872865695 [ug/h/km] 31792493506 [ug/h/km]
#>                     V20                   V21                   V22
#> 1 27772900557 [ug/h/km] 18693404826 [ug/h/km] 14084471885 [ug/h/km]
#>                     V23                  V24
#> 1 11541609851 [ug/h/km] 8083705065 [ug/h/km]
head(etm$Street, 1)
#>   i       length width height
#> 1 1 10.87211 [m]     0     30

# to export to CSV which supports bit64
data.table::fwrite(
  etm$Emissions,
  paste0(basename(getwd()), "_Emissions.csv")
)
data.table::fwrite(
  etm$Street,
  paste0(basename(getwd()), "_Streets.csv")
)

# to export to Excel
writexl::write_xlsx(
  x = etm$Emissions,
  path = paste0(basename(getwd()), "_Emissions.xlsx")
)
writexl::write_xlsx(
  x = etm$Street,
  path = paste0(basename(getwd()), "_Streets.xlsx")
)

^{Created on 2020-10-14 by the reprex package (v0.3.0)}