Skip to content

Commit

Permalink
Merge pull request #3143 from moia-oss/gpkg-support
Browse files Browse the repository at this point in the history
proposal: add geopackage support to matsim gis
  • Loading branch information
nkuehnel authored Mar 7, 2024
2 parents f5ba881 + 3f11a92 commit 7ff5887
Show file tree
Hide file tree
Showing 72 changed files with 1,352 additions and 1,001 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@
import org.matsim.core.network.NetworkUtils;
import org.matsim.core.network.algorithms.TransportModeNetworkFilter;
import org.matsim.core.utils.geometry.CoordUtils;
import org.matsim.core.utils.gis.ShapeFileReader;
import org.matsim.core.utils.gis.GeoFileReader;
import org.matsim.facilities.*;
import org.opengis.feature.simple.SimpleFeature;

Expand Down Expand Up @@ -225,7 +225,7 @@ public static final ActivityFacilities createFacilityForEachLink(String facility


public static final ActivityFacilities createFacilityFromBuildingShapefile(String shapeFileName, String identifierCaption, String numberOfHouseholdsCaption) {
ShapeFileReader shapeFileReader = new ShapeFileReader();
GeoFileReader shapeFileReader = new GeoFileReader();
Collection<SimpleFeature> features = shapeFileReader.readFileAndInitialize(shapeFileName);

ActivityFacilities facilities = FacilitiesUtils.createActivityFacilities("DensitiyFacilities");
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -28,17 +28,17 @@
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.triangulate.VoronoiDiagramBuilder;
import org.matsim.contrib.matrixbasedptrouter.utils.BoundingBox;
import org.matsim.core.utils.gis.ShapeFileWriter;
import org.matsim.core.utils.gis.GeoFileWriter;
import org.opengis.feature.simple.SimpleFeature;

/**
* @author dziemke
*/
class VoronoiExample {

public static void main(String[] args) {
GeometryFactory geometryFactory = new GeometryFactory();

Collection<Coordinate> sites = new ArrayList<>();
sites.add(new Coordinate(70, 70));
sites.add(new Coordinate(50, 150));
Expand All @@ -48,16 +48,16 @@ public static void main(String[] args) {
sites.add(new Coordinate(250, 150));
sites.add(new Coordinate(350, 50));
sites.add(new Coordinate(370, 170));

VoronoiDiagramBuilder voronoiDiagramBuilder = new VoronoiDiagramBuilder();
voronoiDiagramBuilder.setSites(sites);
voronoiDiagramBuilder.setSites(sites);

List<Polygon> polygons = voronoiDiagramBuilder.getSubdivision().getVoronoiCellPolygons(geometryFactory);

BoundingBox boundingBox = BoundingBox.createBoundingBox(0, 0, 400, 200);
Polygon boundingPolygon = VoronoiGeometryUtils.createBoundingPolygon(boundingBox);
Collection<Geometry> cutGeometries = VoronoiGeometryUtils.cutPolygonsByBoundary(polygons, boundingPolygon);
Collection<SimpleFeature> features = VoronoiGeometryUtils.createFeaturesFromPolygons(cutGeometries);
ShapeFileWriter.writeGeometries(features, "/Users/dominik/voronoi_test.shp");
GeoFileWriter.writeGeometries(features, "/Users/dominik/voronoi_test.shp");
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@
import org.matsim.core.utils.geometry.CoordinateTransformation;
import org.matsim.core.utils.geometry.geotools.MGC;
import org.matsim.core.utils.geometry.transformations.TransformationFactory;
import org.matsim.core.utils.gis.ShapeFileWriter;
import org.matsim.core.utils.gis.GeoFileWriter;
import org.matsim.facilities.ActivityFacility;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
Expand All @@ -74,7 +74,7 @@ public GeoserverUpdater (String crs, String name, Map<Id<ActivityFacility>, Geom
this.pushing2Geoserver = pushing2Geoserver;
this.createQGisOutput = createQGisOutput;
}

private Map<Tuple<ActivityFacility, Double>, Map<String,Double>> accessibilitiesMap = new HashMap<>() ;

@Override
Expand All @@ -95,16 +95,16 @@ public void finish() {
SimpleFeatureTypeBuilder featureTypeBuilder = createFeatureTypeBuilder();
SimpleFeatureType featureType = featureTypeBuilder.buildFeatureType();
DefaultFeatureCollection featureCollection = createFeatureCollection(geometryFactory, featureType);

if (outputDirectory != null) {
File file = new File(outputDirectory);
file.mkdirs();
}

if (createQGisOutput) {
ShapeFileWriter.writeGeometries(featureCollection, outputDirectory + "/result.shp");
GeoFileWriter.writeGeometries(featureCollection, outputDirectory + "/result.shp");
}

if (pushing2Geoserver) {
updateOnGeoserver(featureType, featureCollection);
}
Expand All @@ -130,7 +130,7 @@ private DefaultFeatureCollection createFeatureCollection(GeometryFactory geometr
LOG.info("Start creating features from accessibility data.");
DefaultFeatureCollection featureCollection = new DefaultFeatureCollection("internal", featureType);
SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(featureType);

CoordinateTransformation transformation = TransformationFactory.getCoordinateTransformation(this.crs, TransformationFactory.WGS84);

for (Entry<Tuple<ActivityFacility, Double>, Map<String, Double>> entry : accessibilitiesMap.entrySet()) {
Expand All @@ -145,10 +145,10 @@ private DefaultFeatureCollection createFeatureCollection(GeometryFactory geometr
i++;
}
featureBuilder.add(geometryFactory.createPolygon(transformedCoordinates));

featureBuilder.add(Integer.parseInt(entry.getKey().getFirst().getId().toString()));
featureBuilder.add(entry.getKey().getSecond());

for (Modes4Accessibility modeEnum : Modes4Accessibility.values()) {
String mode = modeEnum.toString();
Double accessibility = entry.getValue().get(mode);
Expand Down Expand Up @@ -188,14 +188,14 @@ private void updateOnGeoserver(SimpleFeatureType featureType, DefaultFeatureColl
// There have been errors with the data store if the dependency "gt-jdbc-postgis", version 13.0 was missing!
DataStore dataStore = DataStoreFinder.getDataStore(params);
LOG.info("dataStore = " + dataStore);

// Remove schema in case it already exists
try {
dataStore.removeSchema(name);
} catch (IllegalArgumentException e) {
LOG.warn("Could not remove schema. Probably, it has not existed yet.");
}

dataStore.createSchema(featureType);
SimpleFeatureStore featureStore = (SimpleFeatureStore) dataStore.getFeatureSource(name);
featureStore.addFeatures(featureCollection);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@
import org.matsim.core.utils.geometry.CoordinateTransformation;
import org.matsim.core.utils.geometry.geotools.MGC;
import org.matsim.core.utils.geometry.transformations.TransformationFactory;
import org.matsim.core.utils.gis.ShapeFileReader;
import org.matsim.core.utils.gis.GeoFileReader;
import org.opengis.feature.simple.SimpleFeature;

/**
Expand All @@ -59,26 +59,26 @@ public class AccidentsNetworkModification {
public AccidentsNetworkModification(Scenario scenario) {
this.scenario = scenario;
}

public Network setLinkAttributsBasedOnOSMFile(String landuseOsmFile, String osmCRS, String[] tunnelLinkIDs, String[] planfreeLinkIDs) throws MalformedURLException, IOException {

AccidentsConfigGroup accidentsCfg = (AccidentsConfigGroup) scenario.getConfig().getModules().get(AccidentsConfigGroup.GROUP_NAME);

Map<String, SimpleFeature> landUseFeaturesBB = new HashMap<>();
Map<String, String> landUseDataBB = new HashMap<>();


log.info("Initializing all link-specific information...");

if (landuseOsmFile == null) {
log.warn("Landuse shape file is null. Using default values...");
} else {
SimpleFeatureSource ftsLandUseBB;
if (!landuseOsmFile.startsWith("http")) {
ftsLandUseBB = ShapeFileReader.readDataFile(landuseOsmFile);
ftsLandUseBB = GeoFileReader.readDataFile(landuseOsmFile);
} else {
ftsLandUseBB = FileDataStoreFinder.getDataStore(new URL(landuseOsmFile)).getFeatureSource();
}
}
try (SimpleFeatureIterator itLandUseBB = ftsLandUseBB.getFeatures().features()) {
while (itLandUseBB.hasNext()) {
SimpleFeature ftLandUseBB = itLandUseBB.next();
Expand All @@ -95,57 +95,57 @@ public Network setLinkAttributsBasedOnOSMFile(String landuseOsmFile, String osmC
e.printStackTrace();
}
}


int linkCounter = 0;
for (Link link : this.scenario.getNetwork().getLinks().values()) {

if (linkCounter % 100 == 0) {
log.info("Link #" + linkCounter + " (" + (int) ((double) linkCounter / this.scenario.getNetwork().getLinks().size() * 100) + "%)");
}
linkCounter++;
linkCounter++;

link.getAttributes().putAttribute(accidentsCfg.getAccidentsComputationMethodAttributeName(), AccidentsConfigGroup.AccidentsComputationMethod.BVWP.toString());

ArrayList<Integer> bvwpRoadType = new ArrayList<>();

// 'plangleich', 'planfrei' or tunnel?
bvwpRoadType.add(0, 1);

for(int j=0; j < planfreeLinkIDs.length; j++){
if(planfreeLinkIDs[j].equals(String.valueOf(link.getId()))){
bvwpRoadType.set(0, 0); // Change to Plan free
log.info(link.getId() + " Changed to Plan free!");
break;
}
}

for(int i=0; i < tunnelLinkIDs.length; i++){
if(tunnelLinkIDs[i].equals(String.valueOf(link.getId()))){
bvwpRoadType.set(0, 2); // Change to Tunnel
log.info(link.getId() + " Changed to Tunnel!");
break;
}
}

// builtup or not builtup area?
String osmLandUseFeatureBBId = getOSMLandUseFeatureBBId(link, landUseFeaturesBB, osmCRS);
if (osmLandUseFeatureBBId == null) {

if (osmLandUseFeatureBBId == null) {
log.warn("No area type found for link " + link.getId() + ". Using default value: not built-up area.");
if (link.getFreespeed() > 16.) {
bvwpRoadType.add(1, 0);
} else {
bvwpRoadType.add(1, 2);
}

} else {
String landUseTypeBB = landUseDataBB.get(osmLandUseFeatureBBId);
if (landUseTypeBB.matches("commercial|industrial|recreation_ground|residential|retail")) { //built-up area
if (link.getFreespeed() > 16.) {
bvwpRoadType.add(1, 1);
} else {
bvwpRoadType.add(1, 3);
bvwpRoadType.add(1, 3);
}
} else {
if (link.getFreespeed() > 16.) {
Expand All @@ -163,85 +163,85 @@ public Network setLinkAttributsBasedOnOSMFile(String landuseOsmFile, String osmC
numberOfLanesBVWP = (int) link.getNumberOfLanes();
}
bvwpRoadType.add(2, numberOfLanesBVWP);

link.getAttributes().putAttribute( AccidentsConfigGroup.BVWP_ROAD_TYPE_ATTRIBUTE_NAME, bvwpRoadType.get(0) + "," + bvwpRoadType.get(1) + "," + bvwpRoadType.get(2));
}
log.info("Initializing all link-specific information... Done.");
return scenario.getNetwork();
}

private String getOSMLandUseFeatureBBId(Link link, Map<String, SimpleFeature> landUseFeaturesBB, String osmCRS) {

if (landUseFeaturesBB == null || landUseFeaturesBB.isEmpty()) return null;

CoordinateTransformation ctScenarioCRS2osmCRS = TransformationFactory.getCoordinateTransformation(this.scenario.getConfig().global().getCoordinateSystem(), osmCRS);

Coord linkCoordinateTransformedToOSMCRS = ctScenarioCRS2osmCRS.transform(link.getCoord()); // this Method gives the middle point of the link back
Point pMiddle = MGC.xy2Point(linkCoordinateTransformedToOSMCRS.getX(), linkCoordinateTransformedToOSMCRS.getY());

Coord linkStartCoordinateTransformedToOSMCRS = ctScenarioCRS2osmCRS.transform(link.getFromNode().getCoord());
Point pStart = MGC.xy2Point(linkStartCoordinateTransformedToOSMCRS.getX(), linkStartCoordinateTransformedToOSMCRS.getY());

Coord linkEndCoordinateTransformedToOSMCRS = ctScenarioCRS2osmCRS.transform(link.getToNode().getCoord());
Point pEnd = MGC.xy2Point(linkEndCoordinateTransformedToOSMCRS.getX(), linkEndCoordinateTransformedToOSMCRS.getY());

String osmLandUseFeatureBBId = null;

for (SimpleFeature feature : landUseFeaturesBB.values()) {
if (((Geometry) feature.getDefaultGeometry()).contains(pMiddle)) {
return osmLandUseFeatureBBId = feature.getAttribute("osm_id").toString();
}
}

for (SimpleFeature feature : landUseFeaturesBB.values()) {
if (((Geometry) feature.getDefaultGeometry()).contains(pStart)) {
return osmLandUseFeatureBBId = feature.getAttribute("osm_id").toString();
}
}

for (SimpleFeature feature : landUseFeaturesBB.values()) {
if (((Geometry) feature.getDefaultGeometry()).contains(pEnd)) {
return osmLandUseFeatureBBId = feature.getAttribute("osm_id").toString();
}
}

// look around the link

GeometryFactory geoFac = new GeometryFactory();
CoordinateTransformation cTosmCRSToGK4 = TransformationFactory.getCoordinateTransformation(osmCRS, "EPSG:31468");
double distance = 10.0;

double distance = 10.0;

while (osmLandUseFeatureBBId == null && distance <= 500) {
Coord coordGK4 = cTosmCRSToGK4.transform(MGC.coordinate2Coord(pMiddle.getCoordinate()));
Point pGK4 = geoFac.createPoint(MGC.coord2Coordinate(coordGK4));

Point pRightGK4 = geoFac.createPoint(new Coordinate(pGK4.getX() + distance, pGK4.getY()));
Point pRight = transformPointFromGK4ToOSMCRS(pRightGK4, osmCRS);

Point pDownGK4 = geoFac.createPoint(new Coordinate(pGK4.getX(), pGK4.getY() - distance));
Point pDown = transformPointFromGK4ToOSMCRS(pDownGK4, osmCRS);

Point pLeftGK4 = geoFac.createPoint(new Coordinate(pGK4.getX() - distance, pGK4.getY()));
Point pLeft = transformPointFromGK4ToOSMCRS(pLeftGK4, osmCRS);

Point pUpGK4 = geoFac.createPoint(new Coordinate(pGK4.getX(), pGK4.getY() + distance));
Point pUp = transformPointFromGK4ToOSMCRS(pUpGK4, osmCRS);

Point pUpRightGK4 = geoFac.createPoint(new Coordinate(pGK4.getX() + distance, pGK4.getY() + distance));
Point pUpRight = transformPointFromGK4ToOSMCRS(pUpRightGK4, osmCRS);

Point pDownRightGK4 = geoFac.createPoint(new Coordinate(pGK4.getX() + distance, pGK4.getY() - distance));
Point pDownRight = transformPointFromGK4ToOSMCRS(pDownRightGK4, osmCRS);

Point pDownLeftGK4 = geoFac.createPoint(new Coordinate(pGK4.getX() - distance, pGK4.getY() - distance));
Point pDownLeft = transformPointFromGK4ToOSMCRS(pDownLeftGK4, osmCRS);

Point pUpLeftGK4 = geoFac.createPoint(new Coordinate(pGK4.getX() - distance, pGK4.getY() + distance));
Point pUpLeft = transformPointFromGK4ToOSMCRS(pUpLeftGK4, osmCRS);

for (SimpleFeature feature : landUseFeaturesBB.values()) {

if (((Geometry) feature.getDefaultGeometry()).contains(pRight)) {
osmLandUseFeatureBBId = feature.getAttribute("osm_id").toString();
return osmLandUseFeatureBBId;
Expand All @@ -268,17 +268,17 @@ private String getOSMLandUseFeatureBBId(Link link, Map<String, SimpleFeature> la
return osmLandUseFeatureBBId;
}
}

distance += 10.0;
}

log.warn("No area type found. Returning null...");
return null;
}

private Point transformPointFromGK4ToOSMCRS(Point pointGK4, String osmCRS) {
CoordinateTransformation ctGK4toOSMCRS = TransformationFactory.getCoordinateTransformation("EPSG:31468", osmCRS);

Coord coordGK4 = MGC.coordinate2Coord(pointGK4.getCoordinate());
Point pointOSMCRS = new GeometryFactory().createPoint(MGC.coord2Coordinate(ctGK4toOSMCRS.transform(coordGK4)));
return pointOSMCRS;
Expand Down
Loading

0 comments on commit 7ff5887

Please sign in to comment.