integrate survey data for workplace assignment

Makefile

@@ -134,9 +134,10 @@ $p/berlin-$V-counts-vmz.xml.gz: $p/berlin-$V-network.xml.gz
 	 --target-crs $(CRS)\
 	 --counts-mapping input/counts_mapping.csv
 
-$p/berlin-$V-facilities.xml.gz: $p/berlin-$V-network.xml.gz input/facilities.gpkg
+$p/berlin-$V-facilities.xml.gz: $p/berlin-$V-network.xml.gz input/facilities.gpkg $(berlin)/input/shp/Planungsraum_EPSG_25833.shp
 	$(sc) prepare facilities --network $< --shp $(word 2,$^)\
 	 --facility-mapping input/facility_mapping.json\
+	 --zones-shp $(word 3,$^)\
 	 --output $@
 
 $p/berlin-only-$V-100pct.plans.xml.gz: input/PLR_2013_2020.csv $(berlin)/input/shp/Planungsraum_EPSG_25833.shp input/facilities.gpkg
@@ -189,6 +190,7 @@ $p/berlin-initial-$V-25pct.plans.xml.gz: $p/berlin-activities-$V-25pct.plans.xml
 	 --network $(word 3,$^)\
 	 --shp $(germany)/vg5000/vg5000_ebenen_0101/VG5000_GEM.shp\
 	 --commuter $(germany)/regionalstatistik/commuter.csv\
+	 --berlin-commuter src/main/python/berlin-work-commuter.csv
 
 	# For debugging and visualization
 	$(sc) prepare downsample-population $@\

src/main/java/org/matsim/prepare/facilities/AttributedActivityFacility.java

@@ -9,8 +9,7 @@
 
 import java.util.Map;
 
-import static org.matsim.prepare.population.Attributes.ATTRACTION_OTHER;
-import static org.matsim.prepare.population.Attributes.ATTRACTION_WORK;
+import static org.matsim.prepare.population.Attributes.*;
 
 /**
  * Wraps any {@link ActivityFacility} and adds cached attributes because the normal attributes are terrible slow.
@@ -31,7 +30,7 @@ public AttributedActivityFacility(ActivityFacility facility) {
 		this.workAttraction = (double) facility.getAttributes().getAttribute(ATTRACTION_WORK);
 		this.otherAttraction = (double) facility.getAttributes().getAttribute(ATTRACTION_OTHER);
 		this.location = (String) facility.getAttributes().getAttribute("location");
-		this.zone = (String) facility.getAttributes().getAttribute("zone");
+		this.zone = (String) facility.getAttributes().getAttribute(ZONE);
 	}
 
 	public double getWorkAttraction() {

src/main/java/org/matsim/prepare/facilities/CreateMATSimFacilities.java

@@ -25,6 +25,7 @@
 import org.matsim.core.utils.geometry.geotools.MGC;
 import org.matsim.facilities.*;
 import org.matsim.prepare.population.Attributes;
+import org.matsim.run.OpenBerlinScenario;
 import picocli.CommandLine;
 
 import java.math.BigDecimal;
@@ -35,19 +36,17 @@
 import java.util.stream.Collectors;
 
 @CommandLine.Command(
-		name = "facilities",
-		description = "Creates MATSim facilities from shape-file and network"
+	name = "facilities",
+	description = "Creates MATSim facilities from shape-file and network"
 )
 public class CreateMATSimFacilities implements MATSimAppCommand {
 
-	private static final Logger log = LogManager.getLogger(CreateMATSimFacilities.class);
-
 	/**
 	 * Filter link types that don't have a facility associated.
 	 */
 	public static final Set<String> IGNORED_LINK_TYPES = Set.of("motorway", "trunk",
-			"motorway_link", "trunk_link", "secondary_link", "primary_link");
-
+		"motorway_link", "trunk_link", "secondary_link", "primary_link");
+	private static final Logger log = LogManager.getLogger(CreateMATSimFacilities.class);
 	@CommandLine.Option(names = "--network", required = true, description = "Path to car network")
 	private Path network;
 
