diff --git a/src/main/R/commuter-analysis.R b/src/main/R/commuter-analysis.R new file mode 100644 index 0000000..08405ff --- /dev/null +++ b/src/main/R/commuter-analysis.R @@ -0,0 +1,52 @@ +library(matsim) +library(tidyverse) +library(readr) +library(sf) + +#FILES +FILE_DIR = "../../shared-svn/projects/DiTriMo/data/commuters-by-town" +SIM <- paste0(FILE_DIR, "lausitz-v1.0-commuter.csv") +GEMEINDE <- paste0(FILE_DIR, "pgemeinden.csv") +COMMUTER <- paste0(FILE_DIR, "commuter.csv") +SHP <- paste0(FILE_DIR, "VG5000_GEM/VG5000_GEM.shp") +LAUSITZ.SHP <- paste0(FILE_DIR, "network-area/network-area.shp") + +shp <- st_read(SHP) +lausitz.shp <- st_read(LAUSITZ.SHP) + +sim <- read_csv(file = SIM) +gemeinden <- read_csv( file = GEMEINDE) %>% + mutate(code = str_remove(string = code, pattern = "P")) + +commuter <- read_csv(file = COMMUTER) %>% + mutate(key = paste0(from, "-", to)) + +sim.1 <- sim %>% + filter(!is.na(from) & !is.na(to)) %>% + left_join(gemeinden, by = c("from" = "code")) %>% + left_join(gemeinden, by = c("from" = "code"), suffix = c("_from", "_to")) %>% + mutate(key = paste0(from, "-", to)) %>% + select(-c(from, to)) %>% + left_join(commuter, by = "key", suffix = c("_sim", "_real")) %>% + filter(!is.na(from)) %>% +# pivot_longer(cols = starts_with("n_"), names_to = "src", values_to = "n", names_prefix = "n_") %>% + arrange(desc(n_real)) + +breaks <- c(-Inf, 0.8, 1.2, Inf) +labels <- c("less", "exakt", "more") + +sim.2 <- sim.1 %>% + filter(n_sim > 10) %>% + select(from, to, starts_with("n_"), starts_with("name_")) %>% + mutate(n_rel = n_sim / n_real, + quality = cut(n_rel, breaks = breaks, labels = labels)) + +ggplot(sim.2, aes(x = n_real, y = n_sim, col = quality)) + + + geom_point() + + + scale_x_log10() + + + scale_y_log10() + + + theme_bw() diff --git a/src/main/R/counts.R b/src/main/R/counts.R index 46adb56..82a1354 100644 --- a/src/main/R/counts.R +++ b/src/main/R/counts.R @@ -4,10 +4,10 @@ devtools::install_github("matsim-vsp/matsim-r",ref="counts") library(matsim) library(tidyverse) -COUNTS <- "Y:/matsim-lausitz/input/v1.0/lausitz-v1.0-counts-car-bast.xml.gz" -NETWORK <- "Y:/matsim-lausitz/input/v1.0/lausitz-v1.0-network-with-pt.xml.gz" +COUNTS <- "../../public-svn/matsim/scenarios/countries/de/lausitz/input/v1.0/lausitz-v1.0-counts-car-bast.xml.gz" +NETWORK <- "../../public-svn/matsim/scenarios/countries/de/lausitz/input/v1.0/lausitz-v1.0-network-with-pt.xml.gz" -linkstats <- readLinkStats(runId = "v1.0-uncalibrated", file = "Y:/matsim-lausitz/qa/output/lausitz-25pct.output_linkstats.csv.gz") +linkstats <- readLinkStats(runId = "v1.0-uncalibrated", file = "Y:/matsim-lausitz/qa/output/lausitz-25pct.output_linkstats.csv.gz", sampleSize = 1) counts <- readCounts(COUNTS) network <- loadNetwork(NETWORK) @@ -17,7 +17,7 @@ join <- mergeCountsAndLinks(counts = counts, linkStats = list(linkstats), netwo #### VIA-styled scatterplot #### -FILE_DIR = "C:/Users/ACER/Desktop/Uni/VSP/Lausitz-Plots/" +FILE_DIR = "../../shared-svn/projects/DiTriMo/data/commuters-by-town" createCountScatterPlot(joinedFrame = join) ggsave(filename = paste0(FILE_DIR, "Traffic_Count_Scatterplot_with_freight.jpg")) @@ -43,7 +43,7 @@ rm(join.dtv.distribution) #### Analysis of Estimation Quality #### -join.est.quality <- processDtvEstimationQuality(joinedFrame = join, aggr = T) %>% +join.est.quality <- processDtvEstimationQuality(joinedFrame = join, aggr = F) %>% filter(!type %in% c("residential", "unclassified", NA)) ggplot(join.est.quality, aes(estimation, share, fill = type)) + @@ -59,4 +59,26 @@ ggplot(join.est.quality, aes(estimation, share, fill = type)) + theme(legend.position = "none", axis.text.x = element_text(angle = 90)) rm(join.est.quality) -ggsave(filename = paste0(FILE_DIR, "Estimation_quality_by_road_type_with_freight.jpg")) \ No newline at end of file +ggsave(filename = paste0(FILE_DIR, "Estimation_quality_by_road_type_with_freight.jpg")) + + +#### network plot #### + +library(tmap) +library(tmaptools) +library(OpenStreetMap) +library(sf) + +link.geom <- join %>% + left_join(network$links, by = c("loc_id" = "id")) %>% + mutate(geom = sprintf("LINESTRING(%s %s, %s %s)", x.from, y.from, x.to, y.to)) %>% + st_as_sf(crs = 25832, wkt = "geom") %>% + transmute(loc_id, type.x, rel_vol = volume / count, geom) + +tmap_mode("view") + +tm_shape(shp = link.geom) + + tm_lines(col = "estimation", style = "cont", lwd = 3.5, palette = c("red", "green", "blue")) + +tm_shape(shp = link.geom) + + tm_lines(col = "rel_vol", style = "cont", lwd = 5, palette = c("red", "yellow", "green"), breaks = c(0, 0.05, 0.8, 2)) diff --git a/src/main/java/org/matsim/run/LausitzScenario.java b/src/main/java/org/matsim/run/LausitzScenario.java index c5b02b7..4a16f7c 100644 --- a/src/main/java/org/matsim/run/LausitzScenario.java +++ b/src/main/java/org/matsim/run/LausitzScenario.java @@ -34,7 +34,9 @@ import org.matsim.core.population.PopulationUtils; import org.matsim.core.replanning.strategies.DefaultPlanStrategiesModule; import org.matsim.core.scoring.functions.ScoringParametersForPerson; +import org.matsim.run.analysis.CommunityFilter; import org.matsim.run.analysis.CommuterAnalysis; +import org.matsim.run.analysis.DistanceMatrix; import org.matsim.run.prepare.PrepareNetwork; import org.matsim.run.prepare.PreparePopulation; import org.matsim.simwrapper.SimWrapperConfigGroup; @@ -59,7 +61,7 @@ SplitActivityTypesDuration.class, CreateCountsFromBAStData.class, PreparePopulation.class, CleanPopulation.class, PrepareNetwork.class }) @MATSimApplication.Analysis({ - LinkStats.class, CheckPopulation.class, CommuterAnalysis.class, + LinkStats.class, CheckPopulation.class, CommuterAnalysis.class, CommunityFilter.class, DistanceMatrix.class }) public class LausitzScenario extends MATSimApplication { diff --git a/src/main/java/org/matsim/run/analysis/CommunityFilter.java b/src/main/java/org/matsim/run/analysis/CommunityFilter.java new file mode 100644 index 0000000..1bd1b90 --- /dev/null +++ b/src/main/java/org/matsim/run/analysis/CommunityFilter.java @@ -0,0 +1,80 @@ +package org.matsim.run.analysis; + +import org.apache.commons.csv.CSVPrinter; +import org.geotools.api.feature.simple.SimpleFeature; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.Point; +import org.matsim.application.MATSimAppCommand; +import org.matsim.application.options.CsvOptions; +import org.matsim.core.utils.collections.Tuple; +import org.matsim.core.utils.gis.GeoFileReader; +import picocli.CommandLine; + +import java.nio.file.Path; +import java.util.*; + +@CommandLine.Command(name = "community-filter", description = "creates a csv with commuity keys within shape") +public class CommunityFilter implements MATSimAppCommand { + + @CommandLine.Option(names = "--community-shp", description = "path to VG5000_GEM.shp", required = true) + private static Path communityShapePath; + + @CommandLine.Option(names = "--dilution-area", description = "path to area-of-interest-shape file", required = true) + private static Path dilutionAreaShapePath; + + @CommandLine.Option(names = "--output", description = "output csv filepath", required = true) + private static Path output; + + @CommandLine.Mixin + CsvOptions csvOptions = new CsvOptions(); + private final Map> filtered = new HashMap<>(); + + public static void main(String[] args) { + new CommunityFilter().execute(args); + } + + @Override + public Integer call() throws Exception { + + Collection communities = readShapeFile(communityShapePath); + Collection dilutionArea = readShapeFile(dilutionAreaShapePath); + + List geometries = dilutionArea.stream().map(feature -> (Geometry) feature.getDefaultGeometry()).toList(); + + for(var community: communities){ + + String attribute = (String) community.getAttribute("ARS"); + + if(geometries.stream() + .anyMatch(geometry -> geometry.covers((Geometry) community.getDefaultGeometry()))){ + + Point centroid = ((Geometry) community.getDefaultGeometry()).getCentroid(); + Tuple coordinates = Tuple.of(centroid.getX(), centroid.getY()); + + filtered.put(attribute, coordinates); + } + + } + + try (CSVPrinter printer = csvOptions.createPrinter(output)) { + printer.print("ars"); + printer.print("x"); + printer.print("y"); + printer.println(); + + for (Map.Entry> entry : filtered.entrySet()) { + printer.print(entry.getKey()); + printer.print(entry.getValue().getFirst()); + printer.print(entry.getValue().getSecond()); + printer.println(); + } + } + + return 0; + } + + private static Collection readShapeFile(Path filepath){ + + return GeoFileReader.getAllFeatures(filepath.toString()); + } +} diff --git a/src/main/java/org/matsim/run/analysis/DistanceMatrix.java b/src/main/java/org/matsim/run/analysis/DistanceMatrix.java new file mode 100644 index 0000000..3be4a0d --- /dev/null +++ b/src/main/java/org/matsim/run/analysis/DistanceMatrix.java @@ -0,0 +1,108 @@ +package org.matsim.run.analysis; + +import org.apache.commons.csv.CSVPrinter; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; +import org.geotools.api.feature.simple.SimpleFeature; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.Point; +import org.matsim.api.core.v01.Coord; +import org.matsim.application.MATSimAppCommand; +import org.matsim.application.options.CsvOptions; +import org.matsim.application.options.ShpOptions; +import org.matsim.core.utils.geometry.CoordUtils; +import org.matsim.core.utils.geometry.geotools.MGC; +import org.matsim.core.utils.gis.GeoFileReader; +import picocli.CommandLine; + +import java.nio.file.Path; +import java.util.ArrayList; +import java.util.List; +import java.util.function.Predicate; + +@CommandLine.Command(name = "distance-matrix", description = "creates a csv with commuity keys within shape") +public class DistanceMatrix implements MATSimAppCommand{ + + @CommandLine.Option(names = "--output", description = "output csv filepath", required = true) + private static Path output; + + @CommandLine.Option(names = "--dilution-area", description = "shape to filter zones", required = false) + private static Path dilutionArea; + + @CommandLine.Mixin + CsvOptions csvOptions = new CsvOptions(); + + @CommandLine.Mixin + ShpOptions shp = new ShpOptions(); + + private final List distances = new ArrayList<>(); + private static final Logger logger = LogManager.getLogger(DistanceMatrix.class); + + public static void main(String[] args) { + new DistanceMatrix().execute(args); + } + + @Override + public Integer call() throws Exception { + + logger.info("Read features."); + List communities = shp.readFeatures(); + + //to prevent RuntimeExceptions + ArrayList copy = new ArrayList<>(communities); + + Predicate filter = getFilter(dilutionArea); + + logger.info("Calculate distance matrix."); + String delimiter = csvOptions.getFormat().getDelimiterString(); + for(var community: communities){ + + if(!filter.test((Geometry) community.getDefaultGeometry())) + continue; + String nameFrom = (String) community.getAttribute("ARS"); + Point centroid = ((Geometry) community.getDefaultGeometry()).getCentroid(); + Coord from = MGC.point2Coord(centroid); + for(var target: copy){ + + String nameTo = (String) target.getAttribute("ARS"); + + Point centroid2 = ((Geometry) target.getDefaultGeometry()).getCentroid(); + Coord to = MGC.point2Coord(centroid2); + double distance = CoordUtils.calcEuclideanDistance(from, to); + + String distanceString = String.valueOf(distance).replace('.', ','); + + distances.add(nameFrom + delimiter + nameTo + delimiter + distanceString); + } + } + + logger.info("Print results to {}", output); + try (CSVPrinter printer = csvOptions.createPrinter(output)) { + printer.print("from"); + printer.print("to"); + printer.print("distance"); + printer.println(); + + for (String entry : distances) { + for (String col : entry.split(csvOptions.getFormat().getDelimiterString())) + printer.print(col); + printer.println(); + } + } + + logger.info("Done!"); + return 0; + } + + private Predicate getFilter(Path path){ + + if(path == null) + return community -> true; + + List geometries = GeoFileReader.getAllFeatures(path.toString()).stream() + .map(feature -> (Geometry) feature.getDefaultGeometry()) + .toList(); + + return community -> geometries.stream().anyMatch(geometry -> geometry.covers(community)); + } +}