Merge pull request #3167 from moia-oss/h3

Add support for H3 spatial indexing in DRT zonal systems

contribs/common/pom.xml

@@ -18,5 +18,15 @@
             <groupId>org.assertj</groupId>
             <artifactId>assertj-core</artifactId>
         </dependency>
+		<dependency>
+			<groupId>one.util</groupId>
+			<artifactId>streamex</artifactId>
+			<version>0.8.2</version>
+		</dependency>
+		<dependency>
+			<groupId>com.uber</groupId>
+			<artifactId>h3</artifactId>
+			<version>4.1.1</version>
+		</dependency>
     </dependencies>
 </project>

contribs/common/src/main/java/org/matsim/contrib/common/zones/Zone.java

@@ -0,0 +1,19 @@
+package org.matsim.contrib.common.zones;
+
+import org.locationtech.jts.geom.prep.PreparedGeometry;
+import org.matsim.api.core.v01.BasicLocation;
+import org.matsim.api.core.v01.Coord;
+import org.matsim.api.core.v01.Identifiable;
+import org.matsim.api.core.v01.network.Link;
+
+import javax.annotation.Nullable;
+import java.util.List;
+
+public interface Zone extends BasicLocation, Identifiable<Zone> {
+	@Nullable
+	PreparedGeometry getPreparedGeometry();
+
+	Coord getCentroid();
+
+	List<Link> getLinks();
+}

contribs/common/src/main/java/org/matsim/contrib/common/zones/ZoneImpl.java

@@ -0,0 +1,65 @@
+package org.matsim.contrib.common.zones;
+
+import org.locationtech.jts.geom.prep.PreparedGeometry;
+import org.matsim.api.core.v01.Coord;
+import org.matsim.api.core.v01.Id;
+import org.matsim.api.core.v01.network.Link;
+import org.matsim.core.utils.geometry.geotools.MGC;
+
+import javax.annotation.Nullable;
+import java.util.List;
+
+public class ZoneImpl implements Zone {
+
+	private final Id<Zone> id;
+	@Nullable
+	private final PreparedGeometry preparedGeometry; //null for virtual/dummy zones
+	private final List<Link> links;
+	private final Coord centroid;
+
+	public ZoneImpl(Id<Zone> id, PreparedGeometry preparedGeometry, List<Link> links) {
+		this(id, preparedGeometry, links, MGC.point2Coord(preparedGeometry.getGeometry().getCentroid()));
+	}
+
+	private ZoneImpl(Id<Zone> id, @Nullable PreparedGeometry preparedGeometry, List<Link> links, Coord centroid) {
+		this.id = id;
+		this.preparedGeometry = preparedGeometry;
+		this.links = links;
+		this.centroid = centroid;
+	}
+
+	@Override
+	public Id<Zone> getId() {
+		return id;
+	}
+
+	@Override
+	public Coord getCoord() {
+		return centroid;
+	}
+
+	@Override
+	@Nullable
+	public PreparedGeometry getPreparedGeometry() {
+		return preparedGeometry;
+	}
+
+	@Override
+	public Coord getCentroid() {
+		return centroid;
+	}
+
+	@Override
+	public List<Link> getLinks() {
+		return links;
+	}
+
+	boolean isDummy() {
+		return preparedGeometry == null;
+	}
+
+	public static ZoneImpl createDummyZone(Id<Zone> id, List<Link> links, Coord centroid) {
+		return new ZoneImpl(id, null, links, centroid);
+	}
+
+}

contribs/common/src/main/java/org/matsim/contrib/common/zones/ZoneSystem.java

@@ -0,0 +1,14 @@
+package org.matsim.contrib.common.zones;
+
+import org.matsim.api.core.v01.Id;
+import org.matsim.api.core.v01.network.Link;
+
+import javax.annotation.Nullable;
+import java.util.Map;
+
+public interface ZoneSystem {
+	@Nullable
+	Zone getZoneForLinkId(Id<Link> linkId);
+
+	Map<Id<Zone>, Zone> getZones();
+}

contribs/common/src/main/java/org/matsim/contrib/common/zones/ZoneSystemImpl.java