@@ -57,6 +56,9 @@ public class CreateMATSimFacilities implements MATSimAppCommand {
 	@CommandLine.Option(names = "--facility-mapping", description = "Path to facility napping json", required = true)
 	private Path mappingPath;
 
+	@CommandLine.Option(names = "--zones-shp", description = "Path to shp file with zonal system", required = false)
+	private Path zonesPath;
+
 	@CommandLine.Mixin
 	private ShpOptions shp;
 
@@ -75,13 +77,24 @@ public static Id<ActivityFacility> generateId(ActivityFacilities facilities, Spl
 		byte[] bytes = new byte[3];
 		do {
 			rnd.nextBytes(bytes);
-			id = Id.create( "f" + HexFormat.of().formatHex(bytes), ActivityFacility.class);
+			id = Id.create("f" + HexFormat.of().formatHex(bytes), ActivityFacility.class);
 
 		} while (facilities.getFacilities().containsKey(id));
 
 		return id;
 	}
 
+	private static double round(double f) {
+		return BigDecimal.valueOf(f).setScale(4, RoundingMode.HALF_UP).doubleValue();
+	}
+
+	private static boolean hasAttribute(SimpleFeature ft, String name) {
+		return ft.getAttribute(name) != null &&
+			(Boolean.TRUE.equals(ft.getAttribute(name)) || "1".equals(ft.getAttribute(name)) ||
+				(ft.getAttribute(name) instanceof Number number && number.intValue() > 0)
+			);
+	}
+
 	@Override
 	public Integer call() throws Exception {
 
@@ -92,6 +105,9 @@ public Integer call() throws Exception {
 
 		config = new ObjectMapper().readerFor(MappingConfig.class).readValue(mappingPath.toFile());
 
+		ShpOptions zoneShp = zonesPath != null ? new ShpOptions(zonesPath, "EPSG:25833", null) : null;
+		ShpOptions.Index zones = zoneShp != null ? zoneShp.createIndex(OpenBerlinScenario.CRS, "SCHLUESSEL") : null;
+
 		Network completeNetwork = NetworkUtils.readNetwork(this.network.toString());
 		TransportModeNetworkFilter filter = new TransportModeNetworkFilter(completeNetwork);
 		Network carOnlyNetwork = NetworkUtils.createNetwork();
@@ -120,8 +136,8 @@ public Integer call() throws Exception {
 		double workUpper = work.getPercentile(75) + 1.5 * workIQR;
 		double otherUpper = other.getPercentile(75) + 1.5 * otherIQR;
 
-		double workLower = Math.max(work.getPercentile(1), work.getPercentile(25) - 1.5 * workIQR);
-		double otherLower = Math.max(other.getPercentile(1), other.getPercentile(25) - 1.5 * otherIQR);
+		double workLower = Math.max(Math.max(work.getPercentile(1), 0.1), work.getPercentile(25) - 1.5 * workIQR);
+		double otherLower = Math.max(Math.max(other.getPercentile(1), 0.1), other.getPercentile(25) - 1.5 * otherIQR);
 
 		log.info("Work IQR: {} Upper: {} Lower: {}", workIQR, workUpper, workLower);
 		log.info("Other IQR: {} Upper: {} Lower: {}", otherIQR, otherUpper, otherLower);
@@ -151,12 +167,20 @@ public Integer call() throws Exception {
 
 			// Filter outliers from the attraction
 			facility.getAttributes().putAttribute(Attributes.ATTRACTION_WORK,
-			 	round(Math.min(Math.max(h.attractionWork, workLower), workUpper))
+				round(Math.min(Math.max(h.attractionWork, workLower), workUpper))
 			);
 			facility.getAttributes().putAttribute(Attributes.ATTRACTION_OTHER,
 				round(Math.min(Math.max(h.attractionOther, otherLower), otherUpper))
 			);
 
+			if (zones != null) {
+				String zone = zones.query(facility.getCoord());
+				if (zone != null) {
+					facility.getAttributes().putAttribute(Attributes.LOR, Integer.parseInt(zone));
+					facility.getAttributes().putAttribute(Attributes.ZONE, zone.substring(0, 2));
+				}
+			}
+
 			facilities.addActivityFacility(facility);
 		}
 
@@ -168,10 +192,6 @@ public Integer call() throws Exception {
 		return 0;
 	}
 
-	private static double round(double f) {
-		return BigDecimal.valueOf(f).setScale(4, RoundingMode.HALF_UP).doubleValue();
-	}
-
 	/**
 	 * Sample points and choose link with the nearest points.
 	 */
@@ -186,8 +206,8 @@ private Holder processFeature(SimpleFeature ft, Network network) {
 		List<Id<Link>> links = coords.stream().map(coord -> NetworkUtils.getNearestLinkExactly(network, coord).getId()).toList();
 
 		Map<Id<Link>, Long> map = links.stream()
-				.filter(l -> !IGNORED_LINK_TYPES.contains(NetworkUtils.getType(network.getLinks().get(l))))
-				.collect(Collectors.groupingBy(Function.identity(), Collectors.counting()));
+			.filter(l -> !IGNORED_LINK_TYPES.contains(NetworkUtils.getType(network.getLinks().get(l))))
+			.collect(Collectors.groupingBy(Function.identity(), Collectors.counting()));
 
 		// Everything could be filtered and map empty
 		if (map.isEmpty())
@@ -203,7 +223,7 @@ private Holder processFeature(SimpleFeature ft, Network network) {
 		double area = (double) ft.getAttribute("area");
 
 		List<Map.Entry<Id<Link>, Long>> counts = map.entrySet().stream().sorted(Map.Entry.comparingByValue())
-				.toList();
+			.toList();
 
 		// The "main" link of the facility
 		Id<Link> link = counts.get(counts.size() - 1).getKey();
@@ -237,8 +257,8 @@ private List<Coord> samplePoints(MultiPolygon geometry, int n) {
 		for (int i = 0; i < max && result.size() < n; i++) {
 
 			Coord coord = CoordUtils.round(new Coord(
-					bbox.getMinX() + (bbox.getMaxX() - bbox.getMinX()) * rnd.nextDouble(),
-					bbox.getMinY() + (bbox.getMaxY() - bbox.getMinY()) * rnd.nextDouble()
+				bbox.getMinX() + (bbox.getMaxX() - bbox.getMinX()) * rnd.nextDouble(),
+				bbox.getMinY() + (bbox.getMaxY() - bbox.getMinY()) * rnd.nextDouble()
 			));
 
 			try {
@@ -271,13 +291,6 @@ private Set<String> activities(SimpleFeature ft) {
 		return act;
 	}
 
-	private static boolean hasAttribute(SimpleFeature ft, String name) {
-		return ft.getAttribute(name) != null &&
-			(Boolean.TRUE.equals(ft.getAttribute(name)) || "1".equals(ft.getAttribute(name)) ||
-				(ft.getAttribute(name) instanceof Number number && number.intValue() > 0)
-			);
-	}
-
 	/**
 	 * Temporary data holder for facilities.
 	 */

src/main/java/org/matsim/prepare/population/Attributes.java

@@ -21,18 +21,21 @@ public final class Attributes {
 	 * Gemeinde code.
 	 */
 	public static final String GEM = "gem";
-	public static final String ARS = "ars";
 
+	/**
+	 * Amtliche Regionalschlüssel (ARS).
+	 */
+	public static final String ARS = "ars";
 
 	/**
-	 * LOR for Berlin.
+	 * Lebensweltlich orientierte Räume (LOR) for Berlin. / Zonal system of Berlin.
 	 */
 	public static final String LOR = "lor";
 
 	/**
-	 * District of Berlin.
+	 * Zonal code, in Berlin this is the district number.
 	 */
-	public static final String DISTRICT = "district";
+	public static final String ZONE = "zone";
 
 	public static final String RegioStaR7 = "RegioStaR7";
 

src/main/java/org/matsim/prepare/population/CommuterAssignment.java

@@ -1,16 +1,20 @@
 package org.matsim.prepare.population;
 
+import it.unimi.dsi.fastutil.ints.Int2DoubleMap;
+import it.unimi.dsi.fastutil.ints.Int2DoubleOpenHashMap;
+import it.unimi.dsi.fastutil.ints.Int2ObjectMap;
+import it.unimi.dsi.fastutil.ints.Int2ObjectOpenHashMap;
 import it.unimi.dsi.fastutil.longs.*;
 import org.apache.commons.csv.CSVFormat;
 import org.apache.commons.csv.CSVParser;
 import org.apache.commons.csv.CSVRecord;
 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.application.options.CsvOptions;
 import org.matsim.facilities.ActivityFacility;
-import org.geotools.api.feature.simple.SimpleFeature;
 
 import java.io.IOException;
 import java.io.UncheckedIOException;
@@ -28,14 +32,19 @@ public class CommuterAssignment {
 	private final Map<Long, SimpleFeature> zones;
 
 	/**
-	 * Outgoing commuter from ars to ars.
+	 * Outgoing commuter from ars to ars. This is german wide with quite large zones.
 	 */
 	private final Long2ObjectMap<Long2DoubleMap> commuter;
 
+	/**
+	 * Maps home district to probabilities of commuting to other districts.
+	 */
+	private final Int2ObjectMap<Int2DoubleMap> berlinCommuter;
+
 	private final CsvOptions csv = new CsvOptions(CSVFormat.Predefined.Default);
 	private final double sample;
 
-	public CommuterAssignment(Long2ObjectMap<SimpleFeature> zones, Path commuterPath, double sample) {
+	public CommuterAssignment(Long2ObjectMap<SimpleFeature> zones, Path commuterPath, Path berlinCommuterPath, double sample) {
 		this.sample = sample;
 
 		// outgoing commuters
@@ -62,9 +71,32 @@ public CommuterAssignment(Long2ObjectMap<SimpleFeature> zones, Path commuterPath
 			throw new UncheckedIOException(e);
 		}
 
+		this.berlinCommuter = new Int2ObjectOpenHashMap<>();
+
+		try (CSVParser parser = csv.createParser(berlinCommuterPath)) {
+
+			for (CSVRecord row : parser) {
+				int home = Integer.parseInt(row.get("home"));
+				int work = Integer.parseInt(row.get("work"));
+
+				berlinCommuter.computeIfAbsent(home, k -> new Int2DoubleOpenHashMap())
+					.put(work, Double.parseDouble(row.get("n")));
+			}
+
+			// Normalize probabilities
+			for (Int2ObjectMap.Entry<Int2DoubleMap> kv : berlinCommuter.int2ObjectEntrySet()) {
+				Int2DoubleMap m = kv.getValue();
+				double sum = m.values().doubleStream().sum();
+				m.replaceAll((k, v) -> v / sum);
+			}
+
+		} catch (IOException e) {
+			throw new UncheckedIOException(e);
+		}
+
 	}
 
-	/**
+		/**
 	 * Select and return a commute target.
 	 *
 	 * @param f   sampler producing target locations
@@ -126,6 +158,13 @@ public ActivityFacility selectTarget(SplittableRandom rnd, long ars, double dist
 		return null;
 	}
 
+	/**
+	 * Returns the weight of commuting from homeZone to targetZone.
+	 */
+	public double getZoneWeight(String homeZone, String targetZone) {
+		return berlinCommuter.get(Integer.parseInt(homeZone)).get(Integer.parseInt(targetZone));
+	}
+
 	/**
 	 * Sample locations from specific zone.
 	 */

src/main/java/org/matsim/prepare/population/CreateBerlinPopulation.java

@@ -230,7 +230,7 @@ private void processLOR(CSVRecord row) throws ParseException {
 			person.getAttributes().putAttribute(Attributes.GEM, 11000000);
 			person.getAttributes().putAttribute(Attributes.ARS, 110000000000L);
 			person.getAttributes().putAttribute(Attributes.LOR, Integer.parseInt(raumID));
-			person.getAttributes().putAttribute(Attributes.DISTRICT, Integer.parseInt(raumID.substring(0, 2)));
+			person.getAttributes().putAttribute(Attributes.ZONE, raumID.substring(0, 2));
 
 			Plan plan = f.createPlan();
 			plan.addActivity(f.createActivityFromCoord("home", coord));