default case can be any of "brasil", "brazil", "brazil_bu", "brasil_bu", they are -the same
+the same +Projects for Ecuador are in development. +In any case, if you find any error, please, send a pull request in github or gitlab.if (FALSE) { diff --git a/docs/reference/index.html b/docs/reference/index.html index 68fad3fe..9c8168ee 100755 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -309,12 +309,6 @@ef_emfac() - -
- EMFAC2017 emission factors for Statewide California, Calendar Year 2020
- +diff --git a/docs/reference/vein-deprecated.html b/docs/reference/vein-deprecated.html index e9250dc6..bcc88cb8 100644 --- a/docs/reference/vein-deprecated.html +++ b/docs/reference/vein-deprecated.html @@ -52,6 +52,9 @@ units g. This class represents the daily emissions presented by Mellios G and Ntziachristos (2016) Gasoline evaporation, Tier 2. Eventually it will be incorporated the techniques of Tier 3. +ef_emfac returns emission factors reflecting California, US, +conditions. If the user enter speeds more emission factors are returned. +See details. emis_wrf returns a dataframes with columns lat, long, id, pollutants, local time and GMT time. This dataframe has the proper format to be used with WRF assimilation system: "ASimilation System 4 WRF (AS4WRF Vera-Vala et al (2016)) @@ -153,7 +156,7 @@ @@ -163,6 +166,9 @@DEPRECATED Construction function for class "EmissionsList"
- Source:R/EmissionsList.R,R/Evaporative.R,R/emis_wrf.R, and 3 more + Source:R/EmissionsList.R,R/Evaporative.R,R/ef_emfac.R, and 4 morevein-deprecated.RdDEPRECATED Construction function for class "EmissionsList"
units g. This class represents the daily emissions presented by Mellios G and Ntziachristos (2016) Gasoline evaporation, Tier 2. Eventually it will be incorporated the techniques of Tier 3. +
ef_emfacreturns emission factors reflecting California, US, +conditions. If the user enter speeds more emission factors are returned. +See details.@@ -176,6 +182,8 @@
emis_wrfreturns a dataframes with columns lat, long, id, pollutants, local time and GMT time. This dataframe has the proper format to be used with WRF assimilation system: "ASimilation System 4 WRF (AS4WRF Vera-Vala et al (2016))DEPRECATED Construction function for class "EmissionsList"
Evaporative(x, ...) +ef_emfac(x, ...) + emis_wrf(sdf, nr = 1, dmyhm, tz, crs = 4326, islist) hot_soak(x, carb, p, eshotc, eswarmc, eshotfi) @@ -226,12 +234,37 @@Arg
object Object with class "Evaporative"
+ +veh ++ Character; "one of the 40 vehicle categories shown below.
+ +fuel ++ Character; "Diesel", "Gasoline", "Electricity" or "Natural Gas"
+ +mph ++ Numeric; Speed in miles per hour (optional).
+ +pol ++ Character;
+ +season ++ Character: "winter" or "summer".
+ full ++ Logical: To return the whole data.table or not.
Value
Objects of class "EmissionsList" and numeric elements as "units"
Objects of class "Evaporative" or "units"
+data.table with emission factors.
@@ -243,7 +276,8 @@
EvaporativeFor
running_losses, useemis_evap.References
-Mellios G and Ntziachristos 2016. Gasoline evaporation. In: +
https://arb.ca.gov/emfac/emissions-inventory
+Mellios G and Ntziachristos 2016. Gasoline evaporation. In: EEA, EMEP. EEA air pollutant emission inventory guidebook-2009. European Environment Agency, Copenhagen, 2009
Mellios G and Ntziachristos 2016. Gasoline evaporation. In: @@ -253,6 +287,7 @@
See a
vein-deprecated+
vein-deprecated
vein-deprecated
vein-deprecatedExamples
@@ -266,6 +301,10 @@Examp } if (FALSE) { # do not run +# DEPRECATED +} +if (FALSE) { +# do not run # Use library(eixport) # eixport::to_as4wrf() } diff --git a/docs/sitemap.xml b/docs/sitemap.xml index 059117af..9749535e 100644 --- a/docs/sitemap.xml +++ b/docs/sitemap.xml @@ -75,9 +75,6 @@
- http://atmoschem.github.io/vein//reference/ef_eea.html - http://atmoschem.github.io/vein//reference/ef_emfac.html -diff --git a/man/ef_emfac.Rd b/man/ef_emfac.Rd deleted file mode 100644 index 869a5aa6..00000000 --- a/man/ef_emfac.Rd +++ /dev/null @@ -1,76 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/ef_emfac.R -\name{ef_emfac} -\alias{ef_emfac} -\title{EMFAC2017 emission factors for Statewide California, Calendar Year 2020} -\usage{ -ef_emfac( - veh = "LDT1", - fuel = "Gasoline", - mph, - pol = "CO_RUNEX", - season = "winter", - full = FALSE -) -} -\arguments{ -\item{veh}{Character; "one of the 40 vehicle categories shown below.} - -\item{fuel}{Character; "Diesel", "Gasoline", "Electricity" or "Natural Gas"} - -\item{mph}{Numeric; Speed in miles per hour (optional).} - -\item{pol}{Character; -\itemize{ -\item{if the user enter mph}{"NOx_RUNEX", "PM2.5_RUNEX", "PM10_RUNEX", -"CO2_RUNEX", "CH4_RUNEX", "N2O_RUNEX", "ROG_RUNEX", "TOG_RUNEX", "CO_RUNEX", -"SOx_RUNEX"} -\item{if the user do not enter mph}{ -"NOx_RUNEX", "NOx_IDLEX", "NOx_STREX", -"PM2.5_RUNEX", "PM2.5_IDLEX", "PM2.5_STREX", -"PM2.5_PMTW" , "PM2.5_PMBW", "PM10_RUNEX", -"PM10_IDLEX", "PM10_STREX", "PM10_PMTW", -"PM10_PMBW", "CO2_RUNEX", "CO2_IDLEX", -"CO2_STREX", "CH4_RUNEX", "CH4_IDLEX", -"CH4_STREX", "N2O_RUNEX", "N2O_IDLEX", -"N2O_STREX", "ROG_RUNEX", "ROG_IDLEX", -"ROG_STREX", "ROG_HOTSOAK", "ROG_RUNLOSS", -"ROG_RESTLOSS", "ROG_DIURN", "TOG_RUNEX", -"TOG_IDLEX", "TOG_STREX", "TOG_HOTSOAK", -"TOG_RUNLOSS", "TOG_RESTLOSS", "TOG_DIURN", -"CO_RUNEX", "CO_IDLEX", "CO_STREX", -"SOx_RUNEX", "SOx_IDLEX", "SOx_STREX" -} -}} - -\item{season}{Character: "winter" or "summer".} - -\item{full}{Logical: To return the whole data.table or not.} -} -\value{ -data.table with emission factors. -} -\description{ -\code{\link{ef_emfac}} returns emission factors reflecting California, US, -conditions. If the user enter speeds more emission factors are returned. -See details. -} -\examples{ -\dontrun{ -#do not run -pols <- c("CO_RUNEX", "NOx_RUNEX") -dfef2 <- ef_emfac(full = TRUE) -colplot(df = dfef2, - x = dfef2$Model_Year, - cols = pols, - main = "EF from LDT1 with Gasoline on Winter", - ylab = units(dfef2[[pols[1]]][1])) -} -} -\references{ -https://arb.ca.gov/emfac/emissions-inventory -} -\keyword{emfac} -\keyword{emission} -\keyword{factors} -\keyword{speed} diff --git a/man/get_project.Rd b/man/get_project.Rd index 4b1695a7..9e270fe7 100644 --- a/man/get_project.Rd +++ b/man/get_project.Rd @@ -23,9 +23,10 @@ get_project(directory, case = "brasil", url) amazon2014 \tab Top-down Amazon\tab CETESB+tunnel\tab csv and.rds\cr curitiba \tab Bottom-down +GTFS\tab CETESB+tunnel\tab csv and.rds\cr masp2020 \tab Bottom-down\tab CETESB+tunnel\tab csv and.rds\cr + ecuador_td \tab Top-down\tab EEA\tab csv and.rds\cr ecuador_td_hot \tab Top-down\tab EEA\tab csv and.rds\cr ecuador_td_hot_month \tab Top-down\tab EEA\tab csv and.rds\cr - moves_bu \tab Bottom-down\tab US/EPA MOVES \tab csv and.rds (requires MOVES on Windows)\cr + moves \tab Bottom-up\tab US/EPA MOVES \tab csv and.rds (requires MOVES >=3.0 on Windows)\cr }} \item{url}{String, with the URL to download VEIN project} @@ -37,6 +38,8 @@ The projects are available on Github.com/atmoschem/vein/projects \note{ default case can be any of "brasil", "brazil", "brazil_bu", "brasil_bu", they are the same +Projects for Ecuador are in development. +In any case, if you find any error, please, send a pull request in github or gitlab. } \examples{ \dontrun{ diff --git a/man/vein-deprecated.Rd b/man/vein-deprecated.Rd index 836eaad6..dbfd0373 100644 --- a/man/vein-deprecated.Rd +++ b/man/vein-deprecated.Rd @@ -1,6 +1,6 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/EmissionsList.R, R/Evaporative.R, R/emis_wrf.R, -% R/hot_soak.R, R/running_losses.R, R/vein-deprecated.R +% Please edit documentation in R/EmissionsList.R, R/Evaporative.R, R/ef_emfac.R, +% R/emis_wrf.R, R/hot_soak.R, R/running_losses.R, R/vein-deprecated.R \name{vein-deprecated} \alias{vein-deprecated} \alias{EmissionsList} @@ -11,6 +11,7 @@ \alias{print.Evaporative} \alias{summary.Evaporative} \alias{plot.Evaporative} +\alias{ef_emfac} \alias{emis_wrf} \alias{hot_soak} \alias{running_losses} @@ -20,6 +21,8 @@ EmissionsList(x, ...) Evaporative(x, ...) +ef_emfac(x, ...) + emis_wrf(sdf, nr = 1, dmyhm, tz, crs = 4326, islist) hot_soak(x, carb, p, eshotc, eswarmc, eshotfi) @@ -49,11 +52,25 @@ transform the coordinates of the output} \item{islist}{logical value to indicate if sdf is a list or not} \item{object}{Object with class "Evaporative"} + +\item{veh}{Character; "one of the 40 vehicle categories shown below.} + +\item{fuel}{Character; "Diesel", "Gasoline", "Electricity" or "Natural Gas"} + +\item{mph}{Numeric; Speed in miles per hour (optional).} + +\item{pol}{Character;} + +\item{season}{Character: "winter" or "summer".} + +\item{full}{Logical: To return the whole data.table or not.} } \value{ Objects of class "EmissionsList" and numeric elements as "units" Objects of class "Evaporative" or "units" + +data.table with emission factors. } \description{ \code{EmissionsList} returns a tranformed object with class "EmissionsList". @@ -63,6 +80,10 @@ units g. This class represents the daily emissions presented by Mellios G and Ntziachristos (2016) Gasoline evaporation, Tier 2. Eventually it will be incorporated the techniques of Tier 3. +\code{\link{ef_emfac}} returns emission factors reflecting California, US, +conditions. If the user enter speeds more emission factors are returned. +See details. + \code{emis_wrf} returns a dataframes with columns lat, long, id, pollutants, local time and GMT time. This dataframe has the proper format to be used with WRF assimilation system: "ASimilation System 4 WRF (AS4WRF Vera-Vala et al (2016)) @@ -94,6 +115,10 @@ For \code{running_losses}, use \code{\link{emis_evap}}. } \dontrun{ # do not run +# DEPRECATED +} +\dontrun{ +# do not run # Use library(eixport) # eixport::to_as4wrf() } @@ -106,6 +131,8 @@ For \code{running_losses}, use \code{\link{emis_evap}}. } \references{ +https://arb.ca.gov/emfac/emissions-inventory + Mellios G and Ntziachristos 2016. Gasoline evaporation. In: EEA, EMEP. EEA air pollutant emission inventory guidebook-2009. European Environment Agency, Copenhagen, 2009 @@ -119,6 +146,8 @@ Environment Agency, Copenhagen, 2009 \code{\link{vein-deprecated}} +\code{\link{vein-deprecated}} + \code{\link{vein-deprecated}} } \keyword{internal} diff --git a/projects/ecuador_td_hot_country.tar.gz b/projects/ecuador_td.tar.gz similarity index 76% rename from projects/ecuador_td_hot_country.tar.gz rename to projects/ecuador_td.tar.gz index 435efa30..4c488238 100644 Binary files a/projects/ecuador_td_hot_country.tar.gz and b/projects/ecuador_td.tar.gz differ diff --git a/projects/ecuador_td_hot_country/config/clean.R b/projects/ecuador_td/config/clean.R similarity index 100% rename from projects/ecuador_td_hot_country/config/clean.R rename to projects/ecuador_td/config/clean.R diff --git a/projects/ecuador_td_hot_country/config/config.R b/projects/ecuador_td/config/config.R similarity index 99% rename from projects/ecuador_td_hot_country/config/config.R rename to projects/ecuador_td/config/config.R index ae57321b..70dda5c4 100755 --- a/projects/ecuador_td_hot_country/config/config.R +++ b/projects/ecuador_td/config/config.R @@ -130,6 +130,7 @@ if (!"region" %in% names(fuel)) { setDT(fuel) fuel <- fuel[Year == year] #all the months! and all the provinces +fuel$id <- 1:nrow(fuel) saveRDS(metadata, "config/metadata.rds") saveRDS(mileage, "config/mileage.rds") diff --git a/projects/ecuador_td/config/inventory_ecuador.xlsx b/projects/ecuador_td/config/inventory_ecuador.xlsx new file mode 100755 index 00000000..0d9c93e3 Binary files /dev/null and b/projects/ecuador_td/config/inventory_ecuador.xlsx differ diff --git a/projects/ecuador_td_hot_country/main.R b/projects/ecuador_td/main.R similarity index 100% rename from projects/ecuador_td_hot_country/main.R rename to projects/ecuador_td/main.R diff --git a/projects/ecuador_td_hot_country/main.Rproj b/projects/ecuador_td/main.Rproj similarity index 100% rename from projects/ecuador_td_hot_country/main.Rproj rename to projects/ecuador_td/main.Rproj diff --git a/projects/ecuador_td_hot_country/network/ecuador_provincias.gpkg b/projects/ecuador_td/network/ecuador_provincias.gpkg similarity index 100% rename from projects/ecuador_td_hot_country/network/ecuador_provincias.gpkg rename to projects/ecuador_td/network/ecuador_provincias.gpkg diff --git a/projects/ecuador_td_hot_country/network/lights.tif b/projects/ecuador_td/network/lights.tif similarity index 100% rename from projects/ecuador_td_hot_country/network/lights.tif rename to projects/ecuador_td/network/lights.tif diff --git a/projects/ecuador_td_hot_country/network/net.gpkg b/projects/ecuador_td/network/net.gpkg similarity index 100% rename from projects/ecuador_td_hot_country/network/net.gpkg rename to projects/ecuador_td/network/net.gpkg diff --git a/projects/ecuador_td_hot_country/network/net.rds b/projects/ecuador_td/network/net.rds similarity index 100% rename from projects/ecuador_td_hot_country/network/net.rds rename to projects/ecuador_td/network/net.rds diff --git a/projects/ecuador_td_hot_country/network/roads.rds b/projects/ecuador_td/network/roads.rds similarity index 100% rename from projects/ecuador_td_hot_country/network/roads.rds rename to projects/ecuador_td/network/roads.rds diff --git a/projects/ecuador_td_hot_country/notes/README_2021-09-27_225047.txt b/projects/ecuador_td/notes/README_2021-09-29_221845.txt similarity index 86% rename from projects/ecuador_td_hot_country/notes/README_2021-09-27_225047.txt rename to projects/ecuador_td/notes/README_2021-09-29_221845.txt index 56284915..28ba4d18 100644 --- a/projects/ecuador_td_hot_country/notes/README_2021-09-27_225047.txt +++ b/projects/ecuador_td/notes/README_2021-09-29_221845.txt @@ -4,13 +4,13 @@ Ecuador province 2019 Directory: /home/sergio/models/vein/projects/ecuador_td_hot_country -Local Time: 2021-09-27 22:50:47 +Local Time: 2021-09-29 22:18:46 Inventory compiler: sergio ======================================== sysname = Linux -release = 5.11.0-36-generic -version = #40-Ubuntu SMP Fri Sep 17 18:15:22 UTC 2021 +release = 5.11.0-37-generic +version = #41-Ubuntu SMP Mon Sep 20 16:39:20 UTC 2021 nodename = sergio-pc machine = x86_64 user = sergio @@ -79,9 +79,10 @@ loaded via a namespace (and not attached): [19] fields_12.5 dotCall64_1.0-1 munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0 raster_3.4-10 [25] codetools_0.2-18 labeling_0.4.2 parallel_4.0.4 class_7.3-18 fansi_0.5.0 Rcpp_1.0.7 [31] KernSmooth_2.23-18 scales_1.1.1 classInt_0.4-3 farver_2.1.0 gridExtra_2.3 digest_0.6.27 -[37] dplyr_1.0.7 ncdf4_1.17 grid_4.0.4 rgdal_1.5-23 tools_4.0.4 magrittr_2.0.1 -[43] maps_3.4.0 proxy_0.4-26 tibble_3.1.2 crayon_1.4.1 pkgconfig_2.0.3 ellipsis_0.3.2 -[49] assertthat_0.2.1 viridis_0.6.1 R6_2.5.0 units_0.7-2 compiler_4.0.4 +[37] dplyr_1.0.7 ncdf4_1.17 grid_4.0.4 rgdal_1.5-23 cli_3.0.0 tools_4.0.4 +[43] magrittr_2.0.1 maps_3.4.0 proxy_0.4-26 tibble_3.1.2 crayon_1.4.1 pkgconfig_2.0.3 +[49] ellipsis_0.3.2 assertthat_0.2.1 rstudioapi_0.13 viridis_0.6.1 R6_2.5.0 units_0.7-2 +[55] compiler_4.0.4 ======================================== diff --git a/projects/ecuador_td_hot_country/scripts/ef_eea.R b/projects/ecuador_td/scripts/ef_eea.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/ef_eea.R rename to projects/ecuador_td/scripts/ef_eea.R diff --git a/projects/ecuador_td_hot_country/scripts/evaporatives.R b/projects/ecuador_td/scripts/evaporatives.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/evaporatives.R rename to projects/ecuador_td/scripts/evaporatives.R diff --git a/projects/ecuador_td_hot_country/scripts/exhaust.R b/projects/ecuador_td/scripts/exhaust.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/exhaust.R rename to projects/ecuador_td/scripts/exhaust.R diff --git a/projects/ecuador_td_hot_country/scripts/exhaust_eea.R b/projects/ecuador_td/scripts/exhaust_eea.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/exhaust_eea.R rename to projects/ecuador_td/scripts/exhaust_eea.R diff --git a/projects/ecuador_td_hot_country/scripts/fuel_eval.R b/projects/ecuador_td/scripts/fuel_eval.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/fuel_eval.R rename to projects/ecuador_td/scripts/fuel_eval.R diff --git a/projects/ecuador_td_hot_country/scripts/fuel_eval_eea.R b/projects/ecuador_td/scripts/fuel_eval_eea.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/fuel_eval_eea.R rename to projects/ecuador_td/scripts/fuel_eval_eea.R diff --git a/projects/ecuador_td_hot_country/scripts/fuel_eval_eea_year.R b/projects/ecuador_td/scripts/fuel_eval_eea_year.R similarity index 70% rename from projects/ecuador_td_hot_country/scripts/fuel_eval_eea_year.R rename to projects/ecuador_td/scripts/fuel_eval_eea_year.R index 190b100b..5b205434 100644 --- a/projects/ecuador_td_hot_country/scripts/fuel_eval_eea_year.R +++ b/projects/ecuador_td/scripts/fuel_eval_eea_year.R @@ -14,76 +14,77 @@ for (i in seq_along(metadata$vehicles)) { "\n", metadata$vehicles[i], rep("", max(nchar(metadata$vehicles) + 1) - nchar(metadata$vehicles[i])) ) - + x <- readRDS(paste0("veh/", metadata$vehicles[i], ".rds")) - + # euro cate <- as.character(as.roman(gsub("Euro ", "", euro[[metadata$vehicles[i]]]))) cate[is.na(cate)] <- "PRE" for (j in seq_along(pol)) { cat(pol[j], " ") - + if(metadata$v_eea_old[i] %in% c("PC", "LCV", "Motorcycle")) { - + rbindlist(lapply(seq_along(net[[id_region]]), function(k) { - ef_ldv_speed(v = metadata$v_eea_old[i], - t = metadata$t_eea_old[i], - cc = metadata$cc_eea_old[i], - f = metadata$fuel_eea_old[i], - p = "FC", - eu = cate, - speed = Speed(rep(metadata$speed[i], 12))) + ef_ldv_speed(v = metadata$v_eea_old[i], + t = metadata$t_eea_old[i], + cc = metadata$cc_eea_old[i], + f = metadata$fuel_eea_old[i], + p = "FC", + eu = cate, + speed = Speed(rep(metadata$speed[i], 12))) }))-> ef } else { rbindlist(lapply(seq_along(net[[id_region]]), function(k) { - ef <- ef_hdv_speed(v = metadata$v_eea_old[i], - t = metadata$t_eea_old[i], - g = metadata$cc_eea_old[i], - eu = cate, - gr = 0, - l = 0.5, - p = "FC", - speed = Speed(rep(metadata$speed[i], 12))) - - }))-> ef + ef <- ef_hdv_speed(v = metadata$v_eea_old[i], + t = metadata$t_eea_old[i], + g = metadata$cc_eea_old[i], + eu = cate, + gr = 0, + l = 0.5, + p = "FC", + speed = Speed(rep(metadata$speed[i], 12))) + + }))-> ef } nrow(x) == nrow(ef) ef$speed <- NULL - array_x <- emis_hot_td( - veh = x, - lkm = mileage[[metadata$vehicles[i]]], - ef = ef[, 1:ncol(x)], - fortran = TRUE, - nt = check_nt() / 2, - # pro_month = pmonth[[metadata$vehicles[i]]], - verbose = verbose, - params = list( - veh = metadata$vehicles[i], - size = metadata$size[i], - fuel = metadata$fuel[i], - pollutant = pol[j], - type_emi = "Exhaust", - subtype_emi = "Exhaust", - baseyear = year, - month = 1:12 - ) - ) - array_x$region <- rep(net[[id_region]], each = 12) + eff <- ef[, 1:ncol(x)] + lkm <- mileage[[metadata$vehicles[i]]] - fwrite(array_x, "emi/FC.csv", append = TRUE) + rbindlist(lapply(1:ncol(x), function(k) { + emi = (as.numeric(x[[k]])*eff[[k]])*lkm[k] + + dt <- data.table(emissions = emi, + rows = seq_along(emi)) + dt$age <- k + dt$veh = metadata$vehicles[i] + dt$size = metadata$size[i] + dt$fuel = metadata$fuel[i] + dt$pollutant = pol[j] + type_emi = "Exhaust" + subtype_emi = "Exhaust" + baseyear = year + dt$month <- rep(1:12, each = length(unique(net[[id_region]]))) + dt$region <- rep(net[[id_region]], 12) + dt + })) -> array_x + + fwrite(array_x, "emi/FC.csv", append = TRUE) + } } switch(language, - "portuguese" = message("\nArquivos em:"), - "english" = message("\nFiles in:"), - "spanish" = message("\nArchivos en:") + "portuguese" = message("\nArquivos em:"), + "english" = message("\nFiles in:"), + "spanish" = message("\nArchivos en:") ) cat(paste0(getwd(), "/emi/*\n")) @@ -98,8 +99,8 @@ dt$t <- units::set_units(dt$g, t) dt$fuel_month_region <- paste0(dt$fuel, "_", dt$month, "_", dt$region) dt0 <- dt[pollutant == "FC", - round(sum(t), 2), - by = .(fuel_month_region) + round(sum(t), 2), + by = .(fuel_month_region) ] data.table::setkey(dt0, "fuel_month_region") @@ -116,7 +117,12 @@ dtf$consumption_lt <- units::set_units(dtf$consumption_lt, "l") dtf$consumption_m3 <- units::set_units(dtf$consumption_lt, "m^3") dtf$consumption_t <- dtf$consumption_m3 * dtf$density_tm3 dtf$estimation_consumption <- dtf$estimation_t / dtf$consumption_t -print(dtf[, c("region", "Month", "fuel", "estimation_t", "consumption_t", "estimation_consumption")]) +print(dtf[, c("region", + "Month", + "fuel", + "estimation_t", + "consumption_t", + "estimation_consumption")]) # calibrate k #### k_D <- as.numeric(1/dtf[fuel == "D"]$estimation_consumption) @@ -190,29 +196,30 @@ for (i in seq_along(metadata$vehicles)) { nrow(x) == nrow(ef) ef$speed <- NULL + eff <- ef[, 1:ncol(x)] + lkm <- mileage[[metadata$vehicles[i]]] + - array_x <- emis_hot_td( - veh = x, - lkm = mileage[[metadata$vehicles[i]]], - ef = ef[, 1:ncol(x)], - fortran = TRUE, - nt = check_nt() / 2, - # pro_month = pmonth[[metadata$vehicles[i]]], - verbose = verbose, - params = list( - veh = metadata$vehicles[i], - size = metadata$size[i], - fuel = metadata$fuel[i], - pollutant = pol[j], - type_emi = "Exhaust", - subtype_emi = "Exhaust", - baseyear = year, - month = 1:12 - ) - ) - array_x$region <- rep(net[[id_region]], each = 12) + rbindlist(lapply(1:ncol(x), function(k) { + emi = (as.numeric(x[[k]])*eff[[k]])*lkm[k] + + dt <- data.table(emissions = emi, + rows = seq_along(emi)) + dt$age <- k + dt$veh = metadata$vehicles[i] + dt$size = metadata$size[i] + dt$fuel = metadata$fuel[i] + dt$pollutant = pol[j] + type_emi = "Exhaust" + subtype_emi = "Exhaust" + baseyear = year + dt$month <- rep(1:12, each = length(unique(net[[id_region]]))) + dt$region <- rep(net[[id_region]], 12) + dt + })) -> array_x fwrite(array_x, "emi/FC.csv", append = TRUE) + } } @@ -255,9 +262,9 @@ dtf$estimation_consumption <- dtf$estimation_t / dtf$consumption_t print(dtf[, c("region", "Month", "fuel", "estimation_t", "consumption_t", "estimation_consumption")]) switch(language, - "portuguese" = message("Limpando..."), - "english" = message("Cleaning..."), - "spanish" = message("Limpiando...") + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") ) # suppressWarnings(rm(i, j, pol, dt, dt0, dtf, factor_emi, fuel)) # suppressWarnings(rm( @@ -275,4 +282,4 @@ switch(language, # k_G, k_D, k_E, # )) ls() -invisible(gc()) \ No newline at end of file +invisible(gc()) diff --git a/projects/ecuador_td_hot_country/scripts/mech.R b/projects/ecuador_td/scripts/mech.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/mech.R rename to projects/ecuador_td/scripts/mech.R diff --git a/projects/ecuador_td_hot_country/scripts/mech2.R b/projects/ecuador_td/scripts/mech2.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/mech2.R rename to projects/ecuador_td/scripts/mech2.R diff --git a/projects/ecuador_td_hot_country/scripts/net.R b/projects/ecuador_td/scripts/net.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/net.R rename to projects/ecuador_td/scripts/net.R diff --git a/projects/ecuador_td_hot_country/scripts/osm.R b/projects/ecuador_td/scripts/osm.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/osm.R rename to projects/ecuador_td/scripts/osm.R diff --git a/projects/ecuador_td_hot_country/scripts/plots.R b/projects/ecuador_td/scripts/plots.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/plots.R rename to projects/ecuador_td/scripts/plots.R diff --git a/projects/ecuador_td_hot_country/scripts/post.R b/projects/ecuador_td/scripts/post.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/post.R rename to projects/ecuador_td/scripts/post.R diff --git a/projects/ecuador_td_hot_country/scripts/ressuspensao.R b/projects/ecuador_td/scripts/ressuspensao.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/ressuspensao.R rename to projects/ecuador_td/scripts/ressuspensao.R diff --git a/projects/ecuador_td_hot_country/scripts/traffic.R b/projects/ecuador_td/scripts/traffic.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/traffic.R rename to projects/ecuador_td/scripts/traffic.R diff --git a/projects/ecuador_td_hot_country/scripts/traffic_month.R b/projects/ecuador_td/scripts/traffic_month.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/traffic_month.R rename to projects/ecuador_td/scripts/traffic_month.R diff --git a/projects/ecuador_td_hot_country/scripts/traffic_year.R b/projects/ecuador_td/scripts/traffic_year.R similarity index 99% rename from projects/ecuador_td_hot_country/scripts/traffic_year.R rename to projects/ecuador_td/scripts/traffic_year.R index 04e8b33c..b1764875 100644 --- a/projects/ecuador_td_hot_country/scripts/traffic_year.R +++ b/projects/ecuador_td/scripts/traffic_year.R @@ -53,6 +53,7 @@ n_BUS <- nveh[grep(pattern = "BUS", x = nveh)] n_MC <- nveh[grep(pattern = "MC", x = nveh)] n_PC <- c(n_PC, n_TAXI) + # PC and TAXI #### # apply survival functions for (i in seq_along(metadata$vehicles)) { diff --git a/projects/ecuador_td_hot_country/scripts/wrf.R b/projects/ecuador_td/scripts/wrf.R similarity index 100% rename from projects/ecuador_td_hot_country/scripts/wrf.R rename to projects/ecuador_td/scripts/wrf.R diff --git a/projects/ecuador_td_hot.tar.gz b/projects/ecuador_td_hot.tar.gz old mode 100755 new mode 100644 index 5d76569c..cc595f1f Binary files a/projects/ecuador_td_hot.tar.gz and b/projects/ecuador_td_hot.tar.gz differ diff --git a/projects/ecuador_td_hot_country/config/inventory_ecuador.xlsx b/projects/ecuador_td_hot_country/config/inventory_ecuador.xlsx deleted file mode 100755 index 2d7b89be..00000000 Binary files a/projects/ecuador_td_hot_country/config/inventory_ecuador.xlsx and /dev/null differ diff --git a/projects/ecuador_td_hot_month.tar.gz b/projects/ecuador_td_hot_month.tar.gz index e5e6dc0f..4db9fc9a 100644 Binary files a/projects/ecuador_td_hot_month.tar.gz and b/projects/ecuador_td_hot_month.tar.gz differ diff --git a/projects/moves.tar.gz b/projects/moves.tar.gz new file mode 100644 index 00000000..21aa28b7 Binary files /dev/null and b/projects/moves.tar.gz differ diff --git a/projects/moves/.vscode/settings.json b/projects/moves/.vscode/settings.json new file mode 100755 index 00000000..e3bf952e --- /dev/null +++ b/projects/moves/.vscode/settings.json @@ -0,0 +1,9 @@ +{ + "workbench.settings.openDefaultKeybindings": true, + "terminal.integrated.commandsToSkipShell": [ + "{ \"key\": \"alt+x\", \"command\": \"workbench.action.terminal.sendSequence\", \"args\": { \"text\": \"\\u001b[A\\u000d\" } }," + ], + "terminal.integrated.automationShell.windows": "", + "terminal.integrated.experimentalLinkProvider": false, + "terminal.integrated.sendKeybindingsToShell": true +} \ No newline at end of file diff --git a/projects/moves/config/clean.R b/projects/moves/config/clean.R new file mode 100755 index 00000000..e6ae64c0 --- /dev/null +++ b/projects/moves/config/clean.R @@ -0,0 +1,16 @@ +# apagando dados +a <- list.files(path = "config", pattern = ".rds", full.names = T) +file.remove(a) + +unlink("csv", recursive = T) +unlink("emi", recursive = T) +unlink("images", recursive = T) +unlink("notes", recursive = T) +unlink("post", recursive = T) +unlink("wrf/wrfc*") +unlink("veh", recursive = T) +file.remove(list.files("network/", ".rds", full.names = T)) + +system(paste0("tar -caf ", basename(getwd()), ".tar.gz .")) +system(paste0("mv ", basename(getwd()), ".tar.gz ../")) +file.remove(".Rhistory") diff --git a/projects/moves/config/config.R b/projects/moves/config/config.R new file mode 100755 index 00000000..4f78ab28 --- /dev/null +++ b/projects/moves/config/config.R @@ -0,0 +1,457 @@ +# apagando dados #### +a <- list.files(path = "config", pattern = ".rds", full.names = T) +file.remove(a) + +# configuracao #### +metadata <- as.data.frame(metadata) +mileage <- as.data.frame(mileage) +mileage[, metadata$vehicles] <- add_miles(mileage[, metadata$vehicles]) +tfs <- as.data.frame(tfs) +veh <- as.data.frame(veh) +vmt_age <- as.data.frame(vmt_age) +fuel <- as.data.frame(fuel) +fuel_type <- as.data.frame(fuel_type) +met <- as.data.frame(met) + +setDT(met) +setDT(fuel) + +# checkar metadata$vehicles #### +switch(language, + "portuguese" = cat("Metadata$Vehicles é:\n"), + "english" = cat("Metadata$Vehicles is:\n"), + "spanish" = cat("Metadata$Vehicles es:\n") +) + +# cat( "Metadata$Vehicles é:\n") +print(metadata$vehicles) + +# checar nomes mileage #### +if (!length(intersect(metadata$vehicles, names(mileage))) == length(metadata$vehicles)) { + switch(language, + "portuguese" = stop( + "Precisa adicionar coluna ", + setdiff(metadata$vehicles, names(mileage)), + " em `mileage`" + ), + "english" = stop( + "You need to add column ", + setdiff(metadata$vehicles, names(mileage)), + " in `mileage`" + ), + "spanish" = stop( + "Necesitas agregar la columna ", + setdiff(metadata$vehicles, names(mileage)), + " en `mileage`" + ) + ) +} + +# checar nomes tfs #### +if (!length(intersect(metadata$vehicles, names(tfs))) == length(metadata$vehicles)) { + switch(language, + "portuguese" = stop( + "Precisa adicionar coluna ", + setdiff(metadata$vehicles, names(mileage)), + " em `tfs`" + ), + "english" = stop( + "You need to add column ", + setdiff(metadata$vehicles, names(mileage)), + " in `tfs`" + ), + "spanish" = stop( + "Necesitas agregar la columna ", + setdiff(metadata$vehicles, names(mileage)), + " en `tfs`" + ) + ) +} + +# checar nomes veh #### +if (!length(intersect(metadata$vehicles, names(veh))) == length(metadata$vehicles)) { + switch(language, + "portuguese" = stop( + "Precisa adicionar coluna ", + setdiff(metadata$vehicles, names(mileage)), + " em `veh`" + ), + "english" = stop( + "You need to add column ", + setdiff(metadata$vehicles, names(mileage)), + " in `veh`" + ), + "spanish" = stop( + "Necesitas agregar la columna ", + setdiff(metadata$vehicles, names(mileage)), + " en `veh`" + ) + ) +} + +# checar Year #### +if (!"Year" %in% names(veh)) { + switch(language, + "portuguese" = stop("Não estou enxergando a coluna 'Year' em `veh`"), + "english" = stop("I'm not seeing column 'Year' in `veh`"), + "spanish" = stop("No estoy viendo la columna 'Year' in `veh`") + ) +} + + +# veh +veh <- veh[veh$Year <= year, ][1:31,] +dim(veh) + +veh[is.na(veh)] <- 0 + +switch(language, + "portuguese" = message("Arquivos em: ", getwd(), "/config/*\n"), + "english" = message("Files in: ", getwd(), "/config/*\n"), + "spanish" = message("Archivos en: ", getwd(), "/config/*\n") +) +#fuel +names(fuel) +head(fuel) +fuel <- fuel[Year == year & Month == month] + + +# met #### + +met <- met[Year == year & Month == month] + +saveRDS(metadata, "config/metadata.rds") +saveRDS(mileage, "config/mileage.rds") +saveRDS(tfs, "config/tfs.rds") +saveRDS(veh, "config/fleet_age.rds") +saveRDS(fuel, "config/fuel.rds") +saveRDS(fuel_type, "config/fuel_type.rds") +saveRDS(met, "config/met.rds") + +# pastas +if (delete_directories) { + choice <- 1 + + if (language == "portuguese") { + # choice <- utils::menu(c("Sim", "Não"), title="Apagar pastas csv, emi, images, notes, post e veh??") + if (choice == 1) { + message("Apagando pastas `emi`, `images`, `notes`, `post` e `veh`") + unlink("csv", recursive = T) + unlink("emi", recursive = T) + unlink("images", recursive = T) + unlink("notes", recursive = T) + unlink("post", recursive = T) + unlink("veh", recursive = T) + } + } else if (language == "english") { + # choice <- utils::menu(c("Yes", "No"), title="Delete folders `csv`, `emi`, `images`, `notes`, `post` e `veh`??") + if (choice == 1) { + message("Deleting folders `emi`, `images`, `notes`, `post` and `veh`") + unlink("csv", recursive = T) + unlink("emi", recursive = T) + unlink("images", recursive = T) + unlink("notes", recursive = T) + unlink("post", recursive = T) + unlink("veh", recursive = T) + } + } else if (language == "spanish") { + # choice <- utils::menu(c("Si", "No"), title="Borrar carpetas `csv`, `emi`, `images`, `notes`, `post` y `veh`??") + if (choice == 1) { + message("Borrando carpetas `emi`, `images`, `notes`, `post` y `veh`") + unlink("csv", recursive = T) + unlink("emi", recursive = T) + unlink("notes", recursive = T) + unlink("images", recursive = T) + unlink("post", recursive = T) + unlink("veh", recursive = T) + } + } +} + +dir.create(path = "csv", showWarnings = FALSE) +dir.create(path = "emi", showWarnings = FALSE) +dir.create(path = "emi/FC", showWarnings = FALSE) +dir.create(path = "emi/emissions_rate_per_distance", showWarnings = FALSE) +dir.create(path = "emi/emissions_rate_per_vehicle", showWarnings = FALSE) + +dir.create(path = "images", showWarnings = FALSE) +dir.create(path = "notes", showWarnings = FALSE) +dir.create(path = "post", showWarnings = FALSE) +dir.create(path = "post/datatable", showWarnings = FALSE) +dir.create(path = "post/streets", showWarnings = FALSE) +dir.create(path = "post/grids", showWarnings = FALSE) +dir.create(path = "veh", showWarnings = FALSE) + +# for (i in seq_along(metadata$vehicles)) dir.create(path = paste0("emi/", metadata$vehicles[i])) + + +pa <- list.dirs(path = "emi", full.names = T, recursive = T) +po <- list.dirs("post", full.names = T, recursive = T) + +switch(language, + "portuguese" = message("Novas pastas:"), + "english" = message("New directories:"), + "spanish" = message("Nuevas carpetas") +) + +message("csv\n") +message("images\n") +message(paste0(po, "\n")) +message(paste0(pa, "\n")) +message("veh\n") + +# names groups #### +n_PC <- metadata$vehicles[grep(pattern = "PC", x = metadata$vehicles)] +n_PT <- metadata$vehicles[grep(pattern = "PT", x = metadata$vehicles)] +n_LCT <- metadata$vehicles[grep(pattern = "LCT", x = metadata$vehicles)] +n_TRUCKS <- metadata$vehicles[grep(pattern = "TRUCKS", x = metadata$vehicles)] +n_BUS <- metadata$vehicles[grep(pattern = "BUS", x = metadata$vehicles)] +n_MC <- metadata$vehicles[grep(pattern = "MC", x = metadata$vehicles)] +n_veh <- list( + PC = n_PC, + PT = n_PT, + LCT = n_LCT, + TRUCKS = n_TRUCKS, + BUS = n_BUS, + MC = n_MC +) +# Fuel #### +switch(language, + "portuguese" = cat("Plotando combustivel \n"), + "english" = cat("Plotting fuel \n"), + "spanish" = cat("Plotando combustible \n") +) + +png("images/FUEL.png", width = 1500, height = 2000, units = "px", res = 300) +barplot( + height = fuel$gallons, + names.arg = fuel$fuel, xlab = "Fuel", + ylab = "gallons", + main = paste0("Fuel consumded during Year ", year, " and month ", month) +) +month +dev.off() + +setDF(veh) +# Fleet #### +switch(language, + "portuguese" = cat("Plotando frota \n"), + "english" = cat("Plotting fleet \n"), + "spanish" = cat("Plotando flota \n") +) + +for (i in seq_along(n_veh)) { + df_x <- as.data.frame(veh[, n_veh[[i]]]) + if(ncol(df_x) == 1) { + names(df_x) <- n_veh[i][[1]] + } + png( + paste0( + "images/FLEET_", + names(n_veh)[i], + ".png" + ), + 2000, 1500, "px", + res = 300 + ) + colplot( + df = df_x, + cols = n_veh[[i]], + xlab = "Age", + ylab = "veh", + main = names(n_veh)[i], + type = "l", + pch = NULL, + lwd = 1, + theme = theme, + spl = 8 + ) + dev.off() +} + +# Plot TFS #### +switch(language, + "portuguese" = cat("Plotando perfis `tfs`\n"), + "english" = cat("Plotting profiles `tfs`\n"), + "spanish" = cat("Plotando perfiles `tfs`\n") +) + +for (i in seq_along(n_veh)) { + df_x <- as.data.frame(tfs[, n_veh[[i]]]) + if(ncol(df_x) == 1) { + names(df_x) <- n_veh[i][[1]] + } + png( + paste0( + "images/TFS_", + names(n_veh)[i], + ".png" + ), + 2000, 1500, "px", + res = 300 + ) + colplot( + df = df_x, + cols = n_veh[[i]], + xlab = "Hour", + ylab = "", + main = paste0("TFS ", names(n_veh)[i]), + type = "l", + pch = NULL, + lwd = 1, + theme = theme, + spl = 8 + ) + dev.off() +} + + +# Plot Mileage #### +switch(language, + "portuguese" = cat("Plotando quilometragem \n"), + "english" = cat("Plotting mileage `tfs`\n"), + "spanish" = cat("Plotando kilometraje `tfs`\n") +) + +for (i in seq_along(n_veh)) { + + df_x <- as.data.frame(mileage[, n_veh[[i]]]) + if(ncol(df_x) == 1) { + names(df_x) <- n_veh[i][[1]] + } + png( + paste0( + "images/MILEAGE_", + names(n_veh)[i], + ".png" + ), + 2000, 1500, "px", + res = 300 + ) + colplot( + df = df_x, + cols = n_veh[[i]], + xlab = "Age of use", + ylab = "[km/year]", + main = paste0("Mileage ", names(n_veh)[i]), + type = "l", + pch = NULL, + lwd = 1, + theme = theme, + spl = 8 + ) + dev.off() +} + + + +# Plot vmt_age #### +switch(language, + "portuguese" = cat("Plotando vmt \n"), + "english" = cat("Plotting vmt \n"), + "spanish" = cat("Plotando vmt \n") +) + +for (i in seq_along(n_veh)) { + + df_x <- as.data.frame(vmt_age[, n_veh[[i]]]) + if(ncol(df_x) == 1) { + names(df_x) <- n_veh[i][[1]] + } + png( + paste0( + "images/VMT_AGE_", + names(n_veh)[i], + ".png" + ), + 2000, 1500, "px", + res = 300 + ) + colplot( + df = df_x, + cols = n_veh[[i]], + xlab = "Age of use", + ylab = "[km/year]", + main = paste0("VMT ", names(n_veh)[i]), + type = "l", + pch = NULL, + lwd = 1, + theme = theme, + spl = 8 + ) + dev.off() +} + + +# Temperature #### +png("images/Temperature.png", + 2000, 1500, "px", + res = 300 +) +colplot( + df = met, + cols = "Temperature", + xlab = "Months", + ylab = "°F", + main = "Temperature", + type = "l", + pch = NULL, + lwd = 1, + theme = theme, + spl = 8 +) +dev.off() + +# Notes #### +switch(language, + "portuguese" = cat("\nFazendo anotações\n"), + "english" = cat("\nTaking some notes\n"), + "spanish" = cat("\nEscribiendo notas\n") +) + +vein_notes( + notes = c("Default notes for vein::get_project"), + file = "notes/README", + title = paste0(basename(getwd()), year), + approach = "Bottom-up", + traffic = "AADT", + composition = "MOVES", + ef = "MOVES ", + cold_start = "Not Applicable", + evaporative = "all", + standards = "US/EPA", + mileage = "US/EPA" +) + +# message #### +switch(language, + "portuguese" = message("\nArquivos em:"), + "english" = message("\nFiles in:"), + "spanish" = message("\nArchivos en:") +) + +message( + "config/metadata.rds\n", + "config/tfs.rds\n", + "config/fleet_age.rds\n", + "config/fuel.rds\n" +) + +switch(language, + "portuguese" = message("\nFiguras em /images\n"), + "english" = message("\nFigures in /image\n"), + "spanish" = message("\nFiguras en /images\n") +) + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + + +suppressWarnings( + rm(i, choice, pa, metadata, po, tfs, veh, mileage, fuel, theme, + n_PC, n_LCV, n_TRUCKS, n_BUS, n_MC, df_x, ef, cores, vkm, ef, a, rota, + delete_directories, met, month, n_veh, scale, path, year) +) \ No newline at end of file diff --git a/projects/moves/config/inventory_MD.xlsx b/projects/moves/config/inventory_MD.xlsx new file mode 100755 index 00000000..1f13d596 Binary files /dev/null and b/projects/moves/config/inventory_MD.xlsx differ diff --git a/projects/moves/config/veh_fuel.R b/projects/moves/config/veh_fuel.R new file mode 100755 index 00000000..fb7c3205 --- /dev/null +++ b/projects/moves/config/veh_fuel.R @@ -0,0 +1,39 @@ +library(data.table) +fuel <- readxl::read_excel(path = "../USdata/annual-motor-vehicle-population-by-type-of-fuel-used.xlsx") + +unique(fuel$type) +fuel$veh <- ifelse( + fuel$type == "Cars", "LDV", + ifelse( + fuel$type == "Goods and Other Vehicles", "HDV", + ifelse( + fuel$type == "Motorcycles", "MC", + NA))) + +unique(fuel$engine) +# https://data.gov.sg/dataset/annual-motor-vehicle-population-by-type-of-fuel-used?view_id=623935d0-3584-43b1-9a03-974b21a77db9&resource_id=39e17d4c-3bd3-4814-9ae7-bf9df39a2f5d +fuel$fuel <- ifelse( + fuel$engine == "Petrol", "G", + ifelse( + fuel$engine == "Diesel", "D", + ifelse( + fuel$engine %in% c("Petrol-Electric", "Petrol-Electric (Plug-In)", "Electric", "Diesel-Electric"), "ELEC", + ifelse( + fuel$engine %in% c("Petrol-CNG", "CNG"), "CNG", + NA)))) +fuel +setDT(fuel) +fuel[!is.na(veh) & + !is.na(fuel), + sum(number), + by = .(year, veh, fuel)] -> xx +setorderv(xx, cols = "year", order = -1) + +dcast.data.table(data = xx, + formula = year ~ veh+fuel, + value.var = "V1") -> dx +setorderv(dx, cols = "year", order = -1) +writexl::write_xlsx(dx, "config/veh_fuel.xlsx") +# edit inventory.xlsx +rm(list = ls()) +gc() diff --git a/projects/moves/config/veh_fuel.xlsx b/projects/moves/config/veh_fuel.xlsx new file mode 100755 index 00000000..2656a905 Binary files /dev/null and b/projects/moves/config/veh_fuel.xlsx differ diff --git a/projects/moves/main.R b/projects/moves/main.R new file mode 100755 index 00000000..37c9c381 --- /dev/null +++ b/projects/moves/main.R @@ -0,0 +1,220 @@ +options(encoding = "UTF-8") +library(vein) # vein +library(sf) # spatial data +library(cptcity) # 7120 colour palettes +library(ggplot2) # plots +library(ggthemes) +library(eixport) # WRF Chem +library(data.table) # blasting speed +sessionInfo() + + +# 0 Configuration +language <- "english" # portuguese spanishx'x' +path <- "config/inventory_MD.xlsx" +readxl::excel_sheets(path) +metadata <- readxl::read_xlsx(path = path, sheet = "metadata") +mileage <- readxl::read_xlsx(path = path, sheet = "mileage") +tfs <- readxl::read_xlsx(path = path, sheet = "tfs_moves") +veh <- readxl::read_xlsx(path = path, sheet = "veh_age") +vmt_age <- readxl::read_xlsx(path = path, sheet = "vmt_age") +fuel <- readxl::read_xlsx(path = path, sheet = "fuel") +fuel_type <- readxl::read_xlsx(path = path, sheet = "fuel_type") +met <- readxl::read_xlsx(path = path, sheet = "met") +year <- 2019 +month <- 8 +theme <- "black" # dark clean ink +delete_directories <- T +source("config/config.R", encoding = "UTF-8") + +# 1) Network #### +# source("network/pre_net_simu_BMC_qgisv2.R") +# source("network/pre_net_GTFS.R") +language <- "english" # english spanish portuguese +path <- "config/inventory_MD.xlsx" +metadata <- readxl::read_xlsx(path = path, sheet = "metadata") +net <- sf::st_read("network/net_baltimore_city.gpkg") +crs <- 6487 +tit <- "Travel demand model outptus from Baltimore 2019 (veh/h)" +names(net) +categories <- unique(metadata$vein_name) # in network/net.gpkg +source("scripts/net.R", encoding = "UTF-8") +rm(list = ls()) +gc() + +# 2) Traffic #### +language <- "english" # portuguese spanish +metadata <- readRDS("config/metadata.rds") +net <- readRDS("network/net.rds") +tfs <- readRDS("config/tfs.rds") +veh <- readRDS("config/fleet_age.rds") +categories <- unique(metadata$vein_name) # in net +verbose <- FALSE +year <- 2019 +theme <- "ink" # dark clean ink +speed <- TRUE +ps <- "CONGSPD" +ffs <- "FFSPEED" +capacity <- "CAPE" +dist <- "miles" +k_D <- k_CNG <- k_G <- k_E85 <- 1 +source("scripts/traffic.R", encoding = "UTF-8") +rm(list = ls()) +gc() + +# 3) Estimation #### +language <- "english" # portuguese english spanish +metadata <- readRDS("config/metadata.rds") +net <- readRDS("network/net.rds") +tfs <- readRDS("config/tfs.rds") +veh <- readRDS("config/fleet_age.rds") +speed <- readRDS("network/speed_net.rds") +speed_bin <- readRDS("network/speed_bin.rds") +fuel <- readRDS("config/fuel.rds") +fuel_type <- readRDS("config/fuel_type.rds") +lkm <- net$lkm +units(net$lkm)$numerator # miles +verbose <- TRUE +year <- 2019 + +# SQL MOVES PARAMETERS +host <- "localhost" +moves_db <- "baltimore_2017_aug_emissionfactor" +user <- "moves" +password <- "moves" +port <- 3307 + +# 3.1 Fuel eval #### +source("scripts/fuel_eval.R", encoding = "UTF-8") +rm(list = ls()) +gc() + +# 3.2 Emissions per distance #### +language <- "english" # portuguese spanish +metadata <- readRDS("config/metadata.rds") +net <- readRDS("network/net.rds") +tfs <- readRDS("config/tfs.rds") +veh <- readRDS("config/fleet_age.rds") +speed <- readRDS("network/speed_net.rds") +speed_bin <- readRDS("network/speed_bin.rds") +fuel_type <- readRDS("config/fuel_type.rds") +lkm <- net$lkm +verbose <- FALSE + +# SQL MOVES PARAMETERS +host <- "localhost" +moves_db <- "baltimore_2017_aug_emissionfactor" # _mi_aft _mo +user <- "moves" +password <- "moves" +port <- 3307 +source("scripts/emissions_per_distance.R", encoding = "UTF-8") +rm(list = ls()) +gc() + + +language <- "english" # english spanish +metadata <- readRDS("config/metadata.rds") +net <- readRDS("network/net.rds") +tfs <- readRDS("config/tfs.rds") +veh <- readRDS("config/fleet_age.rds") +speed <- readRDS("network/speed_net.rds") +speed_bin <- readRDS("network/speed_bin.rds") +fuel_type <- readRDS("config/fuel_type.rds") +lkm <- net$lkm +verbose <- FALSE + +# SQL MOVES PARAMETERS +host <- "localhost" +moves_db <- "baltimore_2017_aug_emissionfactor" # _mi_aft _mo +user <- "moves" +password <- "moves" +port <- 3307 + +source("scripts/emissions_per_distance_NMHC.R", encoding = "UTF-8") +rm(list = ls()) +gc() + +# 3.2 Emissions per vehicle #### +language <- "english" # english spanish +metadata <- readRDS("config/metadata.rds") +net <- readRDS("network/net.rds") +tfs <- readRDS("config/tfs.rds") +veh <- readRDS("config/fleet_age.rds") +speed <- readRDS("network/speed_net.rds") +speed_bin <- readRDS("network/speed_bin.rds") +fuel_type <- readRDS("config/fuel_type.rds") +lkm <- net$lkm +verbose <- TRUE + +# SQL MOVES PARAMETERS +host <- "localhost" +moves_db <- "baltimore_2017_aug_emissionfactor" # _mi_aft _mo +user <- "moves" +password <- "moves" +port <- 3307 + +source("scripts/emissions_per_vehicle.R", encoding = "UTF-8") +rm(list = ls()) +gc() + + +language <- "english" # english spanish +metadata <- readRDS("config/metadata.rds") +net <- readRDS("network/net.rds") +tfs <- readRDS("config/tfs.rds") +veh <- readRDS("config/fleet_age.rds") +speed <- readRDS("network/speed_net.rds") +speed_bin <- readRDS("network/speed_bin.rds") +fuel_type <- readRDS("config/fuel_type.rds") +lkm <- net$lkm +verbose <- TRUE + +# SQL MOVES PARAMETERS +host <- "localhost" +moves_db <- "baltimore_2017_aug_emissionfactor" # _mi_aft _mo +user <- "moves" +password <- "moves" +port <- 3307 + +source("scripts/emissions_per_vehicle_NMHC.R", encoding = "UTF-8") +rm(list = ls()) +gc() + +# 4) Post-estimation #### +data(decoder) +metadata <- readRDS("config/metadata.rds") +language <- "english" # portuguese spanish +net <- readRDS("network/net.rds") +tfs <- readRDS("config/tfs.rds") +grid_width <- 200 # mts +crs <- 6487 +source("scripts/post.R", encoding = "UTF-8") +rm(list = ls()) +gc() + +# plots #### +language <- "english" # portuguese spanish +metadata <- readRDS("config/metadata.rds") +tfs <- readRDS("config/tfs.rds") +hours <- 8 +bg <- "white" +pal <- "colo_lightningmccarl_into_the_night" # find more in ?cptcity::find_cpt +tit <- "Baltimore City" +source("scripts/plots.R") +rm(list = ls()) +gc() + +# MECH #### +# mechs <- c("CB4", "CB05", "CB05opt2", "S99", "S7","CS7", "S7T", "S11", "S11D", +# "S16C","S18B","RADM2", "RACM2","MOZT1" "CBMZ") +mechs <- c("CB05", "CB05opt2") +for (z in seq_along(mechs)) { + print(mechs[z]) + language <- "english" # portuguese spanish + net <- readRDS("network/net.rds") + tfs <- readRDS("config/tfs.rds") + mech <- mechs[z] + source("scripts/mech.R", encoding = "UTF-8") +} +rm(list = ls()) +gc() \ No newline at end of file diff --git a/projects/moves/main.Rproj b/projects/moves/main.Rproj new file mode 100755 index 00000000..8e3c2ebc --- /dev/null +++ b/projects/moves/main.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX diff --git a/projects/moves/network/net_baltimore_city.gpkg b/projects/moves/network/net_baltimore_city.gpkg new file mode 100755 index 00000000..5544cdf7 Binary files /dev/null and b/projects/moves/network/net_baltimore_city.gpkg differ diff --git a/projects/moves/network/net_no_bus.gpkg b/projects/moves/network/net_no_bus.gpkg new file mode 100755 index 00000000..5deccd06 Binary files /dev/null and b/projects/moves/network/net_no_bus.gpkg differ diff --git a/projects/moves/network/net_qgis.gpkg b/projects/moves/network/net_qgis.gpkg new file mode 100755 index 00000000..0d3128a3 Binary files /dev/null and b/projects/moves/network/net_qgis.gpkg differ diff --git a/projects/moves/network/pre_net.R b/projects/moves/network/pre_net.R new file mode 100755 index 00000000..500bb680 --- /dev/null +++ b/projects/moves/network/pre_net.R @@ -0,0 +1,222 @@ +library(data.table) +library(sf) +library(cptcity) +library(ggplot2) +net <- sf::st_read("../maryland_email.gpkg") +net <- st_transform(net, 6487) +names(net) +net$lanes <- as.numeric(as.character(net$SINGLE_UNIT_AADT)) + +metadata <- readxl::read_excel("../excel_inventory/inventory_MD.xlsx", "metadata") +setDT(metadata) +pop <- readxl::read_excel("../excel_inventory/inventory_MD.xlsx", "pop_moves") +setDT(pop) + +# Motorcycle +# Passenger Car +# Passenger Truck +# Light Commercial Truck +# Intercity Bus +# Transit Bus +# School Bus +# Refuse Truck +# Single Unit Short-haul Truck +# Single Unit Long-haul Truck +# Motor Home +# Combination Short-haul Truck +# Combination Long-haul Truck +# LDV <- PC + PT +class(net$MOTORCYCLE_AADT) +net$MC <- net$MOTORCYCLE_AADT + +LDV <- pop[sourceTypeID %in% c(21, 32), sum(sourceTypePopulation)] +per_PC <- pop[sourceTypeID %in% c(21), sum(sourceTypePopulation)]/LDV +per_PT <- pop[sourceTypeID %in% c(31), sum(sourceTypePopulation)]/LDV +class(net$CAR_AADT) +net$PC <- net$CAR_AADT*per_PC +net$PT <- net$CAR_AADT*per_PT + +class(net$LIGHT_TRUCK_AADT) +net$LCT <- net$LIGHT_TRUCK_AADT + +BUS <- pop[sourceTypeID %in% 41:43, sum(sourceTypePopulation)] +per_IBUS <- pop[sourceTypeID %in% c(41), sum(sourceTypePopulation)]/BUS +per_TBUS <- pop[sourceTypeID %in% c(42), sum(sourceTypePopulation)]/BUS +per_SBUS <- pop[sourceTypeID %in% c(43), sum(sourceTypePopulation)]/BUS +class(net$BUS_AADT) +net$BUS_INTERCITY <- net$BUS_AADT*per_IBUS +net$BUS_TRANSIT <- net$BUS_AADT*per_TBUS +net$BUS_SCHOOL <- net$BUS_AADT*per_SBUS + + +TRUCKS_SU <- pop[sourceTypeID %in% 51:54, sum(sourceTypePopulation)] +per_T_REFUSE <- pop[sourceTypeID %in% c(51), sum(sourceTypePopulation)]/TRUCKS_SU +per_T_SU_SW <- pop[sourceTypeID %in% c(52), sum(sourceTypePopulation)]/TRUCKS_SU +per_T_SU_LW <- pop[sourceTypeID %in% c(53), sum(sourceTypePopulation)]/TRUCKS_SU +per_T_MH <- pop[sourceTypeID %in% c(54), sum(sourceTypePopulation)]/TRUCKS_SU +class(net$SINGLE_UNIT_AADT) +net$TRUCKS_REFUSE <- net$SINGLE_UNIT_AADT*per_T_REFUSE +net$TRUCKS_SU_SW <- net$SINGLE_UNIT_AADT*per_T_SU_SW +net$TRUCKS_SU_LW <- net$SINGLE_UNIT_AADT*per_T_SU_LW +net$TRUCKS_MH <- net$SINGLE_UNIT_AADT*per_T_MH + +class(net$COMBINATION_UNIT_AADT) +net$TRUCKS_CU_SH <- net$COMBINATION_UNIT_AADT/2 +net$TRUCKS_CU_LH <- net$COMBINATION_UNIT_AADT/2 + +class(net$AADT_2018) +net$AADT_2018 <- net$AADT_2018 + +net$id <- 1:nrow(net) + +ids <- net[is.na(net$PC), ]$id + +pbapply::pblapply(seq_along(ids), function(i) { + x <- st_buffer(net[net$id == ids[i], ], 1000) + xt <- st_intersection(net, x) + xt <- st_set_geometry(xt, NULL) + TOT <- sum(xt[, c("MC", + "PC", "PT", + "LCT", + "BUS_INTERCITY", "BUS_TRANSIT", "BUS_SCHOOL", + "TRUCKS_REFUSE", "TRUCKS_SU_SW", "TRUCKS_SU_LW", "TRUCKS_MH", + "TRUCKS_CU_SH", "TRUCKS_CU_LH")], + na.rm = T) + df <- data.frame(MC = sum(xt$MC, na.rm = T)/TOT, + PC = sum(xt$PC, na.rm = T)/TOT, + PT = sum(xt$PT, na.rm = T)/TOT, + LCT = sum(xt$LCT, na.rm = T)/TOT, + BUS_INTERCITY = sum(xt$BUS_INTERCITY, na.rm = T)/TOT, + BUS_TRANSIT = sum(xt$BUS_TRANSIT, na.rm = T)/TOT, + BUS_SCHOOL = sum(xt$BUS_SCHOOL, na.rm = T)/TOT, + TRUCKS_REFUSE = sum(xt$TRUCKS_REFUSE, na.rm = T)/TOT, + TRUCKS_SU_SW = sum(xt$TRUCKS_SU_SW, na.rm = T)/TOT, + TRUCKS_SU_LW = sum(xt$TRUCKS_SU_LW, na.rm = T)/TOT, + TRUCKS_MH = sum(xt$TRUCKS_MH, na.rm = T)/TOT, + TRUCKS_CU_SH = sum(xt$TRUCKS_CU_SH, na.rm = T)/TOT, + TRUCKS_CU_LH = sum(xt$TRUCKS_CU_LH, na.rm = T)/TOT + ) + df$id <- ids[i] + df +}) -> dx + +dxb <- dx +dx <- rbindlist(dx) +# saveRDS(dx, "dx.rds") + +names(dx)[1:(ncol(dx) - 1)] <- paste0("BUF_", names(dx)[1:(ncol(dx) - 1)]) + +net <- merge(net, dx, by = "id", all.x = T) + +net +geo <- net$geometry + +setDT(net) + +net[, + MC_FINAL := fifelse(is.na(MC), + BUF_MC*AADT_2018, + MC)] +net[, + PC_FINAL := fifelse(is.na(PC), + BUF_PC*AADT_2018, + PC)] + +net[, + PT_FINAL := fifelse(is.na(PT), + BUF_PT*AADT_2018, + PT)] + +net[, + LCT_FINAL := fifelse(is.na(LCT), + BUF_LCT*AADT_2018, + LCT)] + +net[, + BUS_INTERCITY_FINAL := fifelse(is.na(BUS_INTERCITY), + BUF_BUS_INTERCITY*AADT_2018, + BUS_INTERCITY)] + +net[, + BUS_TRANSIT_FINAL := fifelse(is.na(BUS_TRANSIT), + BUF_BUS_TRANSIT*AADT_2018, + BUS_TRANSIT)] + + +net[, + BUS_SCHOOL_FINAL := fifelse(is.na(BUS_SCHOOL), + BUF_BUS_SCHOOL*AADT_2018, + BUS_SCHOOL)] + + + +net[, + TRUCKS_REFUSE_FINAL := fifelse(is.na(TRUCKS_REFUSE), + BUF_TRUCKS_REFUSE*AADT_2018, + TRUCKS_REFUSE)] + + + + +net[, + TRUCKS_SU_SW_FINAL := fifelse(is.na(TRUCKS_SU_SW), + BUF_TRUCKS_SU_SW*AADT_2018, + TRUCKS_SU_SW)] + + +net[, + TRUCKS_SU_LW_FINAL := fifelse(is.na(TRUCKS_SU_LW), + BUF_TRUCKS_SU_LW*AADT_2018, + TRUCKS_SU_LW)] + + +net[, + TRUCKS_MH_FINAL := fifelse(is.na(TRUCKS_MH), + BUF_TRUCKS_MH*AADT_2018, + TRUCKS_MH)] + + +net[, + TRUCKS_CU_SH_FINAL := fifelse(is.na(TRUCKS_CU_SH), + BUF_TRUCKS_CU_SH*AADT_2018, + TRUCKS_CU_SH)] + +net[, + TRUCKS_CU_LH_FINAL := fifelse(is.na(TRUCKS_CU_LH), + BUF_TRUCKS_CU_LH*AADT_2018, + TRUCKS_CU_LH)] + + +# Motorcycle +# Passenger Car +# Passenger Truck +# Light Commercial Truck +# Intercity Bus +# Transit Bus +# School Bus +# Refuse Truck +# Single Unit Short-haul Truck +# Single Unit Long-haul Truck +# Motor Home +# Combination Short-haul Truck +# Combination Long-haul Truck +# LDV <- PC + PT + +net_final <- net[, c(paste0(c("MC", + "PC", "PT", + "LCT", + "BUS_INTERCITY", "BUS_TRANSIT", "BUS_SCHOOL", + "TRUCKS_REFUSE","TRUCKS_SU_SW", "TRUCKS_SU_LW", "TRUCKS_MH", + "TRUCKS_CU_SH", "TRUCKS_CU_LH"), + "_FINAL"), "NUM_LANES", "AAWDT_2018")] + names(net_final) <- gsub("_FINAL", "", names(net_final)) +net_final <- st_sf(net_final, + geometry = geo) + +file.remove("network/net.gpkg") +st_write(net_final, "network/net.gpkg") +# plot(net_final["PC"], axes = T) +# edit inventory.xlsx +rm(list = ls()) +gc() +´ \ No newline at end of file diff --git a/projects/moves/network/pre_net_GTFS.R b/projects/moves/network/pre_net_GTFS.R new file mode 100755 index 00000000..96435955 --- /dev/null +++ b/projects/moves/network/pre_net_GTFS.R @@ -0,0 +1,109 @@ +library(sf) +library(gtfs2gps) +library(magrittr) +library(data.table) + +x <- st_read("network/net_no_bus.gpkg") +names(x) +x$id <- 1:nrow(x) +crs_original <- st_crs(x) + +# GTFS BUS LOCAL #### +p2 <- "../GTFS/mdotmta_gtfs_localbus.zip" +curi <- gtfs2gps::read_gtfs(p2) %>% + gtfs2gps::filter_day_period(period_start = "08:00", period_end = "09:00") %>% + gtfs2gps::gtfs2gps(parallel = F) %>% + gtfs2gps::gps_as_sfpoints() +curi$bus <- 1 + +curi <- curi[as.numeric(curi$speed) < 100, ] +bmc <- st_read("../../shapefiles/limits/CSA_NSA_Tracts.shp") +bmc <- st_union(bmc) +bmc <- sfheaders::sf_remove_holes(bmc) +plot(bmc) + +# intersectar con red y luego agregar por id +curi <- st_transform(curi, st_crs(bmc)) +x <- st_transform(x, st_crs(bmc)) + +dim(curi) +curi <- st_crop(curi, bmc) +dim(curi) + +curib <- st_buffer(curi, 20) +curib$id <- 1:nrow(curib) + +curit <- st_intersection(curib["bus"], x["id"]) + +dim(curit) +setDT(curit) +xx <- curit[, sum(bus), by = id] +names(xx)[2] <- "bus" +x <- merge(x, xx, by = "id", all.x = T) +plot(x["bus"], axes = T, pal = cptcity::cpt(colorRampPalette = T, rev = T)) + +# comparison pop moves #### +pop <- readxl::read_excel("config/inventory_MD.xlsx", "pop_moves") +setDT(pop) +fraction_bus <- pop[sourceTypeID %in% 42:43, sum(sourceTypePopulation)]/sum(pop$sourceTypePopulation) + +y <- fraction_bus * ( pop[!sourceTypeID %in% 42:43, sum(sourceTypePopulation)]) / sum(x$bus, na.rm = T) + +sum(x$bus * y, na.rm = T) / (pop[!sourceTypeID %in% 42:43, sum(sourceTypePopulation)]) + +x$bus <- x$bus * y +x$BUS_TRANSIT <- x$bus * pop[!sourceTypeID %in% 42, sum(sourceTypePopulation)] / pop[!sourceTypeID %in% 42:43, sum(sourceTypePopulation)] +x$BUS_SCHOOL <- x$bus * pop[!sourceTypeID %in% 43, sum(sourceTypePopulation)] / pop[!sourceTypeID %in% 42:43, sum(sourceTypePopulation)] +net <- x + + +# GTFS BUS INTERCITY #### +p2 <- "../GTFS/mdotmta_gtfs_commuterbus.zip" +curi <- gtfs2gps::read_gtfs(p2) %>% + gtfs2gps::filter_day_period(period_start = "08:00", period_end = "09:00") %>% + gtfs2gps::gtfs2gps(parallel = F) %>% + gtfs2gps::gps_as_sfpoints() +curi$bus <- 1 + +curi <- curi[as.numeric(curi$speed) < 100, ] +bmc <- st_read("../../shapefiles/limits/CSA_NSA_Tracts.shp") +bmc <- st_union(bmc) +bmc <- sfheaders::sf_remove_holes(bmc) +plot(bmc) + +# intersectar con red y luego agregar por id +curi <- st_transform(curi, st_crs(bmc)) +x <- st_transform(x, st_crs(bmc)) + +dim(curi) +curi <- st_crop(curi, bmc) +dim(curi) + +curib <- st_buffer(curi, 20) +curib$id <- 1:nrow(curib) + +curit <- st_intersection(curib["bus"], x["id"]) + +dim(curit) +setDT(curit) +xx <- curit[, sum(bus), by = id] +names(xx)[2] <- "bus2" +x <- merge(x, xx, by = "id", all.x = T) +plot(x["bus2"], axes = T, pal = cptcity::cpt(colorRampPalette = T, rev = T)) + +# comparison pop moves #### +pop <- readxl::read_excel("config/inventory_MD.xlsx", "pop_moves") +setDT(pop) +fraction_bus <- pop[sourceTypeID %in% 41, sum(sourceTypePopulation)]/sum(pop$sourceTypePopulation) + +y <- fraction_bus * ( pop[!sourceTypeID %in% 41, sum(sourceTypePopulation)]) / sum(x$bus2, na.rm = T) + +sum(x$bus2 * y, na.rm = T) / (pop[!sourceTypeID %in% 41, sum(sourceTypePopulation)]) + +x$bus2 <- x$bus2 * y +x$BUS_INTERCITY <- x$bus2 + +x +st_write(x, "network/net_baltimore_city.gpkg", delete_layer = T) +rm(list = ls()) +gc() diff --git a/projects/moves/network/pre_net_simu_BMC.R b/projects/moves/network/pre_net_simu_BMC.R new file mode 100755 index 00000000..06d45c4d --- /dev/null +++ b/projects/moves/network/pre_net_simu_BMC.R @@ -0,0 +1,205 @@ +library(data.table) +library(sf) +library(cptcity) +library(ggplot2) +# netemail <- sf::st_read("../vein_baltimore/network/net.gpkg") +# netemail <- st_transform(netemail, 4326) +net <- sf::st_read("../BMC_SIMU/HWY/HWY_LDNET_AM2_4.shp") +net <- st_crop(net, md) +names(net) +plot(net$geometry) +# writexl::write_xlsx(net, "config/net.xlsx") +# net <- st_transform(net, 6487) + +# plot(netemail["MC"], axes = T) +# dim(net) +# find_cpt("ocal") +nr <- classInt::classIntervals(var = net$F1[net$F1 > 0], n = 100, style = "sd") + + +# attached zip file contains three shape files - am loaded highway +# network 1) 6:30-7:30, 2) 7:30-8:30, and 3) 8:30 - 9:30. +# The horizon year is 2022, which is the closest we have to 2019. +# +# Link variables (VOL1, VOL2...VOL20) are defined below. +# Trip tables are divided into 5 value of time (1 lowest...5 highest). +# +# Will follow up with additional documentation on link variables. + +# PHASE=ILOOP +# ; SOV (+ SOV Air & Drive to Transit in VOT3) +# PATH = LW.COST_SOV1, PENI = 1, VOL[1]=MI.1.1,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_SOV2, PENI = 1, VOL[2]=MI.1.2,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_SOV3, PENI = 1, VOL[3]=MI.1.3+MI.2.6+MI.3.@PER@,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_SOV4, PENI = 1, VOL[4]=MI.1.4,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_SOV5, PENI = 1, VOL[5]=MI.1.5,EXCLUDEGRP=1,4,7 +# +# ; HOV2 (+ HOV Air in VOT3) +# PATH = LW.COST_HOV21, PENI = 1, VOL[6]=MI.1.6,EXCLUDEGRP=4,7 +# PATH = LW.COST_HOV22, PENI = 1, VOL[7]=MI.1.7,EXCLUDEGRP=4,7 +# PATH = LW.COST_HOV23, PENI = 1, VOL[8]=MI.1.8+MI.2.7,EXCLUDEGRP=4,7 +# PATH = LW.COST_HOV24, PENI = 1, VOL[9]=MI.1.9,EXCLUDEGRP=4,7 +# PATH = LW.COST_HOV25, PENI = 1, VOL[10]=MI.1.10,EXCLUDEGRP=4,7 +# +# ; HOV3 +# PATH = LW.COST_HOV31, PENI = 1, VOL[11]=MI.1.11,EXCLUDEGRP=4 +# PATH = LW.COST_HOV32, PENI = 1, VOL[12]=MI.1.12,EXCLUDEGRP=4 +# PATH = LW.COST_HOV33, PENI = 1, VOL[13]=MI.1.13,EXCLUDEGRP=4 +# PATH = LW.COST_HOV34, PENI = 1, VOL[14]=MI.1.14,EXCLUDEGRP=4 +# PATH = LW.COST_HOV35, PENI = 1, VOL[15]=MI.1.15,EXCLUDEGRP=4 +# +# ; TRUCKS +# PATH = LW.COST_CV, PENI = 1, VOL[16]=MI.2.1,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_MT, PENI = 1, VOL[17]=MI.2.2,EXCLUDEGRP=1,3,4,7 +# PATH = LW.COST_HT, PENI = 1, VOL[18]=MI.2.3,EXCLUDEGRP=1,3,4,7 +# +# ; XXPC +# PATH = LW.COST_SOVXX, PENI = 1, VOL[19]=MI.2.4 EXCLUDEGRP=1,2,4,7 +# PATH = LW.COST_HOVXX, PENI = 1, VOL[20]=MI.2.5 EXCLUDEGRP=2,4 +# ENDPROCESS + +# F1-F13 is 24 hour counts...but i am not sure on year - Brian? +# +# V_1 is the sum of volumes 1-20 +# V1_1 through V1_20 matches the list from the previous emails. VOL 1 SOV VOT1 +# +# V1T_1 is the total volume AB V1_1 + BA V1_1 +# +# Capacity is CAPCLASS * 50 or you can use the value in CAPE. +# + +# net$PC_TOTAL <- net$V19T_1 + net$V20T_1 +# names(net)[120] <- "XXPC = V19T_1 + V20T_1" +# net$PC_AB <- net$V19_1 + net$V20_1 +# names(net)[121] <- "XXPC = V19_1 + V20_1" +# +# plot(net[, c("XXPC = V19T_1 + V20T_1", +# "XXPC = V19_1 + V20_1")], +# key.pos = 4, +# pal = cpt(rev = T, colorRampPalette = T), +# axes = T, +# breaks = "sd") +# HOV: High occupancy vehicle https://en.wikipedia.org/wiki/High-occupancy_vehicle_lane +# SOV: Single occupancy vehicle http://www.nwarpc.org/pdf/Transportation/Sustainable_Mobility/Maximizing_Exchange_Minimizing_Travel.pdf + +# F1 -> F13 +# F1 FHWA Vehicle Class 1 motorcycle +# F2 FHWA Vehicle Class 2 passenger car +# F3 FHWA Vehicle Class 3 2 axle, 4 tire single units +# F4 FHWA Vehicle Class 4 buses +# F5 FHWA Vehicle Class 5 2 axle, 6 tire single units +# F6 FHWA Vehicle Class 6 3 axle single units +# F7 FHWA Vehicle Class 7 4 or more axle single units +# F8 FHWA Vehicle Class 8 4 or less axle single trailers +# F9 FHWA Vehicle Class 9 5 axle single trailers +# F10 FHWA Vehicle Class 10 6 plus axle single trailers +# F11 FHWA Vehicle Class 11 5 or less axle multi-trailers +# F12 FHWA Vehicle Class 12 6 axle multi-trailers +# F13 FHWA Vehicle Class 13 7 or more axle multi-trailers + +# goal #### +# The goal is to have volumes by MOVES categories +# Then, need toy data to proceed +# We can only be sure for PC = V19 + v20 and trucks v16 + v17 + v18 +# HOV3 seems to be mostly buses, we can assume that for now v11 + v12 + v13 + v14 + v15 + +# the other volumes can have a total distribution +sapply(st_set_geometry(net[, c("V19T_1", "V20T_1", "F2", "F1")], NULL), sum) + +metadata <- readxl::read_excel("../excel_inventory/inventory_MD.xlsx", "metadata") +setDT(metadata) +pop <- readxl::read_excel("../excel_inventory/inventory_MD.xlsx", "pop_moves") +setDT(pop) + +# Motorcycle = F1 +# Passenger Car = F2 +# Passenger Truck = F3 +# Light Commercial Truck = F3 +# Intercity Bus = F4 +# Transit Bus = F4 +# School Bus = F4 +# Refuse Truck = F5 + F6 + F8 +# Single Unit Short-haul Truck = F5 + F6 + F8 +# Single Unit Long-haul Truck = F7 + F9 + F10 +# Motor Home = F5 + F6 + F8 +# Combination Short-haul Truck = F11 +# Combination Long-haul Truck = F13 + F12 + +# assuming vols 1 to 10 covers all vehicles +# assumes trucks only on SOV, which makes sense +net$ALL <- net$V1T_1 + net$V2T_1 + net$V3T_1 + net$V4T_1 + net$V5T_1 +efes <- st_set_geometry(net[, paste0("F", 1:13)], NULL) +setDT(efes) +efes <- efes[rowSums(efes) > 0] + +col_efes <- colSums(efes)/sum(colSums(efes)) +df <- as.data.table(matrix(col_efes, nrow = 1)) +names(df) <- paste0("F", 1:13) + +df + +unique(pop$sourceTypeID) +MC_PC <- pop[sourceTypeID %in% c(11, 21), sum(sourceTypePopulation)] +per_MC <- pop[sourceTypeID %in% c(11), sum(sourceTypePopulation)]/MC_PC +per_PC <- pop[sourceTypeID %in% c(21), sum(sourceTypePopulation)]/MC_PC +# assumes V19T + V6T ...V10T = MC + PC +# assumes net$ALL*dfF1 = MC +# assumes net$ALL*dfF2 = PC + +net$MC <- (net$V19T_1 + net$V6T_1 +net$V7T_1 +net$V8T_1 +net$V9T_1 +net$V10T_1)*per_MC + net$ALL*df$F1 +net$PC <- (net$V19T_1 + net$V6T_1 +net$V7T_1 +net$V8T_1 +net$V9T_1 +net$V10T_1)*per_PC + net$ALL*df$F2 + + +PT_LCT <- pop[sourceTypeID %in% c(31, 32), sum(sourceTypePopulation)] +per_PT <- pop[sourceTypeID %in% c(31), sum(sourceTypePopulation)]/PT_LCT +per_LCT <- pop[sourceTypeID %in% c(31), sum(sourceTypePopulation)]/PT_LCT + +# assuming V20T are PT +net$PT <- (net$V20T_1 + net$ALL*df$F3)*per_PT +net$LCT <- (net$V20T_1 + net$ALL*df$F3)*per_LCT + +BUS <- pop[sourceTypeID %in% 41:43, sum(sourceTypePopulation)] +per_IBUS <- pop[sourceTypeID %in% c(41), sum(sourceTypePopulation)]/BUS +per_TBUS <- pop[sourceTypeID %in% c(42), sum(sourceTypePopulation)]/BUS +per_SBUS <- pop[sourceTypeID %in% c(43), sum(sourceTypePopulation)]/BUS +net$BUS_INTERCITY <- (net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1)*per_IBUS + net$ALL*df$F4*per_IBUS +per_IBUS*(net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1 + net$ALL*df$F4) + +net$BUS_TRANSIT <- (net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1)*per_TBUS+ net$ALL*df$F4*per_TBUS + +net$BUS_SCHOOL <- (net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1)*per_SBUS+ net$ALL*df$F4*per_SBUS + + +TRUCKS <- pop[sourceTypeID %in% c(51:54, 61, 62), sum(sourceTypePopulation)] +per_T_REFUSE <- pop[sourceTypeID %in% c(51), sum(sourceTypePopulation)]/TRUCKS +per_T_SU_SH <- pop[sourceTypeID %in% c(52), sum(sourceTypePopulation)]/TRUCKS +per_T_SU_LH <- pop[sourceTypeID %in% c(53), sum(sourceTypePopulation)]/TRUCKS +per_T_MH <- pop[sourceTypeID %in% c(54), sum(sourceTypePopulation)]/TRUCKS +per_T_CU_SH <- pop[sourceTypeID %in% c(61), sum(sourceTypePopulation)]/TRUCKS +per_T_CU_LH <- pop[sourceTypeID %in% c(62), sum(sourceTypePopulation)]/TRUCKS +# assumes trucks + sum V15..V18 +# assumes net$ALL (SOV) covers trucks F7...F13 +net$TRUCKS_REFUSE <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_REFUSE + +net$TRUCKS_SU_SH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_SU_SH + +net$TRUCKS_SU_LH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_SU_LH + +net$TRUCKS_MH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_MH + + +net$TRUCKS_CU_SH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_CU_SH +net$TRUCKS_CU_LH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_CU_LH + +file.remove("network/net.gpkg") +st_write(net, "network/net.gpkg") +# plot(net["PC"], axes = T, pal = cpt(colorRampPalette = T, rev = T), key.pos = 1) +# edit inventory.xlsx +rm(list = ls()) +gc() diff --git a/projects/moves/network/pre_net_simu_BMC_qgis.R b/projects/moves/network/pre_net_simu_BMC_qgis.R new file mode 100755 index 00000000..dc315a68 --- /dev/null +++ b/projects/moves/network/pre_net_simu_BMC_qgis.R @@ -0,0 +1,203 @@ +library(data.table) +library(sf) +library(cptcity) +library(ggplot2) +# netemail <- sf::st_read("../vein_baltimore/network/net.gpkg") +# netemail <- st_transform(netemail, 4326) +net <- sf::st_read("network/net_qgis.gpkg") +plot(net$geom) +# writexl::write_xlsx(net, "config/net.xlsx") +# net <- st_transform(net, 6487) + +# plot(netemail["MC"], axes = T) +# dim(net) +# find_cpt("ocal") +nr <- classInt::classIntervals(var = net$F1[net$F1 > 0], n = 100, style = "sd") + + +# attached zip file contains three shape files - am loaded highway +# network 1) 6:30-7:30, 2) 7:30-8:30, and 3) 8:30 - 9:30. +# The horizon year is 2022, which is the closest we have to 2019. +# +# Link variables (VOL1, VOL2...VOL20) are defined below. +# Trip tables are divided into 5 value of time (1 lowest...5 highest). +# +# Will follow up with additional documentation on link variables. + +# PHASE=ILOOP +# ; SOV (+ SOV Air & Drive to Transit in VOT3) +# PATH = LW.COST_SOV1, PENI = 1, VOL[1]=MI.1.1,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_SOV2, PENI = 1, VOL[2]=MI.1.2,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_SOV3, PENI = 1, VOL[3]=MI.1.3+MI.2.6+MI.3.@PER@,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_SOV4, PENI = 1, VOL[4]=MI.1.4,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_SOV5, PENI = 1, VOL[5]=MI.1.5,EXCLUDEGRP=1,4,7 +# +# ; HOV2 (+ HOV Air in VOT3) +# PATH = LW.COST_HOV21, PENI = 1, VOL[6]=MI.1.6,EXCLUDEGRP=4,7 +# PATH = LW.COST_HOV22, PENI = 1, VOL[7]=MI.1.7,EXCLUDEGRP=4,7 +# PATH = LW.COST_HOV23, PENI = 1, VOL[8]=MI.1.8+MI.2.7,EXCLUDEGRP=4,7 +# PATH = LW.COST_HOV24, PENI = 1, VOL[9]=MI.1.9,EXCLUDEGRP=4,7 +# PATH = LW.COST_HOV25, PENI = 1, VOL[10]=MI.1.10,EXCLUDEGRP=4,7 +# +# ; HOV3 +# PATH = LW.COST_HOV31, PENI = 1, VOL[11]=MI.1.11,EXCLUDEGRP=4 +# PATH = LW.COST_HOV32, PENI = 1, VOL[12]=MI.1.12,EXCLUDEGRP=4 +# PATH = LW.COST_HOV33, PENI = 1, VOL[13]=MI.1.13,EXCLUDEGRP=4 +# PATH = LW.COST_HOV34, PENI = 1, VOL[14]=MI.1.14,EXCLUDEGRP=4 +# PATH = LW.COST_HOV35, PENI = 1, VOL[15]=MI.1.15,EXCLUDEGRP=4 +# +# ; TRUCKS +# PATH = LW.COST_CV, PENI = 1, VOL[16]=MI.2.1,EXCLUDEGRP=1,4,7 +# PATH = LW.COST_MT, PENI = 1, VOL[17]=MI.2.2,EXCLUDEGRP=1,3,4,7 +# PATH = LW.COST_HT, PENI = 1, VOL[18]=MI.2.3,EXCLUDEGRP=1,3,4,7 +# +# ; XXPC +# PATH = LW.COST_SOVXX, PENI = 1, VOL[19]=MI.2.4 EXCLUDEGRP=1,2,4,7 +# PATH = LW.COST_HOVXX, PENI = 1, VOL[20]=MI.2.5 EXCLUDEGRP=2,4 +# ENDPROCESS + +# F1-F13 is 24 hour counts...but i am not sure on year - Brian? +# +# V_1 is the sum of volumes 1-20 +# V1_1 through V1_20 matches the list from the previous emails. VOL 1 SOV VOT1 +# +# V1T_1 is the total volume AB V1_1 + BA V1_1 +# +# Capacity is CAPCLASS * 50 or you can use the value in CAPE. +# + +# net$PC_TOTAL <- net$V19T_1 + net$V20T_1 +# names(net)[120] <- "XXPC = V19T_1 + V20T_1" +# net$PC_AB <- net$V19_1 + net$V20_1 +# names(net)[121] <- "XXPC = V19_1 + V20_1" +# +# plot(net[, c("XXPC = V19T_1 + V20T_1", +# "XXPC = V19_1 + V20_1")], +# key.pos = 4, +# pal = cpt(rev = T, colorRampPalette = T), +# axes = T, +# breaks = "sd") +# HOV: High occupancy vehicle https://en.wikipedia.org/wiki/High-occupancy_vehicle_lane +# SOV: Single occupancy vehicle http://www.nwarpc.org/pdf/Transportation/Sustainable_Mobility/Maximizing_Exchange_Minimizing_Travel.pdf + +# F1 -> F13 +# F1 FHWA Vehicle Class 1 motorcycle +# F2 FHWA Vehicle Class 2 passenger car +# F3 FHWA Vehicle Class 3 2 axle, 4 tire single units +# F4 FHWA Vehicle Class 4 buses +# F5 FHWA Vehicle Class 5 2 axle, 6 tire single units +# F6 FHWA Vehicle Class 6 3 axle single units +# F7 FHWA Vehicle Class 7 4 or more axle single units +# F8 FHWA Vehicle Class 8 4 or less axle single trailers +# F9 FHWA Vehicle Class 9 5 axle single trailers +# F10 FHWA Vehicle Class 10 6 plus axle single trailers +# F11 FHWA Vehicle Class 11 5 or less axle multi-trailers +# F12 FHWA Vehicle Class 12 6 axle multi-trailers +# F13 FHWA Vehicle Class 13 7 or more axle multi-trailers + +# goal #### +# The goal is to have volumes by MOVES categories +# Then, need toy data to proceed +# We can only be sure for PC = V19 + v20 and trucks v16 + v17 + v18 +# HOV3 seems to be mostly buses, we can assume that for now v11 + v12 + v13 + v14 + v15 + +# the other volumes can have a total distribution +sapply(st_set_geometry(net[, c("V19T_1", "V20T_1", "F2", "F1")], NULL), sum) + +metadata <- readxl::read_excel("../excel_inventory/inventory_MD.xlsx", "metadata") +setDT(metadata) +pop <- readxl::read_excel("../excel_inventory/inventory_MD.xlsx", "pop_moves") +setDT(pop) + +# Motorcycle = F1 +# Passenger Car = F2 +# Passenger Truck = F3 +# Light Commercial Truck = F3 +# Intercity Bus = F4 +# Transit Bus = F4 +# School Bus = F4 +# Refuse Truck = F5 + F6 + F8 +# Single Unit Short-haul Truck = F5 + F6 + F8 +# Single Unit Long-haul Truck = F7 + F9 + F10 +# Motor Home = F5 + F6 + F8 +# Combination Short-haul Truck = F11 +# Combination Long-haul Truck = F13 + F12 + +# assuming vols 1 to 10 covers all vehicles +# assumes trucks only on SOV, which makes sense +net$ALL <- net$V1T_1 + net$V2T_1 + net$V3T_1 + net$V4T_1 + net$V5T_1 +efes <- st_set_geometry(net[, paste0("F", 1:13)], NULL) +setDT(efes) +efes <- efes[rowSums(efes) > 0] + +col_efes <- colSums(efes)/sum(colSums(efes)) +df <- as.data.table(matrix(col_efes, nrow = 1)) +names(df) <- paste0("F", 1:13) + +df + +unique(pop$sourceTypeID) +MC_PC <- pop[sourceTypeID %in% c(11, 21), sum(sourceTypePopulation)] +per_MC <- pop[sourceTypeID %in% c(11), sum(sourceTypePopulation)]/MC_PC +per_PC <- pop[sourceTypeID %in% c(21), sum(sourceTypePopulation)]/MC_PC +# assumes V19T + V6T ...V10T = MC + PC +# assumes net$ALL*dfF1 = MC +# assumes net$ALL*dfF2 = PC + +net$MC <- (net$V19T_1 + net$V6T_1 +net$V7T_1 +net$V8T_1 +net$V9T_1 +net$V10T_1)*per_MC + net$ALL*df$F1 +net$PC <- (net$V19T_1 + net$V6T_1 +net$V7T_1 +net$V8T_1 +net$V9T_1 +net$V10T_1)*per_PC + net$ALL*df$F2 + + +PT_LCT <- pop[sourceTypeID %in% c(31, 32), sum(sourceTypePopulation)] +per_PT <- pop[sourceTypeID %in% c(31), sum(sourceTypePopulation)]/PT_LCT +per_LCT <- pop[sourceTypeID %in% c(31), sum(sourceTypePopulation)]/PT_LCT + +# assuming V20T are PT +net$PT <- (net$V20T_1 + net$ALL*df$F3)*per_PT +net$LCT <- (net$V20T_1 + net$ALL*df$F3)*per_LCT + +BUS <- pop[sourceTypeID %in% 41:43, sum(sourceTypePopulation)] +per_IBUS <- pop[sourceTypeID %in% c(41), sum(sourceTypePopulation)]/BUS +per_TBUS <- pop[sourceTypeID %in% c(42), sum(sourceTypePopulation)]/BUS +per_SBUS <- pop[sourceTypeID %in% c(43), sum(sourceTypePopulation)]/BUS +net$BUS_INTERCITY <- (net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1)*per_IBUS + net$ALL*df$F4*per_IBUS +per_IBUS*(net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1 + net$ALL*df$F4) + +net$BUS_TRANSIT <- (net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1)*per_TBUS+ net$ALL*df$F4*per_TBUS + +net$BUS_SCHOOL <- (net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1)*per_SBUS+ net$ALL*df$F4*per_SBUS + + +TRUCKS <- pop[sourceTypeID %in% c(51:54, 61, 62), sum(sourceTypePopulation)] +per_T_REFUSE <- pop[sourceTypeID %in% c(51), sum(sourceTypePopulation)]/TRUCKS +per_T_SU_SH <- pop[sourceTypeID %in% c(52), sum(sourceTypePopulation)]/TRUCKS +per_T_SU_LH <- pop[sourceTypeID %in% c(53), sum(sourceTypePopulation)]/TRUCKS +per_T_MH <- pop[sourceTypeID %in% c(54), sum(sourceTypePopulation)]/TRUCKS +per_T_CU_SH <- pop[sourceTypeID %in% c(61), sum(sourceTypePopulation)]/TRUCKS +per_T_CU_LH <- pop[sourceTypeID %in% c(62), sum(sourceTypePopulation)]/TRUCKS +# assumes trucks + sum V15..V18 +# assumes net$ALL (SOV) covers trucks F7...F13 +net$TRUCKS_REFUSE <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_REFUSE + +net$TRUCKS_SU_SH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_SU_SH + +net$TRUCKS_SU_LH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_SU_LH + +net$TRUCKS_MH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_MH + + +net$TRUCKS_CU_SH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_CU_SH +net$TRUCKS_CU_LH <- (net$V16T_1 + net$V17T_1 + net$V18T_1+ + net$ALL*(df$F5 + df$F6 + df$F7 + df$F8 + df$F9 + df$F10 + df$F11 + df$F12 + df$F13))*per_T_CU_LH + +file.remove("network/net.gpkg") +st_write(net, "network/net.gpkg") +# plot(net["PC"], axes = T, pal = cpt(colorRampPalette = T, rev = T), key.pos = 1) +# edit inventory.xlsx +rm(list = ls()) +gc() diff --git a/projects/moves/network/pre_net_simu_BMC_qgisv2.R b/projects/moves/network/pre_net_simu_BMC_qgisv2.R new file mode 100755 index 00000000..3a8b1bf0 --- /dev/null +++ b/projects/moves/network/pre_net_simu_BMC_qgisv2.R @@ -0,0 +1,80 @@ +library(data.table) +library(sf) +library(cptcity) +library(ggplot2) +# netemail <- sf::st_read("../vein_baltimore/network/net.gpkg") +# netemail <- st_transform(netemail, 4326) +net <- sf::st_read("network/net_qgis.gpkg") + +# V19 = PC + MC +# V20 = PT +# V16 17 18 TRUCKS +# V1 ... 15 ALL (without buses) + +# the other volumes can have a total distribution +sapply(st_set_geometry(net[, c("V19T_1", "V20T_1", "F2", "F1")], NULL), sum) + +metadata <- readxl::read_excel("config/inventory_MD.xlsx", "metadata") +setDT(metadata) +pop <- readxl::read_excel("config/inventory_MD.xlsx", "pop_moves") +setDT(pop) + +net$ALL <- net$V1T_1 + net$V2T_1 + net$V3T_1 + net$V4T_1 + net$V5T_1 + + net$V6T_1 + net$V7T_1 + net$V8T_1 + net$V9T_1 + net$V10T_1 + + net$V11T_1 + net$V12T_1 + net$V13T_1 + net$V14T_1 + net$V15T_1 + + +tot_no_bus <- pop[HPMSVtypeID != 40, sum(sourceTypePopulation)] +tot_bus <- pop[HPMSVtypeID == 40, sum(sourceTypePopulation)] +per_MC_ALL <- pop[HPMSVtypeID == 10]$sourceTypePopulation/tot_no_bus +per_PC_ALL <- pop[sourceTypeID == 21]$sourceTypePopulation/tot_no_bus +per_PT_ALL <- pop[sourceTypeID == 31]$sourceTypePopulation/tot_no_bus +per_LCT_ALL <- pop[sourceTypeID == 32]$sourceTypePopulation/tot_no_bus +per_REFUSE_TRUCKS_ALL <- pop[sourceTypeID == 51]$sourceTypePopulation/tot_no_bus +per_SUSH_TRUCKS_ALL <- pop[sourceTypeID == 52]$sourceTypePopulation/tot_no_bus +per_SULH_TRUCKS_ALL <- pop[sourceTypeID == 53]$sourceTypePopulation/tot_no_bus +per_MH_TRUCKS_ALL <- pop[sourceTypeID == 54]$sourceTypePopulation/tot_no_bus +per_CUSH_TRUCKS_ALL <- pop[sourceTypeID == 61]$sourceTypePopulation/tot_no_bus +per_CULH_TRUCKS_ALL <- pop[sourceTypeID == 62]$sourceTypePopulation/tot_no_bus + + +# MC PC #### +MC_PC <- pop[sourceTypeID %in% c(11, 21), sum(sourceTypePopulation)] +per_MC <- pop[sourceTypeID %in% c(11), sum(sourceTypePopulation)]/MC_PC +per_PC <- pop[sourceTypeID %in% c(21), sum(sourceTypePopulation)]/MC_PC + +net$MC <- net$V19T_1*per_MC + net$ALL*per_MC_ALL +net$PC <- net$V19T_1*per_PC + net$ALL*per_PC_ALL + +# PT LCT #### +PT_LCT <- pop[sourceTypeID %in% c(31, 32), sum(sourceTypePopulation)] +per_PT <- pop[sourceTypeID %in% c(31), sum(sourceTypePopulation)]/PT_LCT +per_LCT <- pop[sourceTypeID %in% c(32), sum(sourceTypePopulation)]/PT_LCT + +net$PT <- net$V20T_1*per_PT + net$ALL*per_PT_ALL +net$LCT <- net$V20T_1*per_LCT + net$ALL*per_LCT_ALL + +# TRUCKS #### +TRUCKS <- pop[sourceTypeID %in% c(51:54, 61, 62), sum(sourceTypePopulation)] +per_T_REFUSE <- pop[sourceTypeID %in% c(51), sum(sourceTypePopulation)]/TRUCKS +per_T_SU_SH <- pop[sourceTypeID %in% c(52), sum(sourceTypePopulation)]/TRUCKS +per_T_SU_LH <- pop[sourceTypeID %in% c(53), sum(sourceTypePopulation)]/TRUCKS +per_T_MH <- pop[sourceTypeID %in% c(54), sum(sourceTypePopulation)]/TRUCKS +per_T_CU_SH <- pop[sourceTypeID %in% c(61), sum(sourceTypePopulation)]/TRUCKS +per_T_CU_LH <- pop[sourceTypeID %in% c(62), sum(sourceTypePopulation)]/TRUCKS + +net$TRUCKS_REFUSE <- (net$V16T_1 + net$V17T_1 + net$V18T_1)*per_T_REFUSE + net$ALL*per_REFUSE_TRUCKS_ALL +net$TRUCKS_SU_SH <- (net$V16T_1 + net$V17T_1 + net$V18T_1)*per_T_SU_SH + net$ALL*per_SUSH_TRUCKS_ALL +net$TRUCKS_SU_LH <- (net$V16T_1 + net$V17T_1 + net$V18T_1)*per_T_SU_LH + net$ALL*per_SULH_TRUCKS_ALL +net$TRUCKS_MH <- (net$V16T_1 + net$V17T_1 + net$V18T_1)*per_T_MH + net$ALL*per_MH_TRUCKS_ALL +net$TRUCKS_CU_SH <- (net$V16T_1 + net$V17T_1 + net$V18T_1)*per_T_CU_SH + net$ALL*per_CUSH_TRUCKS_ALL +net$TRUCKS_CU_LH <- (net$V16T_1 + net$V17T_1 + net$V18T_1)*per_T_CU_LH + net$ALL*per_CULH_TRUCKS_ALL + + +file.remove("network/net_no_bus.gpkg") +st_write(net, "network/net_no_bus.gpkg") +# plot(net["PC"], axes = T, pal = cpt(colorRampPalette = T, rev = T), key.pos = 1) +# edit inventory.xlsx +rm(list = ls()) +gc() + diff --git a/projects/moves/scripts/emissions_per_distance.R b/projects/moves/scripts/emissions_per_distance.R new file mode 100755 index 00000000..c99ab59c --- /dev/null +++ b/projects/moves/scripts/emissions_per_distance.R @@ -0,0 +1,191 @@ + +file.remove(list.files(path = "emi/emissions_rate_per_distance/", full.names = T)) + + +# switch(language, +# "portuguese" = cat("Filtrando carro electrico: \n"), +# "english" = cat("Filtering electric car\n"), +# "spanish" = cat("Filtrando autos electricos\n") +# ) +setDT(metadata) +metadata_original <- metadata + +# me parece que la mejor opcion es filtrar por +# contaminante y por proceso. +# 1) por ejemplo, leer EF para pollutantID = 2 +# 2) luego filtrar por processID y relacionar con metadata +# 3) Luego, leer los veiculos con metadata +# 4) y estimar las eisiones para el processID + +# 1) Identify pollutants in data.base + +# estimation #### +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + +RMariaDB::dbGetQuery(con,paste0(" +SELECT DISTINCT pollutantID +FROM `rateperdistance` +WHERE rateperdistance > 0 +")) -> unique_pol_all +RMariaDB::dbDisconnect(con) + +setDT(unique_pol_all) +unique_pol <- unique_pol_all[!pollutantID %in% c(91, #ENERGY + 1, #THC + 118)] #NMHC + +# we keep pollutant 70 in case the user wants to applies a different speciation +# 79 NMHC, for instance, MOVES + +data(decoder) +# 2) Obtener EF + +# 2) Emissions per distance #### + +for(i in 1:nrow(unique_pol)) { + + cat("Detecting pollutants and process in SQL ...\n") + + # pollutantID = i + con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + + RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, processID, +modelYearID, hourID, avgSpeedBinID, ratePerDistance, AVG(ratePerDistance) +FROM `rateperdistance` +WHERE pollutantID = ", unique_pol$pollutantID[i], " AND +rateperdistance > 0 +GROUP BY MOVESRunID, pollutantID, modelYearID, processID, +fuelTypeID, avgSpeedBinID, hourID, sourceTypeID +")) -> ef + + names(ef)[ncol(ef)] <- "EF" + setDT(ef) + + RMariaDB::dbDisconnect(con) + + cat(paste0("Estimating emissions ...\n", + round(i/nrow(unique_pol)*100, 2), " %\n\n" )) + + moves_rpdy_meta( + metadata = metadata, + lkm = lkm, + ef = ef, + fuel_type = fuel_type, + speed_bin = speed_bin, + profile = tfs, + simplify = T) -> emi + + + data.table::fwrite(emi$streets, + paste0("emi/emissions_rate_per_distance/STREETS_pollutantID_", + unique_pol$pollutantID[i], + ".csv.gz"), + append = T) + + data.table::fwrite(emi$veh, "emi/emissions_rate_per_distance/DF.csv.gz", append = T) + rm(ef, emi) + gc() +} + +# 3) tire and break wear emissions for electric cars #### + + +switch(language, + "portuguese" = message("\nEmissões de freio e pneus de veículos elétricos\n"), + "english" = message("Vehicular emissions of electric cars for tire and break wear\n"), + "spanish" = message("Emisiones de vehículos electricos por desgaste de frenos y neumáticos\n") +) + + +metadata <- metadata_original[fuel == "ELEC"] + +for(i in seq_along(metadata$vehicles)) { + + + if(verbose) cat("Detecting pollutants and process in SQL ...\n") + + con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + + RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, processID, +modelYearID, hourID, avgSpeedBinID, ratePerDistance, AVG(ratePerDistance) +FROM `rateperdistance` +WHERE rateperdistance > 0 AND +sourceTypeID = ", metadata$sourceTypeID[i]," AND +fuelTypeID = 1 AND +(processID = 9 OR processID = 10) +GROUP BY MOVESRunID, pollutantID, modelYearID, +fuelTypeID, avgSpeedBinID, hourID, sourceTypeID, processID +")) -> ef + + setDT(ef) + names(ef)[ncol(ef)] <- "EF" + RMariaDB::dbDisconnect(con) + + + + cat(paste0("\nEstimating emissions ...\n", + "sourceTypeID: ", metadata$sourceTypeID[i], "\n", + "fuelTypeID: ", metadata$fuelTypeID[i], "\n", + round(i/nrow(metadata)*100, 2), " %\n" )) + + moves_rpdy_sf( + veh = readRDS(paste0("veh/", metadata$vehicles[i], ".rds")), + lkm = lkm, + ef = ef, + speed_bin = speed_bin, + profile = tfs[[metadata$vehicles[i]]], + source_type_id = metadata$sourceTypeID[i], + vehicle = metadata$vein_name[i], + vehicle_type = metadata$name[i], + fuel_subtype = metadata$fuel[i] + ) -> emi + + data.table::fwrite(emi$streets, + paste0("emi/emissions_rate_per_distance/STREETS_", + metadata$vehicles[i], + ".csv.gz"), + append = T) + data.table::fwrite(emi$veh, "emi/emissions_rate_per_distance/DF_ELEC.csv.gz", append = T) + +} + + +switch(language, + "portuguese" = message("\n\nArquivos em: /emi/*:"), + "english" = message("\nFiles in: /emi/*"), + "spanish" = message("\nArchivos en: /emi/*") +) + + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + +suppressWarnings( + rm( + i, j, pol, + n_PC, n_LCV, n_TRUCKS, n_BUS, n_MC, + ns, ln, p, df, dl, cores + ) +) + +invisible(gc()) diff --git a/projects/moves/scripts/emissions_per_distance_NMHC.R b/projects/moves/scripts/emissions_per_distance_NMHC.R new file mode 100755 index 00000000..ac481bbc --- /dev/null +++ b/projects/moves/scripts/emissions_per_distance_NMHC.R @@ -0,0 +1,139 @@ + +file.remove("emi/emissions_rate_per_distance/STREETS_pollutantID_NMHC_BY_TYPE.csv.gz") + +# switch(language, +# "portuguese" = cat("Filtrando carro electrico: \n"), +# "english" = cat("Filtering electric car\n"), +# "spanish" = cat("Filtrando autos electricos\n") +# ) +setDT(metadata) +metadata_original <- metadata + +# me parece que la mejor opcion es filtrar por +# contaminante y por proceso. +# 1) por ejemplo, leer EF para pollutantID = 2 +# 2) luego filtrar por processID y relacionar con metadata +# 3) Luego, leer los veiculos con metadata +# 4) y estimar las eisiones para el processID + +# select unique fuels #### +cat("Detecting unique fuels for NMHC ...\n") + +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + +RMariaDB::dbGetQuery(con,paste0(" +SELECT DISTINCT fuelTypeID +FROM `rateperdistance` +WHERE pollutantID = 79 AND +rateperdistance > 0 +")) -> fuelTypeID + +names(fuelTypeID)[ncol(fuelTypeID)] <- "EF" +setDT(fuelTypeID) + +RMariaDB::dbDisconnect(con) + +for(i in seq_along(fuelTypeID$EF)) { + +cat(paste0("Detecting NMHC and process for fuelTypeID ", fuelTypeID$EF[i], " in SQL ...\n")) +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + +RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, processID, +modelYearID, hourID, avgSpeedBinID, ratePerDistance, AVG(ratePerDistance) +FROM `rateperdistance` +WHERE pollutantID = 79 AND +fuelTypeID = ", fuelTypeID$EF[i], " AND +rateperdistance > 0 +GROUP BY MOVESRunID, pollutantID, modelYearID, processID, +fuelTypeID, avgSpeedBinID, hourID, sourceTypeID +")) -> ef + +names(ef)[ncol(ef)] <- "EF" +setDT(ef) + +RMariaDB::dbDisconnect(con) + +unique(ef$fuelT) + +moves_rpdy_meta( + metadata = metadata, + lkm = lkm, + ef = ef, + fuel_type = fuel_type, + speed_bin = speed_bin, + profile = tfs, + simplify = F) -> emi + + +metadata$vtype <- c(rep("LDV", 13), rep("HDV", 39-14)) +id_ldv <- unique(metadata[vtype == "LDV"]$sourceTypeID) +id_hdv <- unique(metadata[vtype == "HDV"]$sourceTypeID) + +emi[, family := fifelse( + sourceTypeID %in% id_ldv, "LDV", "HDV" + )] + +emi[, family := fifelse( + fuel %in% "CNG", "ALL", family +)] + +data(decoder) + +names(decoder$emission_process) <- c("processID", "process") + +emi <- merge(emi, + decoder$emission_process, + by = "processID", + all.x = TRUE) + +# we need +setDT(decoder$emission_process) + +pros <- decoder$emission_process[processID %in% unique(emi$processID)] + +# Assuming Refuelling as evap +exhaust <- unique(grep(pattern = "Exhaust", + x = emi$process, + value = T)) + +emi[, process2 := fifelse(process %in% exhaust, "exhaust", "evap")] + + +emi[, + sum(age_total), + by = .(id, hour, family, fuel, process2)] -> emi2 + +streets <- data.table::dcast.data.table(emi2, + formula = id +family + fuel + process2 ~ hour, + value.var = "V1") + +names(streets)[5:ncol(streets)] <- paste0("H", names(streets)[5:ncol(streets)]) + +# Need streets results by + +data.table::fwrite(streets, + "emi/emissions_rate_per_distance/STREETS_pollutantID_NMHC_BY_TYPE.csv.gz", + append = T) + +} + +suppressWarnings( + rm( + i, j, pol, + n_PC, n_LCV, n_TRUCKS, n_BUS, n_MC, + ns, ln, p, df, dl, cores, streets, emi, emi2 + ) +) + +invisible(gc()) diff --git a/projects/moves/scripts/emissions_per_vehicle.R b/projects/moves/scripts/emissions_per_vehicle.R new file mode 100755 index 00000000..cb4f1576 --- /dev/null +++ b/projects/moves/scripts/emissions_per_vehicle.R @@ -0,0 +1,127 @@ + +file.remove(list.files(path = "emi/emissions_rate_per_vehicle/", full.names = T)) + + +# switch(language, +# "portuguese" = cat("Filtrando carro electrico: \n"), +# "english" = cat("Filtering electric car\n"), +# "spanish" = cat("Filtrando autos electricos\n") +# ) +setDT(metadata) +metadata_original <- metadata + +# me parece que la mejor opcion es filtrar por +# contaminante y por proceso. +# 1) por ejemplo, leer EF para pollutantID = 2 +# 2) luego filtrar por processID y relacionar con metadata +# 3) Luego, leer los veiculos con metadata +# 4) y estimar las eisiones para el processID + +# 1) Identify pollutants in data.base + +metadata <- metadata_original[fuel != "ELEC"] + +metadata$km_trip <- units::set_units(metadata$km_trip, "km") +metadata$km_trip <- units::set_units(metadata$km_trip, "miles") + +# estimation #### +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + +RMariaDB::dbGetQuery(con,paste0(" +SELECT DISTINCT pollutantID +FROM `ratepervehicle` +WHERE ratepervehicle > 0 AND +pollutantID != 91 +")) -> unique_pol_all +RMariaDB::dbDisconnect(con) + +setDT(unique_pol_all) +unique_pol <- unique_pol_all[!pollutantID %in% c(91, #ENERGY + 1, #THC + 118)] #Composite + +# we keep pollutant 70 in case the user wants to applies a different speciation +# 79 NMHC, for instance, MOVES + +data(decoder) +# 2) Obtener EF +for(pols in 1:nrow(unique_pol)) { + + cat("Detecting pollutants and process in SQL ...\n") + + # pollutantID = i + con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + + RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, processID, +modelYearID, hourID, ratePerVehicle, AVG(ratePerVehicle) +FROM `rateperVehicle` +WHERE pollutantID = ",unique_pol$pollutantID[pols], " AND +ratepervehicle > 0 +GROUP BY MOVESRunID, pollutantID, modelYearID, processID, +fuelTypeID, hourID, sourceTypeID +")) -> ef + + names(ef)[ncol(ef)] <- "EF" + setDT(ef) + + RMariaDB::dbDisconnect(con) + + cat(paste0("Estimating emissions ...\n", + round(pols/nrow(unique_pol)*100, 2), " %\n" )) + + moves_rpsy_meta( + metadata = metadata, + lkm = rep(1, length(lkm)), + ef = ef, + fuel_type = fuel_type, + profile = tfs, + colkt = F, + simplify = T, + verbose = verbose) -> emi + + + data.table::fwrite(emi$streets, + paste0("emi/emissions_rate_per_vehicle/STREETS_pollutantID_", + unique_pol$pollutantID[pols], + ".csv.gz"), + append = F) + + data.table::fwrite(emi$veh, "emi/emissions_rate_per_vehicle/DF.csv.gz", append = T) + rm(ef, emi) + gc() +} + + +switch(language, + "portuguese" = message("\n\nArquivos em: /emi/*:"), + "english" = message("\nFiles in: /emi/*"), + "spanish" = message("\nArchivos en: /emi/*") +) + + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + +suppressWarnings( + rm( + i, j, pol, + n_PC, n_LCV, n_TRUCKS, n_BUS, n_MC, + ns, ln, p, df, dl, cores + ) +) + +invisible(gc()) diff --git a/projects/moves/scripts/emissions_per_vehicle_NMHC.R b/projects/moves/scripts/emissions_per_vehicle_NMHC.R new file mode 100755 index 00000000..3269ffdb --- /dev/null +++ b/projects/moves/scripts/emissions_per_vehicle_NMHC.R @@ -0,0 +1,166 @@ + +file.remove("emi/emissions_rate_per_vehicle/STREETS_pollutantID_NMHC_BY_TYPE.csv.gz") + +# switch(language, +# "portuguese" = cat("Filtrando carro electrico: \n"), +# "english" = cat("Filtering electric car\n"), +# "spanish" = cat("Filtrando autos electricos\n") +# ) +setDT(metadata) +metadata_original <- metadata + +# me parece que la mejor opcion es filtrar por +# contaminante y por proceso. +# 1) por ejemplo, leer EF para pollutantID = 2 +# 2) luego filtrar por processID y relacionar con metadata +# 3) Luego, leer los veiculos con metadata +# 4) y estimar las eisiones para el processID + +# 1) Identify pollutants in data.base + +metadata <- metadata_original[fuel != "ELEC"] + +# we keep pollutant 70 in case the user wants to applies a different speciation +# 79 NMHC, for instance, MOVES + +data(decoder) +# 2) Obtener EF + + +# select unique fuels #### +cat("Detecting unique fuels for NMHC ...\n") + +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + +RMariaDB::dbGetQuery(con,paste0(" +SELECT DISTINCT fuelTypeID +FROM `rateperdistance` +WHERE pollutantID = 79 AND +rateperdistance > 0 +")) -> fuelTypeID + +names(fuelTypeID)[ncol(fuelTypeID)] <- "EF" +setDT(fuelTypeID) + +RMariaDB::dbDisconnect(con) + +for(i in seq_along(fuelTypeID$EF)) { + +cat("Detecting pollutants and process in SQL ...\n") + +# pollutantID = i +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + +RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, processID, +modelYearID, hourID, ratePerVehicle, AVG(ratePerVehicle) +FROM `rateperVehicle` +WHERE pollutantID = 79 AND +fuelTypeID = ", fuelTypeID$EF[i], " AND +ratepervehicle > 0 +GROUP BY MOVESRunID, pollutantID, modelYearID, processID, +fuelTypeID, hourID, sourceTypeID +")) -> ef + +names(ef)[ncol(ef)] <- "EF" +setDT(ef) + +RMariaDB::dbDisconnect(con) + + +moves_rpsy_meta( + metadata = metadata, + lkm = rep(1, length(lkm)), + ef = ef, + fuel_type = fuel_type, + profile = tfs, + colkt = F, + simplify = F, + verbose = verbose) -> emi + + +#NMHC +metadata <- readRDS("config/metadata.rds") +setDT(metadata) +metadata$vtype <- c(rep("LDV", 13), rep("HDV", 39-14)) +id_ldv <- unique(metadata[vtype == "LDV"]$sourceTypeID) +id_hdv <- unique(metadata[vtype == "HDV"]$sourceTypeID) + +emi[, family := fifelse( + sourceTypeID %in% id_ldv, "LDV", "HDV" +)] + +emi[, family := fifelse( + fuel %in% "CNG", "ALL", family +)] + +data(decoder) + +names(decoder$emission_process) <- c("processID", "process") + +emi <- merge(emi, + decoder$emission_process, + by = "processID", + all.x = TRUE) + +# we need +setDT(decoder$emission_process) +pros <- decoder$emission_process[processID %in% unique(emi$processID)] +# ssuming Refuelling as evap +exhaust <- unique(grep(pattern = "Exhaust", + x = emi$process, + value = T)) + +emi[, process2 := fifelse(process %in% exhaust, "exhaust", "evap")] + + +emi[, + sum(age_total, na.rm = T), + by = .(id, hour, family, fuel, process2)] -> emi2 + +streets <- data.table::dcast.data.table(emi2, + formula = id +family + fuel + process2 ~ hour, + value.var = "V1") + +names(streets)[5:ncol(streets)] <- paste0("H", names(streets)[5:ncol(streets)]) + +# Need streets results by + +data.table::fwrite(streets, + "emi/emissions_rate_per_vehicle/STREETS_pollutantID_NMHC_BY_TYPE.csv.gz", + append = T) + +} + +switch(language, + "portuguese" = message("\n\nArquivos em: /emi/*:"), + "english" = message("\nFiles in: /emi/*"), + "spanish" = message("\nArchivos en: /emi/*") +) + + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + +suppressWarnings( + rm( + i, j, pol, + n_PC, n_LCV, n_TRUCKS, n_BUS, n_MC, + ns, ln, p, df, dl, cores + ) +) + +invisible(gc()) diff --git a/projects/moves/scripts/evaporatives.R b/projects/moves/scripts/evaporatives.R new file mode 100755 index 00000000..94375e69 --- /dev/null +++ b/projects/moves/scripts/evaporatives.R @@ -0,0 +1,359 @@ +suppressWarnings(file.remove("emi/table_evaporatives.csv")) + + +# temperature +te <- met$Temperature +bb <- (0 + 15) / 2 +cc <- (10 + 25) / 2 +dd <- (20 + 35) / 2 + +tem <- ifelse( + te <= bb, "0_15", + ifelse( + te > bb & te <= cc, "10_25", + "20_35" + ) +) + +nmonth <- ifelse(nchar(seq_along(te)) < 2, + paste0(0, seq_along(te)), + seq_along(te) +) + +# filtrando veiculos otto +meta_ev <- metadata[metadata$fuel != "D", ] +veh_ev <- meta_ev$vehicles + +# checks name +type_emis <- c("Diurnal", "Running Losses", "Hot Soak") +name_file_evap <- c("/DIURNAL_", "/RUNNING_LOSSES_", "/HOT_SOAK_") + +ef_d <- paste0("D_", tem) +ef_rl <- paste0("R_", tem) +ef_hs <- paste0("S_", tem) + +# plot +n_PC <- metadata[metadata$family == "PC", ]$vehicles +n_LCV <- metadata[metadata$family == "LCV", ]$vehicles[1:4] +n_MC <- metadata[metadata$family == "MC", ]$vehicles + +ns <- c( + "PC", "LCV", "MC", + "PC", "LCV", "MC", + "PC", "LCV", "MC" +) +ln <- list( + n_PC, n_LCV, n_MC, + n_PC, n_LCV, n_MC, + n_PC, n_LCV, n_MC +) +laby <- c( + "g/day", "g/day", "g/day", + "g/trip", "g/trip", "g/trip", + "g/trip", "g/trip", "g/trip" +) +ev <- c( + "DIURNAL", "DIURNAL", "DIURNAL", + "RUNNING_LOSSES", "RUNNING_LOSSES", "RUNNING_LOSSES", + "HOT_SOAK", "HOT_SOAK", "HOT_SOAK" +) + +# plotting +switch(language, + "portuguese" = cat("Plotando EF\n"), + "english" = cat("Plotting EF\n"), + "spanish" = cat("Plotando EF\n") +) + +for (i in seq_along(ns)) { + dl <- lapply(seq_along(ef_d), function(j) { + data.frame(ef_cetesb( + p = ef_d[j], veh = ln[[i]], year = year, + agemax = 40, verbose = verbose + ), month = nmonth[j]) + }) + + dl <- rbindlist(dl) + df <- wide_to_long( + df = dl, + column_with_data = ln[[i]], + column_fixed = "month" + ) + df$age <- 1:40 + setDT(df) + names(df) <- c("ef", "month", "veh", "age") + p <- ggplot( + df[df$ef > 0, ], + aes(x = age, y = ef, colour = veh) + ) + + geom_line() + + facet_wrap(~month) + + ylim(0, NA) + + labs(y = laby[i], title = ev[i]) + + # scale_y_log10() + + theme_bw() + + png( + filename = paste0("images/EF_", ev[i], "_", ns[i], ".png"), + width = 2100, height = 1500, units = "px", pointsize = 12, + bg = "white", res = 300 + ) + print(p) + dev.off() +} + +switch(language, + "portuguese" = message("\nFiguras em /images\n"), + "english" = message("\nFigures in /image\n"), + "spanish" = message("\nFiguras en /images\n") +) + + + +# Diurnal #### +switch(language, + "portuguese" = cat("\nEmissões evaporativas diurnal\n"), + "english" = cat("\nEvaporative diurnal emissions\n"), + "spanish" = cat("\nEmisiones evaporativas diurnal\n") +) + +for (i in seq_along(veh_ev)) { + cat( + "\n", veh_ev[i], + rep("", max(nchar(veh_ev) + 1) - nchar(veh_ev[i])) + ) + + + for (j in seq_along(te)) { + cat(nmonth[j], " ") + + x <- readRDS(paste0("veh/", veh_ev[i], ".rds")) + + ef <- ef_cetesb( + p = ef_d[j], + veh = veh_ev[i], + year = year, + agemax = ncol(x), + # diurnal: g/day* day/km= g/km + verbose = verbose + ) / (mileage[[veh_ev[i]]] / 365) # mean daily mileage + + # muda NaNaN para 0 + ef[is.na(ef)] <- 0 + + array_x <- emis( + veh = x, + lkm = lkm, + ef = ef, + profile = tfs[[veh_ev[i]]], + fortran = TRUE, + nt = check_nt() / 2, + simplify = TRUE, + verbose = verbose + ) + + x_DF <- emis_post( + arra = array_x, + veh = veh_ev[i], + size = metadata$size[i], + fuel = metadata$fuel[i], + pollutant = "NMHC", + type_emi = "Evaporative", + by = "veh" + ) + x_DF$subtype_emi <- "Diurnal" + + fwrite(x_DF, "emi/table_evaporatives.csv", append = TRUE) + + x_STREETS <- emis_post( + arra = array_x, + pollutant = veh_ev[j], + by = "streets" + ) + saveRDS(x_STREETS, + file = paste0( + "emi/", + veh_ev[i], + name_file_evap[1], + veh_ev[i], "_", + "EVAP", "_", + nmonth[j], "_", + "NMHC_STREETS.rds" + ) + ) + + rm(array_x, ef, x, x_DF, x_STREETS) + } +} + + +# Running Losses #### +switch(language, + "portuguese" = cat("\nEmissões evaporativas running-losses\n"), + "english" = cat("\nEvaporative running-losses emissions\n"), + "spanish" = cat("\nEmisiones evaporativas running-loses\n") +) + +for (i in seq_along(veh_ev)) { + cat( + "\n", veh_ev[i], + rep("", max(nchar(veh_ev) + 1) - nchar(veh_ev[i])) + ) + + + for (j in seq_along(te)) { + cat(nmonth[j], " ") + x <- readRDS(paste0("veh/", veh_ev[i], ".rds")) + + ef <- ef_cetesb( + p = ef_rl[j], + veh = veh_ev[i], + year = year, + agemax = ncol(x), + # g/trip * trip/day * day/km = g/km + verbose = verbose + ) * meta_ev$trips_day[i] / (mileage[[veh_ev[i]]] / 365) + + # muda NaNaN para 0 + ef[is.na(ef)] <- 0 + + + array_x <- emis( + veh = x, + lkm = lkm, + ef = ef, + profile = tfs[[veh_ev[i]]], + fortran = TRUE, + nt = check_nt() / 2, + simplify = TRUE, + verbose = verbose + ) + + x_DF <- emis_post( + arra = array_x, + veh = veh_ev[i], + size = metadata$size[i], + fuel = metadata$fuel[i], + pollutant = "NMHC", + type_emi = "Evaporative", + by = "veh" + ) + x_DF$subtype_emi <- "Running Losses" + + fwrite(x_DF, "emi/table_evaporatives.csv", append = TRUE) + + x_STREETS <- emis_post( + arra = array_x, + pollutant = veh_ev[j], + by = "streets" + ) + saveRDS(x_STREETS, + file = paste0( + "emi/", + veh_ev[i], + name_file_evap[2], + veh_ev[i], "_", + "EVAP", "_", + nmonth[j], "_", + "NMHC_STREETS.rds" + ) + ) + + rm(array_x, ef, x, x_DF, x_STREETS) + } +} + +# Hot Soak #### +switch(language, + "portuguese" = cat("\nEmissões evaporativas hot-soak\n"), + "english" = cat("\nEvaporative hot-soak emissions\n"), + "spanish" = cat("\nEmisiones evaporativas hot-soak\n") +) + +for (i in seq_along(veh_ev)) { + cat( + "\n", veh_ev[i], + rep("", max(nchar(veh_ev) + 1) - nchar(veh_ev[i])) + ) + + + for (j in seq_along(te)) { + cat(nmonth[j], " ") + x <- readRDS(paste0("veh/", veh_ev[i], ".rds")) + + ef <- ef_cetesb( + p = ef_hs[j], + veh = veh_ev[i], + year = year, + agemax = ncol(x), + verbose = verbose + ) * meta_ev$trips_day[i] / (mileage[[veh_ev[i]]] / 365) # quilometragem medio diario + + # muda NaNaN para 0 + ef[is.na(ef)] <- 0 + + + array_x <- emis( + veh = x, + lkm = lkm, + ef = ef, + profile = tfs[[veh_ev[i]]], + fortran = TRUE, + nt = check_nt() / 2, + simplify = TRUE, + verbose = verbose + ) + + x_DF <- emis_post( + arra = array_x, + veh = veh_ev[i], + size = metadata$size[i], + fuel = metadata$fuel[i], + pollutant = "NMHC", + type_emi = "Evaporative", + by = "veh" + ) + x_DF$subtype_emi <- "Hot Soak" + + fwrite(x_DF, "emi/table_evaporatives.csv", append = TRUE) + + x_STREETS <- emis_post( + arra = array_x, + pollutant = veh_ev[j], + by = "streets" + ) + saveRDS(x_STREETS, + file = paste0( + "emi/", + veh_ev[i], + name_file_evap[3], + veh_ev[i], "_", + "EVAP", "_", + nmonth[j], "_", + "NMHC_STREETS.rds" + ) + ) + + rm(array_x, ef, x, x_DF, x_STREETS) + } +} + +switch(language, + "portuguese" = message("\n\nArquivos em: /emi/*:"), + "english" = message("\nFiles in: /emi/*"), + "spanish" = message("\nArchivos en: /emi/*") +) + + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + +suppressWarnings(rm( + i, mileage, meta_ev, veh_ev, year, + diurnal_ef, hot_soak_ef, running_losses_ef +)) + +invisible(gc()) \ No newline at end of file diff --git a/projects/moves/scripts/fuel_eval.R b/projects/moves/scripts/fuel_eval.R new file mode 100755 index 00000000..d4cd4693 --- /dev/null +++ b/projects/moves/scripts/fuel_eval.R @@ -0,0 +1,347 @@ +switch(language, + "portuguese" = cat("Filtrando carro electrico: \n"), + "english" = cat("Filtering electric car\n"), + "spanish" = cat("Filtrando autos electricos\n") +) + +year +month <- fuel$Month[1] +date <- as.Date(ISOdate(year, month, 1,0,0)) + +numberOfDays <- function(date) { + m <- format(date, format="%m") + + while (format(date, format="%m") == m) { + date <- date + 1 + } + return(as.integer(format(date - 1, format="%d"))) +} + +factor_emi <- numberOfDays(date = date) # daily to month + +data(decoder) + +setDT(metadata) +metadata <- metadata[fuel != "ELEC"] + +# vehicles #### +lapply(seq_along(metadata$vehicles), function(kk) { + readRDS(paste0("veh/", metadata$vehicles[kk], ".rds")) +}) -> lv + +names(lv) <- metadata$vehicles + +#SQL EF PER DISTANCE #### +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + + +RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, +modelYearID, hourID, avgSpeedBinID, ratePerDistance, AVG(ratePerDistance) +FROM `rateperdistance` +WHERE pollutantID = 91 AND +processID = 1 AND +rateperdistance > 0 +GROUP BY MOVESRunID, pollutantID, modelYearID, +fuelTypeID, avgSpeedBinID, hourID, sourceTypeID +")) -> ef + +names(ef)[ncol(ef)] <- "EF" +setDT(ef) + +RMariaDB::dbDisconnect(con) + +# moves rates per distance year +moves_rpdy( + veh = lv, + lkm = lkm, + ef = ef, + source_type_id = metadata$sourceTypeID, + fuel_type_id = metadata$fuelTypeID, + pollutant_id = 91, + fuel_type = fuel_type, + speed_bin = speed_bin, + profile = tfs, + vehicle = metadata$vein_name, + vehicle_type = metadata$name, + fuel_subtype = metadata$fuel, + process_id = 1, + verbose = verbose +) -> emi + +#SQL EF PER VEHICLE #### + +cat("Detecting pollutants and process in SQL ...\n") + +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + +RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, processID, +modelYearID, hourID, ratePerVehicle, AVG(ratePerVehicle) +FROM `rateperVehicle` +WHERE pollutantID = 91 AND +ratepervehicle > 0 +GROUP BY MOVESRunID, pollutantID, modelYearID, processID, +fuelTypeID, hourID, sourceTypeID +")) -> ef + +names(ef)[ncol(ef)] <- "EF" +setDT(ef) + +RMariaDB::dbDisconnect(con) + +moves_rpsy_meta( + metadata = metadata, + lkm = rep(1, length(lkm)), + ef = ef, + fuel_type = fuel_type, + profile = tfs, + colkt = F, + simplify = T, + verbose = verbose) -> emiv + + + # merge #### +emiv$veh$value <- emiv$veh$value *units::set_units(1,"miles") +names(emiv$veh) +names(emi$veh) + + +emiv$veh$fuelTypeID <- emiv$veh$fuel + +# emiallback up +bu_emi <- emi$veh +bu_emi$db <- "emissions_per_distance" +bu_emiv <- emiv$veh +bu_emiv$db <- "emissions_per_vehicle" + +emi <- rbind(bu_emi, + bu_emiv[, names(bu_emi), with = F]) + +emi[, sum(value), by = hour][, plot(V1, type = "b")] +fc <- emi[, sum(value)*factor_emi, by = fuelTypeID] + +fuel <- fuel[, c("fuel", "gallons", "Year", "Month")] + + +names(fuel)[1:2] <- c("fuelTypeID","consumption_gallons") + +fuel$consumption_gallons <- units::set_units(fuel$consumption_gallons, "gallons") + +fuel <- merge(fc, fuel, by = "fuelTypeID", all.x = T) +names(fuel)[2] <- "estimated_gallons" + +fuel$estimation_consumption <- as.numeric(fuel$estimated_gallons)/as.numeric(fuel$consumption_gallons) +# dtf$estimation_consumption <- dtf$estimation_t / dtf$consumption_t +fuel <-fuel[fuelTypeID != "E85"] + +print(fuel) + +############################################################################ +# obtaining k #### +############################################################################ +k_D <- as.numeric(1/fuel[fuelTypeID == "D", ]$estimation_consumption) +k_CNG <- as.numeric(1/fuel[fuelTypeID == "CNG", ]$estimation_consumption) +k_G <- as.numeric(1/fuel[fuelTypeID == "G", ]$estimation_consumption) +k_E85 <- 1 + +saveRDS(list(k_D = k_D, + k_CNG = k_CNG, + k_G = k_G), + "veh/list_k_fuel.rds") + +language <- "english" # portuguese spanish +metadata <- readRDS("config/metadata.rds") +net <- readRDS("network/net.rds") +tfs <- readRDS("config/tfs.rds") +veh <- readRDS("config/fleet_age.rds") +categories <- unique(metadata$vein_name) # in net +year <- 2019 +theme <- "dark" # dark clean ink + +verbose <- FALSE +speed <- FALSE +fuel <- readRDS("config/fuel.rds") +source("scripts/traffic.R", encoding = "UTF-8") + +############################################################################ +# Re estimating FC #### +############################################################################ +factor_emi <- numberOfDays(date = date) # daily to month + +setDT(metadata) +metadata <- metadata[fuel != "ELEC"] + +lapply(seq_along(metadata$vehicles), function(kk) { + readRDS(paste0("veh/", metadata$vehicles[kk], ".rds")) +}) -> lv + +names(lv) <- metadata$vehicles + + +#SQL EF PER DISTANCE #### +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + + +RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, +modelYearID, hourID, avgSpeedBinID, ratePerDistance, AVG(ratePerDistance) +FROM `rateperdistance` +WHERE pollutantID = 91 AND +processID = 1 AND +rateperdistance > 0 +GROUP BY MOVESRunID, pollutantID, modelYearID, +fuelTypeID, avgSpeedBinID, hourID, sourceTypeID +")) -> ef + +names(ef)[ncol(ef)] <- "EF" +setDT(ef) + +RMariaDB::dbDisconnect(con) + +# moves rates per distance year +moves_rpdy( + veh = lv, + lkm = lkm, + ef = ef, + source_type_id = metadata$sourceTypeID, + fuel_type_id = metadata$fuelTypeID, + pollutant_id = 91, + fuel_type = fuel_type, + speed_bin = speed_bin, + profile = tfs, + vehicle = metadata$vein_name, + vehicle_type = metadata$name, + fuel_subtype = metadata$fuel, + process_id = 1, + verbose = verbose +) -> emi + + +#SQL EF PER VEHICLE #### + +cat("Detecting pollutants and process in SQL ...\n") + +con <- DBI::dbConnect(RMariaDB::MariaDB(), + host = host, + db = moves_db, + user = user, + password = password, + port = port) + +RMariaDB::dbGetQuery(con,paste0(" +SELECT pollutantID, sourceTypeID, fuelTypeID, processID, +modelYearID, hourID, ratePerVehicle, AVG(ratePerVehicle) +FROM `rateperVehicle` +WHERE pollutantID = 91 AND +ratepervehicle > 0 +GROUP BY MOVESRunID, pollutantID, modelYearID, processID, +fuelTypeID, hourID, sourceTypeID +")) -> ef + +names(ef)[ncol(ef)] <- "EF" +setDT(ef) + +RMariaDB::dbDisconnect(con) + +moves_rpsy_meta( + metadata = metadata, + lkm = rep(1, length(lkm)), + ef = ef, + fuel_type = fuel_type, + profile = tfs, + colkt = F, + simplify = T, + verbose = verbose) -> emiv + + +# merge #### +emiv$veh$value <- emiv$veh$value *units::set_units(1,"miles") +names(emiv$veh) +names(emi$veh) + + +emiv$veh$fuelTypeID <- emiv$veh$fuel + +# emiallback up +bu_emi <- emi$veh +bu_emi$db <- "emissions_per_distance" +bu_emiv <- emiv$veh +bu_emiv$db <- "emissions_per_vehicle" + +emi <- rbind(bu_emi, + bu_emiv[, names(bu_emi), with = F]) + +fc <- emi$veh[, sum(value)*factor_emi, by = fuelTypeID] + +fuel <- fuel[, c("fuel", "gallons", "Year", "Month")] + + +names(fuel)[1:2] <- c("fuelTypeID","consumption_gallons") + +fuel$consumption_gallons <- units::set_units(fuel$consumption_gallons, "gallons") + +fuel <- merge(fc, fuel, by = "fuelTypeID", all.x = T) +names(fuel)[2] <- "estimated_gallons" + +fuel$estimation_consumption <- as.numeric(fuel$estimated_gallons)/as.numeric(fuel$consumption_gallons) +# dtf$estimation_consumption <- dtf$estimation_t / dtf$consumption_t +fuel <-fuel[fuelTypeID != "E85"] + +print(fuel) + + +switch(language, + "portuguese" = message("\nArquivos em: /emi/*:"), + "english" = message("\nFiles in: /emi/*"), + "spanish" = message("\nArchivos en: /emi/*") +) + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + +suppressWarnings(rm(i, j, pol, dt, dt0, dtf, factor_emi, fuel, emi)) + +suppressWarnings( + rm( + kPC, kPC_G, kPC_E, kPC_FG, kPC_FE, + kLCV, kLCV_G, kLCV_E, kLCV_FG, kLCV_FE, kLCV_D, + kTRUCKS, kTRUCKS_SL_D, kTRUCKS_L_D, kTRUCKS_M_D, kTRUCKS_SH_D, kTRUCKS_H_D, + kBUS, kBUS_URBAN_D, kBUS_MICRO_D, kBUS_COACH_D, + kMC, kMC_150_G, kMC_150_500_G, kMC_500_G, + kMC_150_FG, kMC_150_500_FG, kMC_500_FG, + kMC_150_FE, kMC_150_500_FE, kMC_500_FE, + l_PC, l_LCV, l_TRUCKS, l_BUS, l_MC, + i, arquivos, cores, df, f, ff, + n_PC, n_LCV, n_TRUCKS, n_BUS, n_MC, + na, nveh, p, tit, tit2, categories, verbose, x, kf, + k_G, k_D, k_E, + metadata, + net, veh, year, + theme, + speed, + sp_mov, + df, + dfs + ) +) + +invisible(gc()) diff --git a/projects/moves/scripts/mech.R b/projects/moves/scripts/mech.R new file mode 100755 index 00000000..2fa93b68 --- /dev/null +++ b/projects/moves/scripts/mech.R @@ -0,0 +1,339 @@ + +unlink(paste0("post/", mech), recursive = T) + +dir.create(paste0("post/", mech), showWarnings = FALSE) + +fs <- list.files( + path = paste0("post/", mech), + pattern = ".rds", + full.names = TRUE, + recursive = TRUE +) +file.remove(fs) + +# ntfs +ntfs <- paste0("H", 1:nrow(tfs)) + +# NMHC +dfnmhc <- st_set_geometry(readRDS("post/streets/NMHC.rds"), NULL) +setDT(dfnmhc) + +# check +dfnmhc[, unique(cbind(family, fuel, process2))] + +# vein does not have speciation for heavy vehicles using gasoline +# we assumed light vehicles +dfnmhc[family == "HDV" & fuel == "G"]$family <- "LDV" + +# check +dfnmhc[, unique(cbind(family, fuel, process2))] + +dfnmhc[family == "HDV" & fuel == "D", unique(process2)] + +# vein does not have speciation for heavy vehicles evaporated +# we assumed exhaust +dfnmhc[fuel == "D" & process2 == "evap"]$process2 <- "exhaust" + +# check +dfnmhc[, unique(cbind(family, fuel, process2))] -> check + +if(nrow(check) == 7) print("Emissions adjusted") + +# updating the emissions +dfnmhc[, + lapply(.SD, sum, na.rm = T), + .SDcols = ntfs, + by = .(id, family, fuel, process2)] -> dfnmhc + +nrow(dfnmhc)/nrow(net) +dfnmhc[, unique(cbind(family, fuel, process2))] + +# 1) LDV Diesel Exhaust #### +x <- dfnmhc[family == "LDV" & + fuel == "D", + c("id", ntfs), with = F] + +x[is.na(x)] <- 0 +id <- x$id +x$id <- NULL + +# nx!!!! +nx <- ntfs + +dx <- speciate( + x = x, + spec = "nmhc", + fuel = "D", + veh = "LDV", + eu = "all" +) + +dx$id <- rep(id, length(unique(dx$pol))) + +nmhc_ldv_d_exh <- emis_chem2( + df = dx, + mech = mech, + nx = nx, + na.rm = TRUE +) + +# 2) HDV Diesel Exhaust #### +x <- dfnmhc[family == "HDV" & + fuel == "D", + c("id", ntfs), with = F] + +x[is.na(x)] <- 0 +id <- x$id +x$id <- NULL + +# nx!!!! +nx <- ntfs + +dx <- speciate( + x = x, + spec = "nmhc", + fuel = "D", + veh = "HDV", + eu = "all" +) + +dx$id <- rep(id, length(unique(dx$pol))) + +nmhc_hdv_d_exh <- emis_chem2( + df = dx, + mech = mech, + nx = nx, + na.rm = TRUE +) + + +# 3) LDV Gasoline Exhaust #### +x <- dfnmhc[family == "LDV" & + fuel == "G" & + process2 == "exhaust", + c("id", ntfs), with = F] + +x[is.na(x)] <- 0 +id <- x$id +x$id <- NULL + +# nx!!!! +nx <- ntfs + +dx <- speciate( + x = x, + spec = "nmhc", + fuel = "G", + veh = "LDV", + eu = "I" +) + +dx$id <- rep(id, length(unique(dx$pol))) + +nmhc_ldv_g_exh <- emis_chem2( + df = dx, + mech = mech, + nx = nx, + na.rm = TRUE +) + +# 4) LDV Gasoline Evaporative #### +x <- dfnmhc[family == "LDV" & + fuel == "G" & + process2 == "evap", + c("id", ntfs), with = F] + +x[is.na(x)] <- 0 +id <- x$id +x$id <- NULL + +dx <- speciate( + x = x, + spec = "nmhc", + fuel = "G", + veh = "LDV", + eu = "Evaporative" +) + +dx$id <- rep(id, length(unique(dx$pol))) + +nmhc_ldv_g_evap <- emis_chem2( + df = dx, + mech = mech, + nx = nx, + na.rm = TRUE +) + +# 5) LDV E85 Exhaust #### +x <- dfnmhc[family == "LDV" & + fuel == "E85" & + process2 == "exhaust", + c("id", ntfs), with = F] + +x[is.na(x)] <- 0 +id <- x$id +x$id <- NULL + +dx <- speciate( + x = x, + spec = "nmhc", + fuel = "E85", + veh = "LDV", + eu = "Exhaust" +) + +dx$id <- rep(id, length(unique(dx$pol))) + +nmhc_ldv_e85_exh <- emis_chem2( + df = dx, + mech = mech, + nx = nx, + na.rm = TRUE +) + + +# 6) LDV E85 Evaporative #### +x <- dfnmhc[family == "LDV" & + fuel == "E85" & + process2 == "evap", + c("id", ntfs), with = F] + +x[is.na(x)] <- 0 +id <- x$id +x$id <- NULL + +dx <- speciate( + x = x, + spec = "nmhc", + fuel = "E85", + veh = "LDV", + eu = "Evaporative" +) + +dx$id <- rep(id, length(unique(dx$pol))) + +nmhc_ldv_e85_evap <- emis_chem2( + df = dx, + mech = mech, + nx = nx, + na.rm = TRUE +) + + +# 7) CNG #### +x <- dfnmhc[fuel == "CNG", + c("id", ntfs), with = F] + +x[is.na(x)] <- 0 +id <- x$id +x$id <- NULL + +dx <- speciate( + x = x, + spec = "nmhc", + fuel = "E85", + veh = "LDV", + eu = "Evaporative" +) + +dx$id <- rep(id, length(unique(dx$pol))) + +nmhc_cng <- emis_chem2( + df = dx, + mech = mech, + nx = nx, + na.rm = TRUE +) + +# aggregate VOC #### +voc <- rbind(nmhc_ldv_d_exh, + nmhc_hdv_d_exh, + nmhc_ldv_g_exh, + nmhc_ldv_g_evap, + nmhc_ldv_e85_exh, + nmhc_ldv_e85_evap, + nmhc_cng, + use.names = FALSE +) + +dfvoc <- voc[, + lapply(.SD, sum, na.rm = T), + .SDcols = nx, + by = .(id, group) +] + +voc <- split(dfvoc, dfvoc$group) + +names(voc) <- paste0("E_", toupper(names(voc))) + +# saving VOC #### +voc[is.na(voc)] <- 0 +names(voc) + +# +# +# +for (i in seq_along(names(voc))) { + + x <- st_sf(voc[[names(voc)[i]]], geometry= net$geom) + x$group <- NULL + saveRDS(x, file = paste0("post/", mech, "/STREETS_", names(voc)[i], ".rds")) +} + +# +# # other gases #### +# for (j in seq_along(pol)) { +# x <- st_set_geometry(readRDS(paste0("post/grids/", pol[j], ".rds")), NULL) +# x[is.na(x)] <- 0 +# mm_x <- units::set_units(mol[j], "g/mol") # mm: massa molar +# for (i in 2:ncol(x)) x[, i] <- x[, i] * (mm_x)^-1 +# x <- st_sf(x, geometry = g$geometry) +# saveRDS(x, paste0("post/", mech, "/E_", pol[j], ".rds")) +# } +# +# # PM +# gPM <- st_set_geometry(readRDS("post/grids/PM.rds"), NULL) +# names(gPM) +# gPM$id <- NULL +# gPM1 <- speciate(x = gPM, spec = aer, list = T) +# +# for (i in 1:length(names(gPM1))) { +# gPMx <- st_sf(gPM1[[i]], geometry = g$geometry) +# saveRDS(gPMx, paste0("post/", mech, "/", toupper(names(gPM1))[i], ".rds")) +# } +# for (i in 1:ncol(gPM)) gPM[, i] <- units::set_units(gPM[, i], "ug/m^2/s") +# +# gPM10 <- st_set_geometry(readRDS("post/grids/PM10.rds"), NULL) +# names(gPM10) +# gPM10$id <- NULL +# for (i in 1:ncol(gPM10)) gPM10[, i] <- units::set_units(gPM10[, i], "ug/m^2/s") +# +# +# # gPM2510 are in g/km/h +# # need to change units to ug/m2/s +# gPM2510 <- gPM10 - gPM +# gPM2510 <- st_sf(gPM2510, geometry = g$g) +# saveRDS(gPM2510, paste0("post/", mech, "/E_PM_10.rds")) + + +switch(language, + "portuguese" = message(paste0("\nEmissões (mol/h) em npost/", mech)), + "english" = message(paste0("\nEmissions (mol/h) in post/", mech)), + "spanish" = message(paste0("\nEmisiones (mol/h) en post/", mech)) +) +fi <- list.files(paste0("post/", mech), full.names = T) +print(fi) + +switch(language, + "portuguese" = message("\nLembra que tu tem emissões de poluentes criterio em post/streets"), + "english" = message("\nRemember that you have criteria emissions in post/streets"), + "spanish" = message("\nRecuerda que tu tienes mas emisiones de conaminantes criterio en post/streets") +) + +fi <- list.files("post/streets", pattern = "ID",full.names = T) +print(fi) + +ls() +suppressWarnings( + rm("g", "gPM", "gPM1", "gPM10", "gPM2510", "gPMx", "i", "pol") +) \ No newline at end of file diff --git a/projects/moves/scripts/net.R b/projects/moves/scripts/net.R new file mode 100755 index 00000000..948ca9d2 --- /dev/null +++ b/projects/moves/scripts/net.R @@ -0,0 +1,60 @@ +switch(language, + "portuguese" = cat("\nNomes: ", names(net), "\n"), + "english" = cat("\nNames: ", names(net), "\n"), + "spanish" = cat("\nNombres: ", names(net), "\n") +) + + +net <- st_transform(net, crs) +intersect(categories, names(net)) +net$sum_categories <- rowSums(st_set_geometry(net, NULL)[, categories], na.rm = T) + +switch(language, + "portuguese" = cat("Filtrando ruas sem fluxo: \n"), + "english" = cat("Filtering empty roads\n"), + "spanish" = cat("Filtrando calles sin volumenes\n") +) + +net <- net[net$sum_categories > 0, ] +# +for (i in seq_along(categories)) { + png( + filename = paste0("images/NET_", categories[i], ".png"), + width = 2300, height = 1500, units = "px", pointsize = 12, + bg = "white", res = 300 + ) + plot(net[categories[i]], + axes = T, + pal = cpt(colorRampPalette = T, rev = T), + main = paste0(tit, ": ", categories[i]) + ) + dev.off() +} +net$lkm <- sf::st_length(net) + +net$lkm <- units::set_units(x = net$lkm, "miles") + +saveRDS(net, "network/net.rds") + +switch(language, + "portuguese" = message("\nArquivos em: /net:"), + "english" = message("\nFiles in: /net"), + "spanish" = message("\nArchivos en: /net") +) + +switch(language, + "portuguese" = message("Figuras em: /images:"), + "english" = message("Figures in: /images"), + "spanish" = message("Figuras en: /images") +) + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + +suppressWarnings( + rm(i, tit, net, categories, crs) +) +invisible(gc()) diff --git a/projects/moves/scripts/plots.R b/projects/moves/scripts/plots.R new file mode 100755 index 00000000..aed59aa1 --- /dev/null +++ b/projects/moves/scripts/plots.R @@ -0,0 +1,320 @@ +met <- readRDS("config/met.rds") +names(met) +year +year <- met$Year[1] +month <- met$Month[1] +date <- as.Date(ISOdate(year, month, 1, 0, 0)) + +numberOfDays <- function(date) { + m <- format(date, format = "%m") + + while (format(date, format = "%m") == m) { + date <- date + 1 + } + return(as.integer(format(date - 1, format = "%d"))) +} + +factor_emi <- numberOfDays(date = date) # daily to month + +# +# check hours +if (length(hours) > nrow(tfs)) { + switch(language, + "portuguese" = stop("\nEscolhe um número de horas menor ou igual que numero de horas de tfs\n"), + "english" = stop("\nChoose a number of hours less than or equal to the number of hours of tfs\n"), + "spanish" = stop("\nEscoje un numero de horas menor o igual que el numero de horas de tfs\n") + ) +} + +hh <- paste0(ifelse(nchar(hours) < 2, paste0("0", hours), hours), "00") + +# streets #### +switch(language, + "portuguese" = cat("\nPlotando ruas\n"), + "english" = cat("\nPlotting streets\n"), + "spanish" = cat("\nPlotando calles\n") +) + +pol <- list.files(path = "post/streets/", full.names = T) + +npol <- list.files(path = "post/streets/", full.names = F) + +npol <- gsub(pattern = ".rds", replacement = "", x = npol) + +for (i in seq_along(pol)) { + + for (j in seq_along(hours)) { + + x <- readRDS(pol[i]) + + cn <- names(x)[hours + 1] + + png( + filename = paste0("images/STREETS_", npol[i], "_", hh[j], ".png"), + width = 2500, + height = 2500, + units = "px", + pointsize = 12, + bg = "white", + res = 300 + ) + par(bg = "white") + plot(x[as.numeric(x[[cn[j]]]) > 0, ][cn[j]], + axes = TRUE, + bg = bg, + main = paste0(npol[i], ": ", tit, " [g/h] ", hh[j], " LT"), + pal = cptcity::cpt(colorRampPalette = TRUE, rev = FALSE, pal = pal), lwd = 2 + ) + dev.off() + } +} + +# grade #### +switch(language, + "portuguese" = cat("\nPlotando grades\n"), + "english" = cat("\nPlotting grids\n"), + "spanish" = cat("\nPlotando grillas\n") +) + + +pol <- list.files(path = "post/grids/", full.names = T) + +npol <- list.files(path = "post/grids/", full.names = F) + +npol <- gsub(pattern = ".rds", replacement = "", x = npol) + +for (i in seq_along(pol)) { + + for (j in seq_along(hours)) { + + x <- readRDS(pol[i]) + + cn <- names(x)[hours + 1] + + png( + filename = paste0("images/GRIDS_", npol[i], "_", hh[j], ".png"), + width = 2500, + height = 2500, + units = "px", + pointsize = 12, + bg = "white", + res = 300 + ) + par(bg = "white") + plot(x[as.numeric(x[[cn[j]]]) > 0, ][cn[j]], + axes = TRUE, + bg = bg, + lty = 0.3, + main = paste0(npol[i], ": ", tit, " [g/km^2/h] ", hh[j], " LT"), + pal = cptcity::cpt(colorRampPalette = TRUE, rev = TRUE, pal = pal) + ) + dev.off() + } +} + + +# categoria #### +dt <- readRDS("post/datatable/emissions.rds") +dt0 <- dt[, round(sum(t) * factor_emi, 2), by = .(pollutant, process)] + +dt$veh <- as.character(dt$veh) + +# totais x veicles +switch(language, + "portuguese" = cat("\nPlotando categorias por total\n"), + "english" = cat("\nPlotting categories by total\n"), + "spanish" = cat("\nPlotando categorias por total\n") +) + +dt1 <- dt[, + as.numeric(sum(g)), + by = .(pollutant, veh) +] + +dt1$veh <- factor( + x = dt1$veh, + levels = metadata$vehicles +) + +pol <- unique(dt1$pollutant) + +for (i in seq_along(pol)) { + + p <- ggplot( + dt1[pollutant == pol[i]], + aes(x = veh, y = V1, fill = V1) + ) + + geom_bar(stat = "identity", col = "black") + + labs( + y = paste0("g/", nrow(tfs),"h"), + title = paste0(pol[i], ": ", tit) + ) + + scale_fill_gradientn(pol[i], colours = cpt()) + + scale_x_discrete(limits = rev(metadata$vehicles)) + + theme_bw() + + coord_flip() + + png( + filename = paste0("images/TOTAL_", gsub("\\.","_",gsub(" ", "", pol[i])), ".png"), + width = 2500, + height = 2500, + units = "px", + pointsize = 12, + bg = "white", + res = 300 + ) + print(p) + dev.off() +} + +# totals by vehicles +switch(language, + "portuguese" = cat("\nPlotando categorias por type_emi\n"), + "english" = cat("\nPlotting categories by type_emi\n"), + "spanish" = cat("\nPlotando categorias por type_emi\n") +) + +dt1 <- dt[, + as.numeric(sum(g)), + by = .(pollutant, veh, process) +] + +dt1$veh <- factor( + x = dt1$veh, + levels = metadata$vehicles +) + +dt1$pol_te <- as.character(paste0(dt1$pollutant, " - ", dt1$process)) + +pol <- as.character(unique(dt1$pollutant)) + +for (i in seq_along(pol)) { + + p <- ggplot( + dt1[pollutant == pol[i]], + aes(x = veh, + y = V1, + fill = process) + ) + + geom_bar(stat = "identity", col = "black") + + labs( + y = paste0("g/", nrow(tfs),"h"), + title = paste0(pol[i], ": ", tit) + ) + + scale_x_discrete(limits = rev(metadata$vehicles)) + + theme_bw() + + coord_flip() + + png( + filename = paste0("images/TOTAL_", gsub("\\.","_",gsub(" ", "", pol[i])), "_FILL.png"), + width = 2500, + height = 2500, + units = "px", + pointsize = 12, + bg = "white", + res = 300 + ) + print(p) + dev.off() +} + +# totals by hour +switch(language, + "portuguese" = cat("\nPlotando categorias por hora\n"), + "english" = cat("\nPlotting categories by hour\n"), + "spanish" = cat("\nPlotando categorias por hora\n") +) + +dt1 <- dt[, + as.numeric(sum(g)), + by = .(pollutant, family, hour) +] + +for (i in seq_along(pol)) { + + p <- ggplot( + dt1[pollutant == pol[i]], + aes(x = hour, + y = V1, + fill = family) + ) + + geom_bar(stat = "identity", col = "black") + + labs( + y = paste0("g/", nrow(tfs),"h"), + title = paste0(pol[i], ": ", tit) + ) + + theme_bw() + + png( + filename = paste0("images/HOUR_", gsub("\\.","_",gsub(" ", "", pol[i])), "_FILL.png"), + width = 2500, + height = 2500, + units = "px", + pointsize = 12, + bg = "white", + res = 300 + ) + print(p) + dev.off() +} + + +# totais x age +dt1 <- dt[, + as.numeric(sum(g)), + by = .(pollutant, family, age) +] + +for (i in seq_along(pol)) { + p <- ggplot( + dt1[pollutant == pol[i]], + aes(x = age, + y = V1, + fill = family) + ) + + geom_bar(stat = "identity", col = "black") + + labs( + y = paste0("g/", nrow(tfs),"h"), + title = paste0(pol[i], ": ", tit) + ) + + theme_bw() + + png( + filename = paste0("images/AGE_", gsub("\\.","_",gsub(" ", "", pol[i])), "_FILL.png"), + width = 2500, + height = 2500, + units = "px", + pointsize = 12, + bg = "white", + res = 300 + ) + print(p) + dev.off() +} + +switch(language, + "portuguese" = message("\nFiguras em /images\n"), + "english" = message("\nFigures in /image\n"), + "spanish" = message("\nFiguras en /images\n") +) + + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + + +suppressWarnings( + rm( + "bg", "breaks", "cn", "dt", "dt0", "dt1", "ef_cetesb2", "factor_emi", + "hh", "hours", "i", "j", "metadata", "n_BUS", "n_LCV", "n_MC", "n_PC", + "n_TRUCKS", "num_vein", "p", "pal", "pol", "pole", "tfs", "tit", "uv", + "veh", "vein_version", "x", "year" + ) +) + + +ls() +invisible(gc()) \ No newline at end of file diff --git a/projects/moves/scripts/post.R b/projects/moves/scripts/post.R new file mode 100755 index 00000000..fc6dc16b --- /dev/null +++ b/projects/moves/scripts/post.R @@ -0,0 +1,236 @@ +met <- readRDS("config/met.rds") +names(met) + +year <- met$Year[1] +month <- met$Month[1] +date <- as.Date(ISOdate(year, month, 1, 0, 0)) + +numberOfDays <- function(date) { + m <- format(date, format = "%m") + + while (format(date, format = "%m") == m) { + date <- date + 1 + } + return(as.integer(format(date - 1, format = "%d"))) +} + +factor_emi <- numberOfDays(date = date) # daily to month + +# name_hours #### +ntfs <- paste0("H", 1:nrow(tfs)) + +# electric #### +f <- list.files(path = "emi", pattern = "ELEC", full.names = T, recursive = T) +f <- grep(pattern = "STREETS", x = f, value = T) +df_elec <- rbindlist(lapply(f, fread)) +names(df_elec) + +df_elec[, + lapply(.SD, sum, na.rm = T), + .SDcols = ntfs, + by = .(id, pollutantID)] -> emi_elec + +lemi_elec <- split(emi_elec, emi_elec$pollutantID) +names(lemi_elec) + +# pollutants #### +switch(language, + "portuguese" = message("\nDetectando poluentes:\n"), + "english" = message("\nDetecting pollutants:\n"), + "spanish" = message("\nDetectando contaminantes\n") +) + +f <- list.files(path = "emi", pattern = "STREETS", full.names = T, recursive = T) + +ff <- gsub(pattern = "emi/emissions_rate_per_distance/STREETS_pollutantID_", + replacement = "", + f) +ff <- gsub(pattern = "emi/emissions_rate_per_vehicle/STREETS_pollutantID_", + replacement = "", + ff) +ff <- gsub(pattern = ".csv.gz", + replacement = "", + ff) +fff <- unique(suppressWarnings(as.numeric(ff))) +fff <- fff[!is.na(fff)] +fff <- fff[order(fff)] + +dfpol <- decoder$pollutants +setDT(dfpol) + +print(dfpol[ID %in% fff]) + + +# streets #### +switch(language, + "portuguese" = message("\nAgregando emissões por rua...\n"), + "english" = message("\nAgregating emissions by street...\n"), + "spanish" = message("\nAgregando emisiones por calle...\n") +) +f <- list.files(path = "emi", pattern = "STREETS", full.names = T, recursive = T) + +# emissions differnet than break and tire wear +for(i in seq_along(fff)) { + if(fff[i] %in% c(106, 107, 116, 117)) next + cat("Reading pollutantID", fff[i], "\n") + df <- rbindlist(lapply(grep(pattern = paste0("pollutantID_", fff[i]), x = f, value = T), fread)) + df[, + lapply(.SD, sum, na.rm = T), + .SDcols = ntfs, + by = .(id)] -> dfx + + for(j in seq_along(ntfs)) dfx[[ntfs[j]]] <- Emissions(dfx[[ntfs[j]]], mass = "g", time = "h") + + dfx <- st_sf(dfx, geometry = st_geometry(net)) + saveRDS(object = dfx, paste0("post/streets/pollutantID_", fff[i], ".rds")) +} + +# break and tire wear +nl <-names(lemi_elec) +for(i in seq_along(nl)) { + dx <- lemi_elec[[nl[i]]] + dx[, pollutantID := NULL] + df <- rbindlist(lapply(grep(pattern = paste0("pollutantID_", nl[i]), x = f, value = T), fread)) + + df <- rbind(dx, df) + df[, + lapply(.SD, sum, na.rm = T), + .SDcols = ntfs, + by = .(id)] -> dfx + + for(j in seq_along(ntfs)) dfx[[ntfs[j]]] <- Emissions(dfx[[ntfs[j]]], mass = "g", time = "h") + + dfx <- st_sf(dfx, geometry = st_geometry(net)) + + saveRDS(object = dfx, paste0("post/streets/pollutantID_", nl[i], ".rds")) +} + +# nmhc +nmhcs <- grep(pattern = "NMHC", x = f, value = T) +nmhc1 <- fread(nmhcs[1]) +nmhc2 <- fread(nmhcs[2], fill = T) # problem with emissions_per_distance_NMHC + +df <- rbind(nmhc1, nmhc2) + +df[, + lapply(.SD, sum, na.rm = T), + .SDcols = ntfs, + by = .(id, family, fuel, process2)] -> dfx + +for(j in seq_along(ntfs)) dfx[[ntfs[j]]] <- Emissions(dfx[[ntfs[j]]], mass = "g", time = "h") + +dfx <- st_sf(dfx, geometry = st_geometry(net)) + +saveRDS(object = dfx, "post/streets/NMHC.rds") + +# grids #### +switch(language, + "portuguese" = message("\nAgregando emissões por grade\n"), + "english" = message("\nAgregating emissions by grid...\n"), + "spanish" = message("\nAgregando emisiones por grilla...\n") +) + +g <- make_grid(spobj = dfx, width = grid_width, crs = crs) + +lf <- list.files(path = "post/streets", pattern = ".rds", full.names = TRUE) +na <- list.files(path = "post/streets", pattern = ".rds", full.names = F) +na <- gsub(pattern = ".rds", replacement = "", x = na) + +for (i in seq_along(lf)) { + switch(language, + "portuguese" = message("\nGerando grade de ", lf[i], "\n"), + "english" = message("\nGridding emissions of ", lf[i], "\n"), + "spanish" = message("\nGenerando grilla de ", lf[i], "\n") + ) + x <- readRDS(lf[i]) + gx <- emis_grid(spobj = x, g = g, sr = crs) + saveRDS(gx, paste0("post/grids/", na[i], ".rds")) +} + +# datatable #### +switch(language, + "portuguese" = message("\nAgregando emissões em data.table\n"), + "english" = message("\nAgregating emissions in data.table...\n"), + "spanish" = message("\nAgregando emisiones en data.table...\n") +) +f <- list.files(path = "emi", pattern = "DF", full.names = T, recursive = T) +f1 <- fread(f[1]) +f2 <- fread(f[2]) +f3 <- fread(f[3]) +f2$fuelTypeID <- 9 +f2$fuel <- "ELEC" +names(f1) +names(f2) +names(f3) + +dt <- rbind(f1[,names(f3), with = F], + f2[,names(f3), with = F], + f3) +names(dt) +dt$g <- units::set_units(dt$value, g) +dt$t <- units::set_units(dt$g, t) +data("decoder") + +# pollutant +dfpol <- decoder$pollutants +names(dfpol)[1] <- "pollutantID" +dt <- merge(dt, dfpol, by = "pollutantID", all.x = T) + +# process +dfpro <- decoder$emission_process +names(dfpro) <- c("processID", "process") +dt <- merge(dt, dfpro, by = "processID", all.x = T) + +# vehicle +dfv <- metadata[, c("family", "vehicles")] +names(dfv)[2] <- "veh" +dt <- merge(dt, dfv, by = "veh", all.x = T) + +saveRDS(dt, "post/datatable/emissions.rds") + +dt0 <- dt[, round(sum(t), 2), by = .(pollutant)] +print(dt0) + +# emissoes by veh +dt1 <- dt[, sum(t), by = .(pollutant, veh)] +df1 <- long_to_wide(df = dt1, + column_with_new_names = "pollutant", + column_with_data = "V1", + column_fixed = "veh") +saveRDS(df1, "post/datatable/emissions_by_veh.rds") +data.table::fwrite(df1, "csv/emissions_by_veh.csv", row.names = FALSE) + +# emissoes by fuel +dt2 <- dt[, sum(t), by = .(pollutant, fuel)] +df2 <- long_to_wide(df = dt2, + column_with_new_names = "pollutant", + column_with_data = "V1", + column_fixed = "fuel") +saveRDS(df2, "post/datatable/emissions_by_fuel.rds") +data.table::fwrite(df2, "csv/emissions_by_fuel.csv", row.names = FALSE) + +# emissoes by age +dt3 <- dt[, sum(t), by = .(pollutant, age)] +df3 <- long_to_wide(df = dt3, column_with_new_names = "pollutant", column_with_data = "V1", column_fixed = "age") +saveRDS(df3, "post/datatable/emissions_by_age.rds") +data.table::fwrite(df3, "csv/emissions_by_age.csv", row.names = FALSE) + +switch(language, + "portuguese" = message("\n\nArquivos em: /post/*:"), + "english" = message("\nFiles in: /post/*"), + "spanish" = message("\nArchivos en: /post/*") +) + + +switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") +) + +suppressWarnings( + rm( + "df1", "df2", "df3", "dt", "dt0", "dt1", "dt2", "dt3", "factor_emi", + "g", "gx", "i", "lf", "na", "net", "pol", "pols", "x", "crs" + ) +) \ No newline at end of file diff --git a/projects/moves/scripts/traffic.R b/projects/moves/scripts/traffic.R new file mode 100755 index 00000000..975257ce --- /dev/null +++ b/projects/moves/scripts/traffic.R @@ -0,0 +1,661 @@ + + +# n_veh + +n_MC <- metadata$vehicles[grep(pattern = "MC", x = metadata$vehicles)] +n_PC <- metadata$vehicles[grep(pattern = "PC", x = metadata$vehicles)] +n_PT <- metadata$vehicles[grep(pattern = "PT", x = metadata$vehicles)] +n_LCT <- metadata$vehicles[grep(pattern = "LCT", x = metadata$vehicles)] +n_BUS_INTERCITY <- metadata$vehicles[grep(pattern = "BUS_INTERCITY", x = metadata$vehicles)] +n_BUS_TRANSIT <- metadata$vehicles[grep(pattern = "BUS_TRANSIT", x = metadata$vehicles)] +n_BUS_SCHOOL <- metadata$vehicles[grep(pattern = "BUS_SCHOOL", x = metadata$vehicles)] +n_TRUCKS_REFUSE <- metadata$vehicles[grep(pattern = "TRUCKS_REFUSE", x = metadata$vehicles)] +n_TRUCKS_SU_SH <- metadata$vehicles[grep(pattern = "TRUCKS_SU_SH", x = metadata$vehicles)] +n_TRUCKS_SU_LH <- metadata$vehicles[grep(pattern = "TRUCKS_SU_LH", x = metadata$vehicles)] +n_TRUCKS_MH <- metadata$vehicles[grep(pattern = "TRUCKS_MH", x = metadata$vehicles)] +n_TRUCKS_CU_SH <- metadata$vehicles[grep(pattern = "TRUCKS_CU_SH", x = metadata$vehicles)] +n_TRUCKS_CU_LH <- metadata$vehicles[grep(pattern = "TRUCKS_CU_LH", x = metadata$vehicles)] + +n_veh <- list( + MC = n_MC, + PC = n_PC, + PT = n_PT, + LCT = n_LCT, + BUS_INTERCITY = n_BUS_INTERCITY, + BUS_TRANSIT= n_BUS_TRANSIT, + BUS_SCHOOL = n_BUS_SCHOOL, + TRUCKS_REFUSE = n_TRUCKS_REFUSE, + TRUCKS_SU_SH = n_TRUCKS_SU_SH, + TRUCKS_SU_LH = n_TRUCKS_SU_LH, + TRUCKS_MH = n_TRUCKS_MH, + TRUCKS_CU_SH = n_TRUCKS_CU_SH, + TRUCKS_CU_LH = n_TRUCKS_CU_LH +) + +for(i in seq_along(metadata$vehicles)) { + if(identical(names(veh)[i + 1],metadata$vehicles[i]) == FALSE) { + + switch (language, + "portuguese" = cat( "Nomes incompativeis:\n", + "metadata:",metadata$vehicles[i], "\n", + "`veh`:",names(veh)[i + 1], "\n"), + "english" = cat( "Incomptaible names:\n", + "metadata:",metadata$vehicles[i], "\n", + "`veh`:",names(veh)[i + 1], "\n"), + "spanish" = cat( "Nombres incomptatibles:\n", + "metadata:",metadata$vehicles[i], "\n", + "`veh`:",names(veh)[i + 1], "\n")) + + + stop() + } +} + + +# apagando arquivos +switch (language, + "portuguese" = message("Apagando veh/*.rds\n"), + "english" = message("Deleting veh/*.rds\n"), + "spanish" = message("Borrando veh/*.rds\n")) + +arquivos <- list.files(path = "veh", pattern = ".rds", full.names = TRUE) +file.remove(arquivos) + +# fleet age +veh[is.na(veh)] <- 0 + +veh <- remove_units(veh) + +# plotting +switch (language, + "portuguese" = cat("Plotando fluxos\n"), + "english" = cat("Plotting traffic flows\n"), + "spanish" = cat("Plotando flujos\n")) + +# list of vehicles to colplot +lveh <- list() + + +# MC #### +fl <- function(x) { + lapply(seq_along(x), function(i) { + if(length(x) == 1) { + 1 + } else { + df <- veh[, x] + sum(df[i]) / sum(df, na.rm = T) + } + })-> lx + names(lx) <- x + lx +} +fl(n_MC) -> l_mc + +lapply(seq_along(l_mc), function(i) { + cat(names(l_mc)[i], " is ", n_MC[[i]], " of MC\n") +}) -> lx +cat("sum: ", sum(unlist(l_mc)), "\n\n") + +# fuel +metadata[metadata$vein_name == "MC", ]$fuel +kf <- c(k_G) #gasoline + +for (i in seq_along(n_MC)) { + x <- my_age( + x = net[[categories[1]]], + y = veh[[n_MC[i]]], + name = n_MC[i], + k = l_mc[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_MC[i], ".rds")) +} +lMC <- data.frame(MC_G = colSums(x)) + +# PC #### +fl(n_PC) -> l_pc + +lapply(seq_along(l_pc), function(i) { + cat(names(l_pc)[i], " is ", l_pc[[i]], " of PC\n") +}) -> lx +cat("sum: ", sum(unlist(l_pc)), "\n\n") + + +# fuel +metadata[metadata$vein_name == "PC", ]$fuel +kf <- c(k_G, #diesel + 1, #eletric + k_E85, #E85 + k_D) #gasoline + +# colplot +lPC <- list() + +for (i in seq_along(n_PC)) { + x <- my_age( + x = net[[categories[2]]], + y = veh[[n_PC[i]]], + name = n_PC[i], + k = l_pc[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_PC[i], ".rds")) + lPC[[i]] <- colSums(x) +} + +# PT #### +fl(n_PT) -> l_pt + +lapply(seq_along(l_pt), function(i) { + cat(names(l_pt)[i], " is ", l_pt[[i]], " of PT\n") +}) -> lx +cat("sum: ", sum(unlist(l_pt)), "\n\n") + +# fuel +metadata[metadata$family == "PT", ]$fuel +kf <- c(k_G, #diesel + 1, #eletric + k_E85, #E85 + k_D) #gasoline + +# colplot +lPT <- list() + +for (i in seq_along(n_PT)) { + x <- my_age( + x = net[[categories[3]]], + y = veh[[n_PT[i]]], + name = n_PT[i], + k = l_pt[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_PT[i], ".rds")) + lPT[[i]] <- colSums(x) +} + +# LCT #### +fl(n_LCT) -> l_lct + +lapply(seq_along(l_lct), function(i) { + cat(names(l_lct)[i], " is ", l_lct[[i]], " of LCT\n") +}) -> lx +cat("sum: ", sum(unlist(l_lct)), "\n\n") + +# fuel +metadata[metadata$vein_name == "LCT", ]$fuel +kf <- c(k_G, #diesel + 1, #eletric + k_E85, #E85 + k_D) #gasoline + +# colplot +lLCT <- list() + +for (i in seq_along(n_LCT)) { + x <- my_age( + x = net[[categories[4]]], + y = veh[[n_LCT[i]]], + name = n_LCT[i], + k = l_lct[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_LCT[i], ".rds")) + lLCT[[i]] <- colSums(x) +} + +# BUS_INTERCITY #### +fl(n_BUS_INTERCITY) -> l_bus_intercity + +lapply(seq_along(l_bus_intercity), function(i) { + cat(names(l_bus_intercity)[i], " is ", l_bus_intercity[[i]], " of BUS_INTERCITY\n") +}) -> lx +cat("sum: ", sum(unlist(l_bus_intercity)), "\n\n") + +# fuel +metadata[metadata$vein_name == "BUS_INTERCITY", ]$fuel +kf <- c(k_G, + k_D, + k_CNG) + +# colplot +lBI <- list() + +for (i in seq_along(n_BUS_INTERCITY)) { + x <- my_age( + x = net[[categories[5]]], + y = veh[[n_BUS_INTERCITY[i]]], + name = n_BUS_INTERCITY[i], + k = l_bus_intercity[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_BUS_INTERCITY[i], ".rds")) + lBI[[i]] <- colSums(x) +} + + +# BUS_TRANSIT #### +fl(n_BUS_TRANSIT) -> l_bus_transit + +lapply(seq_along(l_bus_transit), function(i) { + cat(names(l_bus_transit)[i], " is ", l_bus_transit[[i]], " of BUS_TRANSIT\n") +}) -> lx +cat("sum: ", sum(unlist(l_bus_transit)), "\n\n") + +# fuel +metadata[metadata$vein_name == "BUS_TRANSIT", ]$fuel +kf <- c(k_G, + k_D, + k_CNG) + +# colplot +lBT <- list() + +for (i in seq_along(n_BUS_TRANSIT)) { + x <- my_age( + x = net[[categories[6]]], + y = veh[[n_BUS_TRANSIT[i]]], + name = n_BUS_TRANSIT[i], + k = l_bus_transit[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_BUS_TRANSIT[i], ".rds")) + lBT[[i]] <- colSums(x) +} + + +# BUS_SCHOOL #### +fl(n_BUS_SCHOOL) -> l_bus_school + +lapply(seq_along(l_bus_school), function(i) { + cat(names(l_bus_school)[i], " is ", l_bus_school[[i]], " of BUS_SCHOOL\n") +}) -> lx +cat("sum: ", sum(unlist(l_bus_school)), "\n\n") + +# fuel +metadata[metadata$vein_name == "BUS_SCHOOL", ]$fuel +kf <- c(k_G, + k_D, + k_CNG) + +# colplot +lBS <- list() + +for (i in seq_along(n_BUS_SCHOOL)) { + x <- my_age( + x = net[[categories[7]]], + y = veh[[n_BUS_SCHOOL[i]]], + name = n_BUS_SCHOOL[i], + k = l_bus_school[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_BUS_SCHOOL[i], ".rds")) + lBS[[i]] <- colSums(x) +} + + +# TRUCKS_REFUSE #### +fl(n_TRUCKS_REFUSE) -> l_trucks_refuse + +lapply(seq_along(l_trucks_refuse), function(i) { + cat(names(l_trucks_refuse)[i], " is ", l_trucks_refuse[[i]], " of TRUCKS_REFUSE\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_refuse)), "\n\n") + +# fuel +metadata[metadata$vein_name == "TRUCKS_REFUSE", ]$fuel +kf <- c(k_G, + k_D, + k_CNG) + +# colplot +lTR <- list() + +for (i in seq_along(n_TRUCKS_REFUSE)) { + x <- my_age( + x = net[[categories[8]]], + y = veh[[n_TRUCKS_REFUSE[i]]], + name = n_TRUCKS_REFUSE[i], + k = l_trucks_refuse[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_REFUSE[i], ".rds")) + lTR[[i]] <- colSums(x) +} + + +# TRUCKS_SU_SH #### +fl(n_TRUCKS_SU_SH) -> l_trucks_su_sh + +lapply(seq_along(l_trucks_su_sh), function(i) { + cat(names(l_trucks_su_sh)[i], " is ", l_trucks_su_sh[[i]], " of TRUCKS_SU_SH\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_su_sh)), "\n\n") + +# fuel +metadata[metadata$vein_name == "TRUCKS_SU_SH", ]$fuel +kf <- c(k_G, + k_D, + k_CNG) + +# colplot +lTSS <- list() + +for (i in seq_along(n_TRUCKS_SU_SH)) { + x <- my_age( + x = net[[categories[9]]], + y = veh[[n_TRUCKS_SU_SH[i]]], + name = n_TRUCKS_SU_SH[i], + k = l_trucks_su_sh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_SU_SH[i], ".rds")) + lTSS[[i]] <- colSums(x) +} + +# TRUCKS_SU_LH #### +fl(n_TRUCKS_SU_LH) -> l_trucks_su_lh + +lapply(seq_along(l_trucks_su_lh), function(i) { + cat(names(l_trucks_su_lh)[i], " is ", l_trucks_su_lh[[i]], " of TRUCKS_SU_LH\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_su_lh)), "\n\n") + +# fuel +metadata[metadata$vein_name == "TRUCKS_SU_LH", ]$fuel +kf <- c(k_G, + k_D, + k_CNG) + +# colplot +lTSL <- list() + +for (i in seq_along(n_TRUCKS_SU_LH)) { + x <- my_age( + x = net[[categories[10]]], + y = veh[[n_TRUCKS_SU_LH[i]]], + name = n_TRUCKS_SU_LH[i], + k = l_trucks_su_lh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_SU_LH[i], ".rds")) + lTSL[[i]] <- colSums(x) +} + + +# TRUCKS_MH #### +fl(n_TRUCKS_MH) -> l_trucks_mh + +lapply(seq_along(l_trucks_mh), function(i) { + cat(names(l_trucks_mh)[i], " is ", l_trucks_mh[[i]], " of TRUCKS_MH\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_mh)), "\n\n") + +# fuel +metadata[metadata$vein_name == "TRUCKS_MH", ]$fuel +kf <- c(k_G, + k_D, + k_CNG) + +# colplot +lMH <- list() + +for (i in seq_along(n_TRUCKS_MH)) { + x <- my_age( + x = net[[categories[11]]], + y = veh[[n_TRUCKS_MH[i]]], + name = n_TRUCKS_MH[i], + k = l_trucks_mh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_MH[i], ".rds")) + lMH[[i]] <- colSums(x) +} + + +# TRUCKS_CU_SH #### +fl(n_TRUCKS_CU_SH) -> l_trucks_cu_sh + +lapply(seq_along(l_trucks_cu_sh), function(i) { + cat(names(l_trucks_cu_sh)[i], " is ", l_trucks_cu_sh[[i]], " of TRUCKS_CU_SH\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_cu_sh)), "\n\n") + +# fuel +metadata[metadata$vein_name == "TRUCKS_CU_SH", ]$fuel +kf <- c(k_G, + k_D, + k_CNG) + +# colplot +lTCS <- list() + +for (i in seq_along(n_TRUCKS_CU_SH)) { + x <- my_age( + x = net[[categories[12]]], + y = veh[[n_TRUCKS_CU_SH[i]]], + name = n_TRUCKS_CU_SH[i], + k = l_trucks_cu_sh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_CU_SH[i], ".rds")) + lTCS[[i]] <- colSums(x) +} + + +# TRUCKS_CU_LH #### +fl(n_TRUCKS_CU_LH) -> l_trucks_cu_lh + +lapply(seq_along(l_trucks_cu_lh), function(i) { + cat(names(l_trucks_cu_lh)[i], " is ", l_trucks_cu_lh[[i]], " of TRUCKS_CU_LH\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_cu_lh)), "\n\n") + +# fuel +metadata[metadata$vein_name == "TRUCKS_CU_LH", ]$fuel +kf <- c(k_D) + +# colplot +lTCL <- list() + +for (i in seq_along(n_TRUCKS_CU_LH)) { + x <- my_age( + x = net[[categories[13]]], + y = veh[[n_TRUCKS_CU_LH[i]]], + name = n_TRUCKS_CU_LH[i], + k = l_trucks_cu_lh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_CU_LH[i], ".rds")) + lTCL[[i]] <- colSums(x) +} + + +# colplot #### + +DF <- cbind(lMC, + do.call("cbind", lPC), + do.call("cbind", lPT), + do.call("cbind", lLCT), + do.call("cbind", lBI), + do.call("cbind", lBT), + do.call("cbind", lBS), + do.call("cbind", lTR), + do.call("cbind", lTSS), + do.call("cbind", lTSL), + do.call("cbind", lMH), + do.call("cbind", lTCS), + do.call("cbind", lTCL)) + +names(DF) <- metadata$vehicles + +for(i in seq_along(n_veh)) { + png(paste0("images/FLOW_", n_veh[[i]], ".png"), 2000, 1500, "px", res = 300) + + colplot( + df = DF, + cols = n_veh[[i]], + xlab = "Age", + ylab = "veh/h", + main = paste0("FLOW ", n_veh[[i]]), + type = "l", + pch = NULL, + theme = theme, + spl = 9 + ) + dev.off() +} + + + +# plots #### +switch(language, + "portuguese" = cat("Plotando composição veicular\n"), + "english" = cat("Plotting vehicular composition\n"), + "spanish" = cat("Plotando composición vehicular\n") +) + + +f <- list.files(path = "veh", pattern = ".rds", full.names = T) +na <- list.files(path = "veh", pattern = ".rds") +na <- gsub(pattern = ".rds", replacement = "", x = na) +ff <- unlist(lapply(seq_along(f), function(i) { + sum(readRDS(f[i])) +})) +df <- data.frame( + per = ff, + Categoria = na +) +df <- df[df$Categoria != "fleet_age", ] +df$per <- df$per / sum(df$per) * 100 +df$Categoria <- factor(df$Categoria, levels = metadata$vehicles) +p <- ggplot(df, aes(x = Categoria, y = per, fill = per)) + + geom_bar(stat = "identity", col = "black") + + labs(y = "[%]", title = "Vehicular composition") + + ggplot2::scale_fill_gradientn("[%]", colours = cptcity::cpt(rev = T)) + + theme_bw() + + coord_flip() + + scale_x_discrete(limits = rev(metadata$vehicles)) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + +png( + filename = paste0("images/TRAFFIC.png"), + width = 2500, height = 2500, units = "px", pointsize = 12, + bg = "white", res = 300 +) +print(p) +dev.off() + + +if(speed) { + + + rbindlist(lapply(seq_along(metadata$vehicles), function(i) { + temp_fact(q = rowSums(readRDS(paste0("veh/", metadata$vehicles[i], ".rds"))), + pro = tfs[[metadata$vehicles[i]]]) -> tf + tf$id <- 1:nrow(tf) + tf + })) -> lx + + lx[, + lapply(.SD, sum, na.rm = T), + .SDcols = paste0("V", 1:length(tfs[[1]])), + by = id] -> ft + + ft$id <- NULL + + switch(language, + "portuguese" = cat("Projetando velocidade\n"), + "english" = cat("Projecting speed\n"), + "spanish" = cat("Proyectando velocidad\n") + ) + + df <- netspeed( + q = ft, + ps = net$CONGSPD, + ffs = net$FFSPEED, + cap = net$CAPE, + lkm = net$lkm, + alpha = 1, + dist = "miles" + ) + saveRDS(df, "network/speed.rds") + + dfs <- netspeed( + q = ft, + ps = net[[ps]], + ffs = net[[ffs]], + cap = net[[capacity]], + lkm = net$lkm, + alpha = 1, + dist = "miles", + net = net + ) + + dfs <- st_transform(dfs, 4326) + saveRDS(dfs, "network/speed_net.rds") + + + sp_mov <- moves_speed(df) + saveRDS(sp_mov, "network/speed_bin.rds") + + + # plots #### + switch(language, + "portuguese" = cat("Plotando velocidade\n"), + "english" = cat("Plotting speed\n"), + "spanish" = cat("Plotando velocidad\n") + ) + + png( + filename = paste0("images/SPEED.png"), + width = 2500, height = 3500, units = "px", pointsize = 12, + bg = "white", res = 300 + ) + plot(df, mai2 =c(2.2, 0.82, 0.82, 0.42)) + dev.off() + + png( + filename = paste0("images/SPEED_HOUR.png"), + width = 2500, height = 3000, units = "px", pointsize = 12, + bg = "white", res = 300 + ) + plot(dfs, + max.plot = 25, + axes = T, + pal = cpt(colorRampPalette = T, rev = T), + key.pos = 4) + + dev.off() + + + png( + filename = paste0("images/SPEEDBIN.png"), + width = 2500, height = 3500, units = "px", pointsize = 12, + bg = "white", res = 300 + ) + plot(Speed(remove_units(sp_mov), dist = units::unitless), pal = "mpl_viridis") + dev.off() + + +} + + switch(language, + "portuguese" = message("\nArquivos em: /veh:"), + "english" = message("\nFiles in: /veh"), + "spanish" = message("\nArchivos en: /veh") + ) + + switch(language, + "portuguese" = message("Figuras em: /images:"), + "english" = message("Figures in: /images"), + "spanish" = message("Figuras en: /images") + ) + + switch(language, + "portuguese" = message("Limpando..."), + "english" = message("Cleaning..."), + "spanish" = message("Limpiando...") + ) + + + \ No newline at end of file diff --git a/projects/moves/scripts/trafficfuel.R b/projects/moves/scripts/trafficfuel.R new file mode 100755 index 00000000..b115ee84 --- /dev/null +++ b/projects/moves/scripts/trafficfuel.R @@ -0,0 +1,632 @@ +k_E85 <- 1 + +# n_veh + +n_MC <- metadata$vehicles[grep(pattern = "MC", x = metadata$vehicles)] +n_PC <- metadata$vehicles[grep(pattern = "PC", x = metadata$vehicles)] +n_PT <- metadata$vehicles[grep(pattern = "PT", x = metadata$vehicles)] +n_LCT <- metadata$vehicles[grep(pattern = "LCT", x = metadata$vehicles)] +n_BUS_INTERCITY <- metadata$vehicles[grep(pattern = "BUS_INTERCITY", x = metadata$vehicles)] +n_BUS_TRANSIT <- metadata$vehicles[grep(pattern = "BUS_TRANSIT", x = metadata$vehicles)] +n_BUS_SCHOOL <- metadata$vehicles[grep(pattern = "BUS_SCHOOL", x = metadata$vehicles)] +n_TRUCKS_REFUSE <- metadata$vehicles[grep(pattern = "TRUCKS_REFUSE", x = metadata$vehicles)] +n_TRUCKS_SU_SH <- metadata$vehicles[grep(pattern = "TRUCKS_SU_SH", x = metadata$vehicles)] +n_TRUCKS_SU_LH <- metadata$vehicles[grep(pattern = "TRUCKS_SU_LH", x = metadata$vehicles)] +n_TRUCKS_MH <- metadata$vehicles[grep(pattern = "TRUCKS_MH", x = metadata$vehicles)] +n_TRUCKS_CU_SH <- metadata$vehicles[grep(pattern = "TRUCKS_CU_SH", x = metadata$vehicles)] +n_TRUCKS_CU_LH <- metadata$vehicles[grep(pattern = "TRUCKS_CU_LH", x = metadata$vehicles)] + +n_veh <- list( + MC = n_MC, + PC = n_PC, + PT = n_PT, + LCT = n_LCT, + BUS_INTERCITY = n_BUS_INTERCITY, + BUS_TRANSIT= n_BUS_TRANSIT, + BUS_SCHOOL = n_BUS_SCHOOL, + TRUCKS_REFUSE = n_TRUCKS_REFUSE, + TRUCKS_SU_SH = n_TRUCKS_SU_SH, + TRUCKS_SU_LH = n_TRUCKS_SU_LH, + TRUCKS_MH = n_TRUCKS_MH, + TRUCKS_CU_SH = n_TRUCKS_CU_SH, + TRUCKS_CU_LH = n_TRUCKS_CU_LH +) + +for(i in seq_along(metadata$vehicles)) { + if(identical(names(veh)[i + 1],metadata$vehicles[i]) == FALSE) { + + switch (language, + "portuguese" = cat( "Nomes incompativeis:\n", + "metadata:",metadata$vehicles[i], "\n", + "`veh`:",names(veh)[i + 1], "\n"), + "english" = cat( "Incomptaible names:\n", + "metadata:",metadata$vehicles[i], "\n", + "`veh`:",names(veh)[i + 1], "\n"), + "spanish" = cat( "Nombres incomptatibles:\n", + "metadata:",metadata$vehicles[i], "\n", + "`veh`:",names(veh)[i + 1], "\n")) + + + stop() + } +} + + +# apagando arquivos +switch (language, + "portuguese" = message("Apagando veh/*.rds\n"), + "english" = message("Deleting veh/*.rds\n"), + "spanish" = message("Borrando veh/*.rds\n")) + +arquivos <- list.files(path = "veh", pattern = ".rds", full.names = TRUE) +file.remove(arquivos) + +# fleet age +veh[is.na(veh)] <- 0 + +veh <- remove_units(veh) + +# plotting +switch (language, + "portuguese" = cat("Plotando fluxos\n"), + "english" = cat("Plotting traffic flows\n"), + "spanish" = cat("Plotando flujos\n")) + +# list of vehicles to colplot +lveh <- list() + +# MC #### +fl <- function(x) { + lapply(seq_along(x), function(i) { + if(length(x) == 1) { + 1 + } else { + df <- veh[, x] + sum(df[i]) / sum(df, na.rm = T) + } + })-> lx + names(lx) <- x + lx +} +fl(n_MC) -> l_mc + +lapply(seq_along(l_mc), function(i) { + cat(names(l_mc)[i], " is ", n_MC[[i]], " of MC\n") +}) -> lx +cat("sum: ", sum(unlist(l_mc)), "\n\n") + +# fuel +kf <- c(k_G) #gasoline + +for (i in seq_along(n_MC)) { + x <- my_age( + x = net[[categories[1]]], + y = veh[[n_MC[i]]], + name = n_MC[i], + k = l_mc[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_MC[i], ".rds")) +} +lMC <- data.frame(MC_G = colSums(x)) + +# PC #### +fl(n_PC) -> l_pc + +lapply(seq_along(l_pc), function(i) { + cat(names(l_pc)[i], " is ", l_pc[[i]], " of PC\n") +}) -> lx +cat("sum: ", sum(unlist(l_pc)), "\n\n") + +# fuel +kf <- c(k_D, #diesel + 1, #eletric + k_E85, #E85 + k_G) #gasoline + +# colplot +lPC <- list() + +for (i in seq_along(n_PC)) { + x <- my_age( + x = net[[categories[2]]], + y = veh[[n_PC[i]]], + name = n_PC[i], + k = l_pc[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_PC[i], ".rds")) + lPC[[i]] <- colSums(x) +} + +# PT #### +fl(n_PT) -> l_pt + +lapply(seq_along(l_pt), function(i) { + cat(names(l_pt)[i], " is ", l_pt[[i]], " of PT\n") +}) -> lx +cat("sum: ", sum(unlist(l_pt)), "\n\n") + +# fuel +kf <- c(k_D, #diesel + 1, #eletric + k_E85, #E85 + k_G) #gasoline + +# colplot +lPT <- list() + +for (i in seq_along(n_PT)) { + x <- my_age( + x = net[[categories[3]]], + y = veh[[n_PT[i]]], + name = n_PT[i], + k = l_pt[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_PT[i], ".rds")) + lPT[[i]] <- colSums(x) +} + +# LCT #### +fl(n_LCT) -> l_lct + +lapply(seq_along(l_lct), function(i) { + cat(names(l_lct)[i], " is ", l_lct[[i]], " of LCT\n") +}) -> lx +cat("sum: ", sum(unlist(l_lct)), "\n\n") + +# fuel +kf <- c(k_D, #diesel + 1, #eletric + k_E85, #E85 + k_G) #gasoline + +# colplot +lLCT <- list() + +for (i in seq_along(n_LCT)) { + x <- my_age( + x = net[[categories[4]]], + y = veh[[n_LCT[i]]], + name = n_LCT[i], + k = l_lct[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_LCT[i], ".rds")) + lLCT[[i]] <- colSums(x) +} + +# BUS_INTERCITY #### +fl(n_BUS_INTERCITY) -> l_bus_intercity + +lapply(seq_along(l_bus_intercity), function(i) { + cat(names(l_bus_intercity)[i], " is ", l_bus_intercity[[i]], " of BUS_INTERCITY\n") +}) -> lx +cat("sum: ", sum(unlist(l_bus_intercity)), "\n\n") + +# fuel +kf <- c(k_D) #diesel + +# colplot +lBI <- list() + +for (i in seq_along(n_BUS_INTERCITY)) { + x <- my_age( + x = net[[categories[5]]], + y = veh[[n_BUS_INTERCITY[i]]], + name = n_BUS_INTERCITY[i], + k = l_bus_intercity[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_BUS_INTERCITY[i], ".rds")) + lBI[[i]] <- colSums(x) +} + + +# BUS_TRANSIT #### +fl(n_BUS_TRANSIT) -> l_bus_transit + +lapply(seq_along(l_bus_transit), function(i) { + cat(names(l_bus_transit)[i], " is ", l_bus_transit[[i]], " of TRANSIT\n") +}) -> lx +cat("sum: ", sum(unlist(l_bus_transit)), "\n\n") + +# fuel +kf <- c(k_CNG, #cng + k_D, #gasoline + k_G) #diesel + +# colplot +lBT <- list() + +for (i in seq_along(n_BUS_TRANSIT)) { + x <- my_age( + x = net[[categories[6]]], + y = veh[[n_BUS_TRANSIT[i]]], + name = n_BUS_TRANSIT[i], + k = l_bus_transit[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_BUS_TRANSIT[i], ".rds")) + lBT[[i]] <- colSums(x) +} + + +# BUS_SCHOOL #### +fl(n_BUS_SCHOOL) -> l_bus_school + +lapply(seq_along(l_bus_school), function(i) { + cat(names(l_bus_school)[i], " is ", l_bus_school[[i]], " of SCHOOL\n") +}) -> lx +cat("sum: ", sum(unlist(l_bus_school)), "\n\n") + +# fuel +kf <- c(k_D, #gasoline + k_G) #diesel + +# colplot +lBS <- list() + +for (i in seq_along(n_BUS_SCHOOL)) { + x <- my_age( + x = net[[categories[7]]], + y = veh[[n_BUS_SCHOOL[i]]], + name = n_BUS_SCHOOL[i], + k = l_bus_school[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_BUS_SCHOOL[i], ".rds")) + lBS[[i]] <- colSums(x) +} + + +# TRUCKS_REFUSE #### +fl(n_TRUCKS_REFUSE) -> l_trucks_refuse + +lapply(seq_along(l_trucks_refuse), function(i) { + cat(names(l_trucks_refuse)[i], " is ", l_trucks_refuse[[i]], " of TRUCKS_SU\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_refuse)), "\n\n") + +# fuel +kf <- c(k_D, #diesel + k_G) #gasoline + +# colplot +lTR <- list() + +for (i in seq_along(n_TRUCKS_REFUSE)) { + x <- my_age( + x = net[[categories[8]]], + y = veh[[n_TRUCKS_REFUSE[i]]], + name = n_TRUCKS_REFUSE[i], + k = l_trucks_refuse[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_REFUSE[i], ".rds")) + lTR[[i]] <- colSums(x) +} + + +# TRUCKS_SU_SH #### +fl(n_TRUCKS_SU_SH) -> l_trucks_su_sh + +lapply(seq_along(l_trucks_su_sh), function(i) { + cat(names(l_trucks_su_sh)[i], " is ", l_trucks_su_sh[[i]], " of TRUCKS_SU\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_su_sh)), "\n\n") + +# fuel +kf <- c(k_D, #diesel + k_G) #gasoline + +# colplot +lTSS <- list() + +for (i in seq_along(n_TRUCKS_SU_SH)) { + x <- my_age( + x = net[[categories[9]]], + y = veh[[n_TRUCKS_SU_SH[i]]], + name = n_TRUCKS_SU_SH[i], + k = l_trucks_su_sh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_SU_SH[i], ".rds")) + lTSS[[i]] <- colSums(x) +} + +# TRUCKS_SU_LH #### +fl(n_TRUCKS_SU_LH) -> l_trucks_su_lh + +lapply(seq_along(l_trucks_su_lh), function(i) { + cat(names(l_trucks_su_lh)[i], " is ", l_trucks_su_lh[[i]], " of TRUCKS_SU\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_su_lh)), "\n\n") + +# fuel +kf <- c(k_D, #diesel + k_G) #gasoline + +# colplot +lTSL <- list() + +for (i in seq_along(n_TRUCKS_SU_LH)) { + x <- my_age( + x = net[[categories[10]]], + y = veh[[n_TRUCKS_SU_LH[i]]], + name = n_TRUCKS_SU_LH[i], + k = l_trucks_su_lh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_SU_LH[i], ".rds")) + lTSL[[i]] <- colSums(x) +} + + +# TRUCKS_MH #### +fl(n_TRUCKS_MH) -> l_trucks_mh + +lapply(seq_along(l_trucks_mh), function(i) { + cat(names(l_trucks_mh)[i], " is ", l_trucks_mh[[i]], " of TRUCKS_SU\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_mh)), "\n\n") + +# fuel +kf <- c(k_D, #diesel + k_G) #gasoline + +# colplot +lMH <- list() + +for (i in seq_along(n_TRUCKS_MH)) { + x <- my_age( + x = net[[categories[11]]], + y = veh[[n_TRUCKS_MH[i]]], + name = n_TRUCKS_MH[i], + k = l_trucks_mh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_MH[i], ".rds")) + lMH[[i]] <- colSums(x) +} + + +# TRUCKS_CU_SH #### +fl(n_TRUCKS_CU_SH) -> l_trucks_cu_sh + +lapply(seq_along(l_trucks_cu_sh), function(i) { + cat(names(l_trucks_cu_sh)[i], " is ", l_trucks_cu_sh[[i]], " of TRUCKS_SU\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_cu_sh)), "\n\n") + +# fuel +kf <- c(k_D, #diesel + k_G) #gasoline + +# colplot +lTCS <- list() + +for (i in seq_along(n_TRUCKS_CU_SH)) { + x <- my_age( + x = net[[categories[12]]], + y = veh[[n_TRUCKS_CU_SH[i]]], + name = n_TRUCKS_CU_SH[i], + k = l_trucks_cu_sh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_CU_SH[i], ".rds")) + lTCS[[i]] <- colSums(x) +} + + +# TRUCKS_CU_LH #### +fl(n_TRUCKS_CU_LH) -> l_trucks_cu_lh + +lapply(seq_along(l_trucks_cu_lh), function(i) { + cat(names(l_trucks_cu_lh)[i], " is ", l_trucks_cu_lh[[i]], " of TRUCKS_SU\n") +}) -> lx +cat("sum: ", sum(unlist(l_trucks_cu_lh)), "\n\n") + +# fuel +kf <- c(k_D) #diesel + +# colplot +lTCL <- list() + +for (i in seq_along(n_TRUCKS_CU_LH)) { + x <- my_age( + x = net[[categories[13]]], + y = veh[[n_TRUCKS_CU_LH[i]]], + name = n_TRUCKS_CU_LH[i], + k = l_trucks_cu_lh[[i]] * kf[i], + verbose = verbose + ) + saveRDS(x, paste0("veh/", n_TRUCKS_CU_LH[i], ".rds")) + lTCL[[i]] <- colSums(x) +} + + +# colplot #### + +DF <- cbind(lMC, + do.call("cbind", lPC), + do.call("cbind", lPT), + do.call("cbind", lLCT), + do.call("cbind", lBI), + do.call("cbind", lBT), + do.call("cbind", lBS), + do.call("cbind", lTR), + do.call("cbind", lTSS), + do.call("cbind", lTSL), + do.call("cbind", lMH), + do.call("cbind", lTCS), + do.call("cbind", lTCL)) + +names(DF) <- metadata$vehicles + +for(i in seq_along(n_veh)) { + png(paste0("images/FLOW_", n_veh[[i]], ".png"), 2000, 1500, "px", res = 300) + + colplot( + df = DF, + cols = n_veh[[i]], + xlab = "Age", + ylab = "veh/h", + main = paste0("FLOW ", n_veh[[i]]), + type = "l", + pch = NULL, + theme = theme, + spl = 9 + ) + dev.off() +} + + + +# plots #### +switch(language, + "portuguese" = cat("Plotando composição veicular\n"), + "english" = cat("Plotting vehicular composition\n"), + "spanish" = cat("Plotando composición vehicular\n") +) + + +f <- list.files(path = "veh", pattern = ".rds", full.names = T) +na <- list.files(path = "veh", pattern = ".rds") +na <- gsub(pattern = ".rds", replacement = "", x = na) +ff <- unlist(lapply(seq_along(f), function(i) { + sum(readRDS(f[i])) +})) +df <- data.frame( + per = ff, + Categoria = na +) +df <- df[df$Categoria != "fleet_age", ] +df$per <- df$per / sum(df$per) * 100 +df$Categoria <- factor(df$Categoria, levels = metadata$vehicles) +p <- ggplot(df, aes(x = Categoria, y = per, fill = per)) + + geom_bar(stat = "identity", col = "black") + + labs(y = "[%]", title = "Vehicular composition") + + ggplot2::scale_fill_gradientn("[%]", colours = cptcity::cpt(rev = T)) + + theme_bw() + + coord_flip() + + scale_x_discrete(limits = rev(metadata$vehicles)) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + +png( + filename = paste0("images/TRAFFIC.png"), + width = 2500, height = 2500, units = "px", pointsize = 12, + bg = "white", res = 300 +) +print(p) +dev.off() + + +if(speed) { + + + rbindlist(lapply(seq_along(metadata$vehicles), function(i) { + temp_fact(q = rowSums(readRDS(paste0("veh/", metadata$vehicles[i], ".rds"))), + pro = tfs[[metadata$vehicles[i]]]) -> tf + tf$id <- 1:nrow(tf) + tf + })) -> lx + + lx[, + lapply(.SD, sum, na.rm = T), + .SDcols = paste0("V", 1:length(tfs[[1]])), + by = id] -> ft + + ft$id <- NULL + + switch(language, + "portuguese" = cat("Projetando velocidade\n"), + "english" = cat("Projecting speed\n"), + "spanish" = cat("Proyectando velocidad\n") + ) + + df <- netspeed( + q = ft, + ps = net$CONGSPD, + ffs = net$FFSPEED, + cap = net$CAPE, + lkm = net$lkm, + alpha = 1, + dist = "miles" + ) + saveRDS(df, "network/speed.rds") + + dfs <- netspeed( + q = ft, + ps = net[[ps]], + ffs = net[[ffs]], + cap = net[[capacity]], + lkm = net$lkm, + alpha = 1, + dist = "miles", + net = net + ) + + dfs <- st_transform(dfs, 4326) + saveRDS(dfs, "network/speed_net.rds") + + + sp_mov <- moves_speed(df) + saveRDS(sp_mov, "network/speed_bin.rds") + + + # plots #### + switch(language, + "portuguese" = cat("Plotando velocidade\n"), + "english" = cat("Plotting speed\n"), + "spanish" = cat("Plotando velocidad\n") + ) + + png( + filename = paste0("images/SPEED.png"), + width = 2500, height = 3500, units = "px", pointsize = 12, + bg = "white", res = 300 + ) + plot(df, mai2 =c(2.2, 0.82, 0.82, 0.42)) + dev.off() + + png( + filename = paste0("images/SPEED_HOUR.png"), + width = 2500, height = 3000, units = "px", pointsize = 12, + bg = "white", res = 300 + ) + plot(dfs, + max.plot = 25, + axes = T, + pal = cpt(colorRampPalette = T, rev = T), + key.pos = 4) + + dev.off() + + + png( + filename = paste0("images/SPEEDBIN.png"), + width = 2500, height = 3500, units = "px", pointsize = 12, + bg = "white", res = 300 + ) + plot(Speed(remove_units(sp_mov), dist = units::unitless), pal = "mpl_viridis") + dev.off() + + +} + +switch(language, + "portuguese" = message("\nArquivos em: /veh:"), + "english" = message("\nFiles in: /veh"), + "spanish" = message("\nArchivos en: /veh") +) + +switch(language, + "portuguese" = message("Figuras em: /images:"), + "english" = message("Figures in: /images"), + "spanish" = message("Figuras en: /images") +) + + +invisible(gc()) diff --git a/projects/moves/scripts/wrf.R b/projects/moves/scripts/wrf.R new file mode 100755 index 00000000..881c659e --- /dev/null +++ b/projects/moves/scripts/wrf.R @@ -0,0 +1,238 @@ + + +# OS +sep <- ifelse(Sys.info()[["sysname"]] == "Windows", "00%3A", ":") + +wrfi <- paste0(dir_wrfinput, "/wrfinput_d0", domain) + +# tempos +ti <- wrf_get(file = wrfi, name = "Times") + +switch(language, + "portuguese" = cat("Primer tempo de WRF:", as.character(ti)[1], "\n"), + "english" = cat("First WRF time:", as.character(ti)[1], "\n"), + "spanish" = cat("Primer tiempo de WRF:", as.character(ti)[1], "\n") +) + +# ltemissions, first monday 00:00 before ti +timepos <- as.POSIXct(as.character(ti), format = "%Y-%m-%d_%H:%M:%S") +oneweek <- timepos - 3600 * 24 * 7 + +df_time <- data.frame( + times = seq.POSIXt(from = oneweek, to = timepos, by = "hour") +) + +df_time$wday <- strftime(df_time$times, "%u") +df_time$hour <- hour(df_time$times) +lt_emissions <- df_time[df_time$wday == 1 & df_time$hour == 0, ]$times[1] + +switch(language, + "portuguese" = cat( + "Segunda-feira 00:00 anterior do primeiro tempo WRF:", + as.character(lt_emissions), "\n" + ), + "english" = cat( + "Monday 00:00, previous of first WRF time:", + as.character(lt_emissions), "\n" + ), + "spanish" = cat( + "Lunes 00:00 antes del primer tiempo WRF:", + as.character(lt_emissions), "\n" + ) +) + + +ti <- as.POSIXct(ti) + +dir_mech <- paste0("post/", mech) + +# creating directory wrf/mech +dir_wrfchemi <- paste0(dir_wrfchemi, "/", mech) +dir.create(path = dir_wrfchemi) + +# Grades +lf <- list.files(path = dir_mech, pattern = ".rds", full.names = TRUE) +na <- list.files(path = dir_mech, pattern = ".rds", full.names = F) +emis_option <- na <- gsub(".rds", "", na) + + +if (mech %in% c("MOZT1", "CBMZ")) { + emis_option <- ifelse(emis_option == "E_ETOH", "E_C2H5OH", emis_option) + na <- ifelse(na == "E_ETOH", "E_C2H5OH", na) +} + +if (io_style_emissions == 1) { + # emissions on 00z / 12z style, create the 12z #### + wrfc <- eixport::wrf_create( + wrfinput_dir = dir_wrfinput, + wrfchemi_dir = dir_wrfchemi, + domains = domain, + io_style_emissions = 1, + day_offset = 0, + variables = emis_option, + verbose = TRUE, + n_aero = n_aero, + return_fn = TRUE + ) + + wrfc <- eixport::wrf_create( + wrfinput_dir = dir_wrfinput, + wrfchemi_dir = dir_wrfchemi, + domains = domain, + io_style_emissions = 1, + day_offset = 0.5, + variables = emis_option, + verbose = TRUE, + n_aero = n_aero, + return_fn = TRUE + ) + + for (i in 1:length(na)) { + print(na[i]) + x <- readRDS(lf[i]) + + xx <- emis_order( + x = x, + lt_emissions = lt_emissions, + start_utc_time = ti, + desired_length = wrf_times, + tz_lt = Sys.timezone(), + seconds = hours * 3600, + verbose = TRUE + ) + + # 0-12 + gx <- GriddedEmissionsArray( + x = xx[, 1:12], + cols = cols, # if towdown, cols + rows = rows, # if towdown, rows + times = 12, + rotate = "cols" + ) + wrf_put( + file = paste0(dir_wrfchemi, "/wrfchemi_00z_d0", domain), + name = na[i], + POL = gx + ) + + # 12-0 + gx <- GriddedEmissionsArray( + x = xx[, 13:24], + cols = cols, # if towdown, cols + rows = rows, # if towdown, rows + times = 12, + rotate = "cols" + ) + wrf_put( + file = paste0(dir_wrfchemi, "/wrfchemi_12z_d0", domain), + name = na[i], + POL = gx + ) + } +} else if (io_style_emissions == 2) { + # emissions for all hours #### + + wrfc <- wrf_create( + wrfinput_dir = dir_wrfinput, + wrfchemi_dir = paste0("wrf/", mech), + io_style_emissions = io_style_emissions, + domains = domain, + frames_per_auxinput5 = wrf_times, + auxinput5_interval_m = 60, + variables = emis_option, + separator = "_", + verbose = TRUE, + return_fn = TRUE + ) + + + for (i in 1:length(na)) { + print(na[i]) + x <- readRDS(lf[i]) + + xx <- emis_order( + x = x, + lt_emissions = lt_emissions, + start_utc_time = ti, + desired_length = wrf_times, + tz_lt = Sys.timezone(), + seconds = hours * 3600, + verbose = TRUE + ) + gx <- GriddedEmissionsArray( + x = xx, + cols = cols, + rows = rows, + times = wrf_times, + rotate = rotate, + ) + wrf_put( + file = wrfc, + name = na[i], + POL = gx + ) + } +} else { + stop("io_style_emissions can have values 1 (0-12z files) or 2 (all hours file)") +} + +message(paste0("WRFCHEMI in wrf/", wrfc)) + + +png( + filename = paste0("WRF_NO0.png"), + width = 2100, height = 1500, units = "px", pointsize = 12, + bg = "white", res = 300 +) +a <- wrf_get(wrfc, "E_NO", as_raster = T) +a <- a[[1]] +a[] <- ifelse(a[] <= 0, NA, a[]) +print(sp::spplot(a, + col.regions = cpt(rev = T), + main = paste(names(a), "mol/km²/h"), + scales = list(Draw = T) +)) +dev.off() + +png( + filename = paste0("WRF_NO12.png"), + width = 2100, height = 1500, units = "px", pointsize = 12, + bg = "white", res = 300 +) +a <- wrf_get(wrfc, "E_NO", as_raster = T) +a <- a[[1]] +a[] <- ifelse(a[] <= 0, NA, a[]) +print(sp::spplot(a, + col.regions = cpt(rev = T), + main = paste(names(a), "mol/km²/h"), + scales = list(Draw = T) +)) +dev.off() + +png( + filename = paste0("WRF_CO0.png"), + width = 2100, height = 1500, units = "px", pointsize = 12, + bg = "white", res = 300 +) +a <- wrf_get(wrfc, "E_CO", as_raster = T) +a <- a[[1]] +a[] <- ifelse(a[] <= 0, NA, a[]) +print(sp::spplot(a, + col.regions = cpt(rev = T), + main = paste(names(a), "mol/km²/h"), + scales = list(Draw = T) +)) +dev.off() + +ls() + +suppressWarnings( + rm( + "cols", "dfwrf", "domain", "dx", "emis_opt", "emis_option", + "fCO", "fHC", "firsthour", "fNOx", "fPM10", "fPM2.5", "g_ch3oh", + "g_eth", "g_hc3", "g_hc5", "g_hc8", "g_iso", "g_ket", "g_ol2", "g_oli", + "g_olt", "g_tol", "g_xyl", "grids", "i", "lasthour", "mm_no2", "mm_x", + "net", "no2_mol", "dir_wrfchemi", "dir_wrfinput", "polss", "rows", "ti", + "voc", "vocB5EX", "vocE100EV", "vocE100EX", "vocE25EV", "vocE25EX", "wrf_times", "wrfi", "x", "xx" + ) +) \ No newline at end of file diff --git a/tests/testthat/test-fuel_corr.R b/tests/testthat/test-fuel_corr.R index 360f4a91..4b31d7c2 100644 --- a/tests/testthat/test-fuel_corr.R +++ b/tests/testthat/test-fuel_corr.R @@ -40,3 +40,9 @@ test_that("fuel_corr works", { expect_equal(fuel_corr(euro = "VIc")$LDVG$CO[[1]], 1) }) + + +test_that("fuel_corr works", { + expect_equal(fuel_corr(euro = c("PRE", "I"))$value[1], + 1) +}) http://atmoschem.github.io/vein//reference/ef_evap.html