@@ -0,0 +1,43 @@
+package org.matsim.contrib.common.zones;
+
+import org.apache.commons.lang3.tuple.Pair;
+import org.matsim.api.core.v01.Id;
+import org.matsim.api.core.v01.IdMap;
+import org.matsim.api.core.v01.network.Link;
+
+import javax.annotation.Nullable;
+import java.util.Collection;
+import java.util.Collections;
+import java.util.Map;
+
+public class ZoneSystemImpl implements ZoneSystem {
+
+	private final Map<Id<Zone>, Zone> zones = new IdMap<>(Zone.class);
+	private final IdMap<Link, Zone> link2zone = new IdMap<>(Link.class);
+
+	public ZoneSystemImpl(Collection<Zone> zones) {
+		zones.forEach(zone -> this.zones.put(zone.getId(), zone));
+		zones.stream()
+			.flatMap(zone -> zone.getLinks().stream().map(link -> Pair.of(link.getId(), zone)))
+			.forEach(idZonePair -> link2zone.put(idZonePair.getKey(), idZonePair.getValue()));
+	}
+
+	/**
+	 * @param linkId
+	 * @return the the {@code DrtZone} that contains the {@code linkId}. If the given link's {@code Coord} borders two or more cells, the allocation to a cell is random.
+	 * Result may be null in case the given link is outside of the service area.
+	 */
+	@Override
+	@Nullable
+	public Zone getZoneForLinkId(Id<Link> linkId) {
+		return link2zone.get(linkId);
+	}
+
+	/**
+	 * @return the zones
+	 */
+	@Override
+	public Map<Id<Zone>, Zone> getZones() {
+		return Collections.unmodifiableMap(zones);
+	}
+}

contribs/common/src/main/java/org/matsim/contrib/common/zones/ZoneSystemUtils.java

@@ -0,0 +1,64 @@
+package org.matsim.contrib.common.zones;
+
+import one.util.streamex.EntryStream;
+import one.util.streamex.StreamEx;
+import org.locationtech.jts.geom.Point;
+import org.locationtech.jts.geom.prep.PreparedGeometry;
+import org.matsim.api.core.v01.Id;
+import org.matsim.api.core.v01.network.Link;
+import org.matsim.api.core.v01.network.Network;
+import org.matsim.core.utils.geometry.geotools.MGC;
+
+import javax.annotation.Nullable;
+import java.util.List;
+import java.util.Map;
+import java.util.Objects;
+
+import static java.util.stream.Collectors.toList;
+
+public final class ZoneSystemUtils {
+
+
+	private ZoneSystemUtils() {}
+
+	public static ZoneSystem createFromPreparedGeometries(Network network, Map<String, PreparedGeometry> geometries) {
+
+		//geometries without links are skipped
+		Map<String, List<Link>> linksByGeometryId = StreamEx.of(network.getLinks().values())
+			.mapToEntry(l -> getGeometryIdForLink(l, geometries), l -> l)
+			.filterKeys(Objects::nonNull)
+			.grouping(toList());
+
+		//the zonal system contains only zones that have at least one link
+		List<Zone> zones = EntryStream.of(linksByGeometryId)
+			.mapKeyValue((id, links) -> new ZoneImpl(Id.create(id, Zone.class), geometries.get(id), links))
+			.collect(toList());
+
+		return new ZoneSystemImpl(zones);
+	}
+
+	/**
+	 * @param link
+	 * @return the the {@code PreparedGeometry} that contains the {@code linkId}.
+	 * If a given link's {@code Coord} borders two or more cells, the allocation to a cell is random.
+	 * Result may be null in case the given link is outside of the service area.
+	 */
+	@Nullable
+	private static String getGeometryIdForLink(Link link, Map<String, PreparedGeometry> geometries) {
+		Point linkCoord = MGC.coord2Point(link.getToNode().getCoord());
+		return geometries.entrySet()
+			.stream()
+			.filter(e -> e.getValue().intersects(linkCoord))
+			.findAny()
+			.map(Map.Entry::getKey)
+			.orElse(null);
+	}
+
+	public static Id<Zone> createZoneId(String id) {
+		return Id.create(id, Zone.class);
+	}
+
+	public static Id<Zone> createZoneId(long id) {
+		return Id.create(id, Zone.class);
+	}
+}

contribs/common/src/main/java/org/matsim/contrib/common/zones/h3/H3GridUtils.java

