Skip to content

Commit

Permalink
Accessibility update
Browse files Browse the repository at this point in the history
  • Loading branch information
@CorinStaves
CorinStaves committed Dec 6, 2023
1 parent bd9ed0a commit 8671193
Show file tree
Hide file tree
Showing 23 changed files with 479 additions and 338 deletions.
91 changes: 32 additions & 59 deletions src/main/java/accessibility/AccessibilityWriter.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
package accessibility;

import org.apache.log4j.Logger;
import org.geotools.data.simple.SimpleFeatureReader;
import org.geotools.feature.DefaultFeatureCollection;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
Expand All @@ -18,20 +17,17 @@
import org.locationtech.jts.geom.Point;
import org.matsim.api.core.v01.Coord;
import org.matsim.api.core.v01.Id;
import org.matsim.api.core.v01.IdMap;
import org.matsim.api.core.v01.network.Network;
import org.matsim.api.core.v01.network.Node;
import org.matsim.core.utils.io.IOUtils;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.referencing.FactoryException;
import resources.Properties;
import resources.Resources;

import java.io.BufferedWriter;
import java.io.File;
import java.io.IOException;
import java.util.Collections;
import java.util.Map;
import java.util.*;
import java.util.stream.Collectors;

/**
* Helper methods to write and read matrices as CSV files (well, actually semi-colon separated files).
Expand All @@ -45,44 +41,38 @@ public final class AccessibilityWriter {
private final static char SEP = ',';
private final static char NL = '\n';

public static void writeNodesAsCsv(IdMap<Node,Double> accessibilityData, String filename) throws IOException {

double min = Collections.min(accessibilityData.values());
double max = Collections.max(accessibilityData.values());
double diff = max-min;

BufferedWriter writer = IOUtils.getBufferedWriter(filename);
writer.append("NODE" + SEP + "ACCESSIBILITY" + SEP + "NORMALISED" + NL);
for (Map.Entry<Id<Node>,Double> e : accessibilityData.entrySet()) {
writer.append(e.getKey().toString());
writer.append(SEP);
writer.append(Double.toString(e.getValue()));
writer.append(SEP);
writer.append(Double.toString((e.getValue() - min)/diff));
writer.append(NL);
}
writer.flush();
}

public static void writeNodesAsGpkg(Map<Id<Node>,Double> accessibilityData, Network network, String filename) throws IOException {
public static void writeNodesAsGpkg(Map<Id<Node>, double[]> accessibilityData, List<String> endLocationDescriptions, Network network, String filename) throws IOException {

final double min = Collections.min(accessibilityData.values());
final double max = Collections.max(accessibilityData.values());
final double diff = max-min;

final SimpleFeatureType TYPE = createNodeFeatureType();
final SimpleFeatureType TYPE = createNodeFeatureType(endLocationDescriptions);
final SimpleFeatureBuilder builder = new SimpleFeatureBuilder(TYPE);
final DefaultFeatureCollection collection = new DefaultFeatureCollection("Nodes",TYPE);

int endLocationCount = endLocationDescriptions.size();

double[] mins = new double[endLocationCount];
double[] maxs = new double[endLocationCount];
double[] diffs = new double[endLocationCount];

for(int i = 0 ; i < endLocationCount ; i++) {
int finalI = i;
Set<Double> result = accessibilityData.values().stream().map(v -> v[finalI]).collect(Collectors.toSet());
mins[i] = Collections.min(result);
maxs[i] = Collections.max(result);
diffs[i] = maxs[i] - mins[i];
}

// Convert map entries to feature data to create a feature collection
for(Map.Entry<Id<Node>,Double> e : accessibilityData.entrySet()) {
for(Map.Entry<Id<Node>,double[]> e : accessibilityData.entrySet()) {
Node node = network.getNodes().get(e.getKey());
Coord c = node.getCoord();
Point p = GEOMETRY_FACTORY.createPoint(new Coordinate(c.getX(),c.getY(),c.getZ()));
builder.add(p); // point geometry
builder.add(Integer.parseInt(e.getKey().toString())); // node ID todo: check this works
builder.add(e.getValue()); // accessibility
builder.add((e.getValue() - min) / diff); // normalised accessibility
builder.add(Integer.parseInt(e.getKey().toString())); // node ID
double[] values = e.getValue();
for(int i = 0 ; i < endLocationCount ; i++) {
builder.add(values[i]);
builder.add((values[i] - mins[i]) / diffs[i]);
}
collection.add(builder.buildFeature(null));
}

Expand All @@ -100,7 +90,7 @@ public static void writeNodesAsGpkg(Map<Id<Node>,Double> accessibilityData, Netw
out.close();
}

private static SimpleFeatureType createNodeFeatureType() {
private static SimpleFeatureType createNodeFeatureType(List<String> endLocationDescriptions) {
SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
builder.setName("Nodes");

Expand All @@ -112,30 +102,13 @@ private static SimpleFeatureType createNodeFeatureType() {
}

// add attributes in order
builder.add("Node", Point.class);
builder.add("Id",Integer.class);
builder.add("Accessibility",Double.class);
builder.add("Normalised",Double.class);

// build the type
return builder.buildFeatureType();
}

private static SimpleFeatureType createZoneFeatureType(String originalFilePath) throws IOException {

// Get base type from original file
GeoPackage geopkg = new GeoPackage(new File(originalFilePath));
SimpleFeatureReader r = geopkg.reader(geopkg.features().get(0), null,null);

SimpleFeatureType schema = r.getFeatureType();

SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
builder.setName(schema.getName());
builder.setSuperType((SimpleFeatureType) schema.getSuper());
builder.addAll(schema.getAttributeDescriptors());
builder.add("accessibility",Double.class);
builder.add("normalised",Double.class);
builder.add("node", Point.class);
builder.add("id",Integer.class);

for(String description : endLocationDescriptions) {
builder.add("accessibility_" + description,Double.class);
builder.add("normalised_" + description,Double.class);
}

// build the type
return builder.buildFeatureType();
Expand Down
120 changes: 67 additions & 53 deletions src/main/java/accessibility/FeatureCalculator.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,7 @@
import org.opengis.feature.simple.SimpleFeature;
import resources.Properties;
import resources.Resources;
import routing.graph.LeastCostPathTree3;
import routing.graph.SpeedyGraph;
import routing.graph.*;

import java.util.*;
import java.util.concurrent.ConcurrentLinkedQueue;
Expand All @@ -36,9 +35,9 @@ public class FeatureCalculator {
private final static Person PERSON = PopulationUtils.getFactory().createPerson(Id.create("thePerson", Person.class));

public static void calculate(Network routingNetwork, SimpleFeatureCollection collection,
Map<String, IdSet<Node>> endNodes, Map<String, Double> endWeights,
Map<Id<Node>,Double> nodeResults, int polygonRadius,
boolean fwd, TravelTime travelTime, TravelDisutility travelDisutility,
List<LocationData> endDataList,
Map<Id<Node>,double[]> nodeResults, int polygonRadius,
Boolean fwd, TravelTime travelTime, TravelDisutility travelDisutility,
Vehicle vehicle, DecayFunction decayFunction) {

int numberOfThreads = Resources.instance.getInt(Properties.NUMBER_OF_THREADS);
Expand Down Expand Up @@ -73,7 +72,7 @@ public static void calculate(Network routingNetwork, SimpleFeatureCollection col
Thread[] threads = new Thread[numberOfThreads];

for (int i = 0; i < numberOfThreads; i++) {
FeatureWorker worker = new FeatureWorker(featuresQueue, polygonRadius, nodeResults,endNodes,endWeights,fwd,nodesPerZone,
FeatureWorker worker = new FeatureWorker(featuresQueue, polygonRadius, nodeResults,endDataList,fwd,nodesPerZone,
routingNetwork,routingGraph, decayFunction,marginalTravelTimes,marginalDisutilities,counter);
threads[i] = new Thread(worker, "PolygonAccessibility-" + i);
threads[i].start();
Expand All @@ -89,23 +88,25 @@ public static void calculate(Network routingNetwork, SimpleFeatureCollection col
}

// normalise results
double min = features.stream().mapToDouble(c -> (double) c.getAttribute("accessibility")).min().orElseThrow();
double max = features.stream().mapToDouble(c -> (double) c.getAttribute("accessibility")).max().orElseThrow();
double diff = max - min;
for (SimpleFeature feature : features) {
double accessibility = (double) feature.getAttribute("accessibility");
feature.setAttribute("normalised",(accessibility-min) / diff);
for(LocationData endData : endDataList) {
String attributeName = "accessibility_" + endData.getDescription();
double min = features.stream().mapToDouble(c -> (double) c.getAttribute(attributeName)).min().orElseThrow();
double max = features.stream().mapToDouble(c -> (double) c.getAttribute(attributeName)).max().orElseThrow();
double diff = max - min;
for (SimpleFeature feature : features) {
double accessibility = (double) feature.getAttribute(attributeName);
feature.setAttribute("normalised_" + endData.getDescription(),(accessibility-min) / diff);
}
}
}

private static class FeatureWorker implements Runnable {

private final ConcurrentLinkedQueue<SimpleFeature> polygons;
private final ConcurrentLinkedQueue<SimpleFeature> features;
private final int zoneRadius;
private final Map<Id<Node>, Double> nodeResults;
private final Map<String, Double> endWeights;
private final Map<String, IdSet<Node>> endNodes;
private final boolean fwd;
private final Map<Id<Node>, double[]> nodeResults;
private final List<LocationData> endDataList;
private final Boolean fwd;
private final Map<Id<Link>,Double> marginalTravelTimes;
private final Map<Id<Link>,Double> marginalDisutilities;
private final Map<SimpleFeature, IdSet<Node>> nodesInPolygons;
Expand All @@ -114,17 +115,16 @@ private static class FeatureWorker implements Runnable {
private final Counter counter;
private final DecayFunction decayFunction;

FeatureWorker(ConcurrentLinkedQueue<SimpleFeature> polygons, int zoneRadius, Map<Id<Node>, Double> nodeResults,
Map<String, IdSet<Node>> endNodes, Map<String, Double> endWeights, boolean fwd,
FeatureWorker(ConcurrentLinkedQueue<SimpleFeature> features, int zoneRadius, Map<Id<Node>, double[]> nodeResults,
List<LocationData> endDataList, Boolean fwd,
Map<SimpleFeature, IdSet<Node>> nodesInPolygons, Network network,
SpeedyGraph graph, DecayFunction decayFunction,
Map<Id<Link>,Double> marginalTravelTimes, Map<Id<Link>,Double> marginalDisutilities,
Counter counter) {
this.polygons = polygons;
this.features = features;
this.zoneRadius = zoneRadius;
this.nodeResults = nodeResults;
this.endWeights = endWeights;
this.endNodes = endNodes;
this.endDataList = endDataList;
this.fwd = fwd;
this.nodesInPolygons = nodesInPolygons;
this.network = network;
Expand All @@ -136,11 +136,19 @@ private static class FeatureWorker implements Runnable {
}

public void run() {
LeastCostPathTree3 lcpTree = new LeastCostPathTree3(this.graph);
LeastCostPathTree3.StopCriterion stopCriterion = decayFunction.getTreeStopCriterion();
PathTree lcpTree;
if(fwd != null) {
log.info("Initialising 1-way least cost path tree in " + (fwd ? " FORWARD " : " REVERSE ") + " direction...");
lcpTree = new LcpTree1Way(this.graph,fwd);
} else {
log.info("Initialising 2-way least cost path tree...");
lcpTree = new LcpTree2Way(this.graph);
}

StopCriterion stopCriterion = decayFunction.getTreeStopCriterion();

while (true) {
SimpleFeature feature = this.polygons.poll();
SimpleFeature feature = this.features.poll();
if (feature == null) {
return;
}
Expand All @@ -156,7 +164,11 @@ public void run() {

// If no nodes fall inside zone, then use to & from node of nearest link
if(nodesWithin > 0) {
feature.setAttribute("accessibility", nodesInside.stream().mapToDouble(nodeResults::get).average().orElseThrow());
for(int i = 0 ; i < endDataList.size() ; i++) {
int finalI = i;
feature.setAttribute("accessibility_" + endDataList.get(finalI).getDescription(),
nodesInside.stream().mapToDouble(n -> nodeResults.get(n)[finalI]).average().orElseThrow());
}
} else {
Coord centroid = new Coord((double) feature.getAttribute("centroid_x"), (double) feature.getAttribute("centroid_y"));
calculateForPoint(feature, centroid, lcpTree, stopCriterion);
Expand All @@ -168,7 +180,7 @@ public void run() {
}
}
}
void calculateForPoint(SimpleFeature feature, Coord coord, LeastCostPathTree3 lcpTree, LeastCostPathTree3.StopCriterion stopCriterion) {
void calculateForPoint(SimpleFeature feature, Coord coord, PathTree lcpTree, StopCriterion stopCriterion) {
Link link = NetworkUtils.getNearestLinkExactly(network, coord);
Id<Link> linkId = link.getId();
double connectorMarginalCost = marginalDisutilities.get(linkId);
Expand All @@ -186,39 +198,41 @@ void calculateForPoint(SimpleFeature feature, Coord coord, LeastCostPathTree3 lc
double timeA = connectorMarginalTime * connectorLengthA;
double timeB = connectorMarginalTime * connectorLengthB;

feature.setAttribute("nodeA",nodeA.getId().toString());
feature.setAttribute("costA",costA);
feature.setAttribute("nodeB",nodeB.getId().toString());
feature.setAttribute("costB",costB);

lcpTree.calculate(
nodeA.getId().index(),costA,timeA,connectorLengthA,
nodeB.getId().index(),costB,timeB,connectorLengthB,
0.,stopCriterion,fwd);

double accessibility = 0.;

for (Map.Entry<String, Double> endWeight : this.endWeights.entrySet()) {

double cost = Double.MAX_VALUE;

for (Id<Node> toNodeId : this.endNodes.get(endWeight.getKey())) {
int toNodeIndex = toNodeId.index();
double nodeDist = lcpTree.getDistance(toNodeIndex);
double nodeTime = lcpTree.getTime(toNodeIndex).orElse(Double.POSITIVE_INFINITY);
if(decayFunction.beyondCutoff(nodeDist, nodeTime)) {
continue;
0.,stopCriterion);

for(LocationData endData: endDataList) {
Map<String, IdSet<Node>> endNodes = endData.getNodes();
Map<String, Double> endWeights = endData.getWeights();

double accessibility = 0.;
for (Map.Entry<String, Double> endWeight : endWeights.entrySet()) {
double cost = Double.MAX_VALUE;
for (Id<Node> toNodeId : endNodes.get(endWeight.getKey())) {
int toNodeIndex = toNodeId.index();
double nodeDist = lcpTree.getDistance(toNodeIndex);
double nodeTime = lcpTree.getTime(toNodeIndex).orElse(Double.POSITIVE_INFINITY);
if(decayFunction.beyondCutoff(nodeDist, nodeTime)) {
continue;
}
double nodeCost = lcpTree.getCost(toNodeIndex);
if (nodeCost < cost) {
cost = nodeCost;
}
}
double nodeCost = lcpTree.getCost(toNodeIndex);
if (nodeCost < cost) {
cost = nodeCost;
if(cost != Double.MAX_VALUE) {
accessibility += decayFunction.getDecay(cost) * endWeight.getValue();
}
}
if(cost != Double.MAX_VALUE) {
accessibility += decayFunction.getDecay(cost) * endWeight.getValue();
}
feature.setAttribute("accessibility_" + endData.getDescription(),accessibility);
}

feature.setAttribute("accessibility",accessibility);
feature.setAttribute("nodeA",nodeA.getId().toString());
feature.setAttribute("costA",costA);
feature.setAttribute("nodeB",nodeB.getId().toString());
feature.setAttribute("costB",costB);
}
}
}
10 changes: 7 additions & 3 deletions src/main/java/accessibility/FeatureData.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

import java.io.File;
import java.io.IOException;
import java.util.List;

// Tools for reading zone system for accessibility analysis

Expand All @@ -25,7 +26,7 @@ public class FeatureData {
private final Geometries geometryType;
private final Integer radius;

public FeatureData(String filePath) throws IOException {
public FeatureData(String filePath, List<String> endLocationDescriptions) throws IOException {
GeoPackage geopkg = new GeoPackage(openFile(filePath));
FeatureEntry entry = geopkg.features().get(0);
this.geometryType = entry.getGeometryType();
Expand Down Expand Up @@ -66,8 +67,11 @@ public FeatureData(String filePath) throws IOException {
builder.add("nodeB",String.class);
builder.add("costA",Double.class);
builder.add("costB",Double.class);
builder.add("accessibility",Double.class);
builder.add("normalised",Double.class);
for(String description : endLocationDescriptions) {
builder.add("accessibility_" + description,Double.class);
builder.add("normalised_" + description,Double.class);
}

SimpleFeatureType newSchema = builder.buildFeatureType();

// Create set of zones with updated feature types
Expand Down
Loading

0 comments on commit 8671193

Please sign in to comment.