ipeaGIT / gtfs2emis

R package to estimate public transport emissions based on GTFS data
https://ipeagit.github.io/gtfs2emis/
Other
28 stars 2 forks source link

Reproducible #17

Closed Joaobazzo closed 4 years ago

Joaobazzo commented 4 years ago
  library(gtfs2emis)
  # gtfs to gps to gpslinestring
  cwb_read <- gtfs2gps::read_gtfs(system.file("extdata/gtfs_cur.zip", package ="gtfs2emis"))
  cwb_gtfs <- gtfs2gps::filter_by_shape_id(gtfs_data = cwb_read,shape_ids = c(1708))
  cwb_gtfs <- gtfs2gps::gtfs2gps(gtfs_data =cwb_gtfs)
  cwb_gpslines <- gps2line(input_file = cwb_gtfs) 
  # fleet of Curitiba
  tempd <- file.path(tempdir(), "fleet") # create tempr dir to save GTFS unzipped files
  unlink(normalizePath(paste0(tempd, "/", dir(tempd)), mustWork = FALSE), recursive = TRUE) # clean tempfiles in that dir
  utils::untar(tarfile = "inst/extdata/cur_fleet.tar.xz",exdir = tempd) # untar files
  unzippedfiles <- list.files(tempd) # list of unzipped files
  cur_fleet <- data.table::fread(paste0(tempd,"/cur_fleet.txt"))
  cur_fleet <- cur_fleet[SHP %in% 1708,]
  cur_fleet <- cur_fleet[sample(x = 1:nrow(cur_fleet),size = 1),] # sample one possible vehicle
  # ---
  # emission factor
  FE_local <- vein::ef_cetesb(p = "CO",veh = "BUS_URBAN_D",year = cur_fleet$Ano_fabricacao,agemax = 1)

  FE_scaled <- gtfs2emis::ef_hdv_scaled(dfcol = FE_local,vel = units::set_units(cwb_gpslines$speed,km/h),
                                        veh = "Urban Buses Standard 15 - 18 t",fuel = "Diesel",
                                        euro = "Euro V",pol = "CO",show.equation = TRUE)
  # emission
  cwb_gpslines[,CO_emissions := units::set_units(dist,km) * FE_scaled]
  # --
  # create grid
  cwb_gpslines <- cwb_gpslines %>% sf::st_as_sf() %>% sf::st_transform(32722)
  grid_cwb <- vein::make_grid(spobj = cwb_gpslines,width = 250) #500 mts
  cwb_gridded <- vein::emis_grid(spobj = cwb_gpslines[c("CO_emissions")], g = grid_cwb)
  cwb_gridded <- cwb_gridded[as.numeric(cwb_gridded$CO_emissions) > 0,]
  # view
  ggplot()+
    geom_sf(data = cwb_gridded,aes(fill=as.numeric(CO_emissions))) + 
    geom_sf(data = cwb_gpslines$geometry)

image