@@ -0,0 +1,112 @@
+package org.matsim.contrib.common.zones.h3;
+
+import com.uber.h3core.AreaUnit;
+import com.uber.h3core.H3Core;
+import com.uber.h3core.LengthUnit;
+import com.uber.h3core.util.LatLng;
+import org.apache.logging.log4j.LogManager;
+import org.apache.logging.log4j.Logger;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.geom.prep.PreparedGeometry;
+import org.locationtech.jts.geom.prep.PreparedGeometryFactory;
+import org.matsim.api.core.v01.Coord;
+import org.matsim.api.core.v01.network.Network;
+import org.matsim.core.network.NetworkUtils;
+import org.matsim.core.utils.geometry.CoordUtils;
+import org.matsim.core.utils.geometry.CoordinateTransformation;
+import org.matsim.core.utils.geometry.transformations.TransformationFactory;
+
+import java.util.*;
+import java.util.stream.Collectors;
+
+/**
+ * @author nkuehnel / MOIA
+ */
+public class H3GridUtils {
+
+	static final Logger log = LogManager.getLogger(H3GridUtils.class);
+
+	public static Map<String, PreparedGeometry> createH3GridFromNetwork(Network network, int resolution, String crs) {
+
+		H3Core h3 = H3Utils.getInstance();
+
+		log.info("start creating H3 grid from network at resolution " + resolution);
+		double hexagonEdgeLengthAvg = h3.getHexagonEdgeLengthAvg(resolution, LengthUnit.m);
+		log.info("Average edge length: " + hexagonEdgeLengthAvg + " meters.");
+		log.info("Average centroid distance: " + hexagonEdgeLengthAvg * Math.sqrt(3) + " meters.");
+		log.info("Average hexagon area: " + h3.getHexagonAreaAvg(resolution, AreaUnit.m2) + " m^2");
+
+		double[] boundingbox = NetworkUtils.getBoundingBox(network.getNodes().values());
+		double minX = boundingbox[0];
+		double maxX = boundingbox[2];
+		double minY = boundingbox[1];
+		double maxY = boundingbox[3];
+
+		GeometryFactory gf = new GeometryFactory();
+		PreparedGeometryFactory preparedGeometryFactory = new PreparedGeometryFactory();
+		Map<String, PreparedGeometry> grid = new HashMap<>();
+		CoordinateTransformation toLatLong = TransformationFactory.getCoordinateTransformation(crs, TransformationFactory.WGS84);
+		CoordinateTransformation fromLatLong = TransformationFactory.getCoordinateTransformation(TransformationFactory.WGS84, crs);
+
+		List<LatLng> boundingBoxPoints = new ArrayList<>();
+
+		Coord bottomLeft = toLatLong.transform(new Coord(minX, minY));
+		Coord topLeft = toLatLong.transform(new Coord(minX, maxY));
+		Coord topRight = toLatLong.transform(new Coord(maxX, maxY));
+		Coord bottomRight = toLatLong.transform(new Coord(maxX, minY));
+
+		boundingBoxPoints.add(coordToLatLng(bottomLeft));
+		boundingBoxPoints.add(coordToLatLng(topLeft));
+		boundingBoxPoints.add(coordToLatLng(topRight));
+		boundingBoxPoints.add(coordToLatLng(bottomRight));
+		boundingBoxPoints.add(coordToLatLng(bottomLeft));
+
+		long millis = System.currentTimeMillis();
+
+		//get cells in a finer resolution to catch links at the border
+		List<String> h3Grid = h3.polygonToCellAddresses(boundingBoxPoints, Collections.emptyList(), Math.min(H3Utils.MAX_RES, resolution));
+		h3Grid = h3Grid
+			.parallelStream()
+			//also include neighbors with distance 1
+			.flatMap(h3Id -> h3.gridDisk(h3Id, 1).stream())
+			.distinct()
+			.toList();
+
+		if(h3Grid.isEmpty()) {
+			// bounding box too small to cover even a single H3 cell for a significant part. Use bounding box coords directly.
+			h3Grid = boundingBoxPoints.stream().map(corner -> h3.latLngToCellAddress(corner.lat, corner.lng, resolution)).distinct().toList();
+		}
+
+		log.info("Obtained " + h3Grid.size() + " H3 cells in " + (System.currentTimeMillis() - millis) + " ms.");
+
+
+		for (String h3Id : h3Grid) {
+			List<Coordinate> coordinateList = h3.cellToBoundary(h3Id)
+				.stream()
+				.map(latLng -> CoordUtils.createGeotoolsCoordinate(fromLatLong.transform(latLngToCoord(latLng))))
+				.collect(Collectors.toList());
+
+			if (!coordinateList.isEmpty()) {
+				coordinateList.add(coordinateList.get(0));
+			}
+
+			Polygon polygon = new Polygon(gf.createLinearRing(coordinateList.toArray(new Coordinate[0])), null, gf);
+			grid.put(h3Id, preparedGeometryFactory.create(polygon));
+		}
+
+		log.info("finished creating H3 grid from network.");
+		return grid;
+	}
+
+	public static LatLng coordToLatLng(Coord coord) {
+		//invert coordinate order
+		return new LatLng(coord.getY(), coord.getX());
+	}
+
+	public static Coord latLngToCoord(LatLng latLng) {
+		//invert coordinate order
+		return new Coord(latLng.lng, latLng.lat);
+	}
+}

contribs/common/src/main/java/org/matsim/contrib/common/zones/h3/H3Utils.java

@@ -0,0 +1,27 @@
+package org.matsim.contrib.common.zones.h3;
+
+import com.uber.h3core.H3Core;
+
+import java.io.IOException;
+
+/**
+ * @author nkuehnel / MOIA
+ */
+public final class H3Utils {
+
+	private static H3Core h3;
+
+	public final static int MAX_RES = 16;
+
+
+	public static H3Core getInstance() {
+		if(h3 == null) {
+			try {
+				h3 = H3Core.newInstance();
+			} catch (IOException e) {
+				throw new RuntimeException(e);
+			}
+		}
+		return h3;
+	}
+}