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

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)}