GIScience / osm-transform

Custom data filtering and enrichment for OpenStreetMap PBF files for a better openrouteservice graph building
3 stars 0 forks source link

Create a strategy to cut road features based on elevation cell sizes #19

Closed TheGreatRefrigerator closed 6 months ago

TheGreatRefrigerator commented 9 months ago

Test ors route: current ors: https://classic-maps.openrouteservice.org/directions?n1=49.75001&n2=8.850356&n3=16&a=49.743549,8.849766,49.751008,8.848221&b=1a&c=0&k1=en-US&k2=km

GeoJSON result: ors-export-linestring.json

Optimal solution (QGIS elevation profile for route): optimal 1m resolution in red ors output violet

grafik
rtroilo commented 9 months ago

Uninterpolated values with srtm and dgm elevation data

(Lon,Lat,Ele); Pixel; DGM Data (elevation)
(8.849782, 49.743523, 197.1); [177,874]; [195.58299255371094]
(8.851082, 49.74387, 196.7); [270,835]; [194.5989990234375]
(8.851668, 49.744073, 196.2); [313,813]; [194.8699951171875]
(8.851882, 49.744145, 196.1); [328,805]; [194.82899475097656]
(8.852762, 49.744462, 195.3); [392,770]; [196.11399841308594]
(8.853751, 49.744773, 194.2); [463,735]; [194.3990020751953]
(8.854025, 49.744864, 193.9); [483,725]; [194.0959930419922]
(8.853816, 49.745385, 194.8); [468,667]; [200.83399963378906]
(8.854237, 49.745284, 195.6); [498,678]; [198.58399963378906]
(8.854402, 49.745302, 195.8); [510,676]; [197.72500610351562]
(8.854549, 49.745374, 196.0); [521,669]; [197.08299255371094]
(8.854947, 49.74586, 197.4); [549,615]; [193.08099365234375]
(8.854356, 49.746088, 198.6); [507,589]; [191.7969970703125]
(8.853302, 49.746502, 199.9); [431,543]; [191.0590057373047]
(8.853299, 49.746525, 200.3); [431,540]; [191.03900146484375]
(8.853234, 49.747193, 201.2); [426,466]; [193.2899932861328]
(8.852727, 49.747171, 201.2); [390,468]; [194.43800354003906]
(8.852308, 49.747095, 201.2); [359,477]; [195.6540069580078]
(8.852244, 49.747526, 201.5); [355,429]; [201.61099243164062]
(8.852169, 49.747884, 201.6); [350,389]; [207.62899780273438]
(8.852345, 49.747938, 201.7); [362,383]; [208.08099365234375]
(8.852324, 49.748482, 205.3); [361,323]; [214.37100219726562]
(8.852295, 49.748764, 208.0); [359,291]; [217.1909942626953]
(8.852209, 49.749014, 209.7); [353,263]; [219.6199951171875]
(8.852182, 49.749293, 212.5); [351,232]; [219.00799560546875]
(8.851084, 49.749328, 219.1); [272,228]; [226.49000549316406]
(8.850858, 49.749367, 219.8); [255,224]; [227.53199768066406]
(8.850641, 49.74946, 220.5); [240,214]; [228.55599975585938]
(8.850473, 49.749574, 221.2); [228,201]; [229.6280059814453]
(8.848634, 49.750824, 231.5); [96,62]; [249.322998046875]
(8.848294, 49.751053, 233.3); [71,36]; [251.2030029296875]
rtroilo commented 9 months ago

javaish/pseudo code for cutting strategy. Thank you @TheGreatRefrigerator @koebi

var wayNodeLocations = nodeLocationOf(way);
for (var i = 1; i < wayNodeLocations.size(); i++) {
    var from = wayNodeLocations.get(i - 1);
    newWayRefs.add(from);
    var wgs84Coords = interpolatedCoordinates(wgs84StepWidth, from, wayNodeLocations.get(i));
    var gridCoords = wgs84ToGrid.apply(wgs84Coords);
    var elevation = elevationData(raster, gridCoords);
    for (var newNodeLocation : filterByTreshold(threshold, wgs84Coords, elevation)) {
        var nodeId = generateNodeId(newNodeLocation);
        writeNode(nodeId, newNodeLocation);
        newWayRefs.add(nodeId);
    }
}
newWayRefs.add(wayNodeLocations.getLast());

double[] interpolatedCoordinates(double step, Coordinate from, Coordinate to) {
    var deltaX = to.x - from.x;
    var deltaY = to.y - from.y;
    var length = sqrt(deltaX * deltaX + deltaY * deltaY);

    var nX = deltaX / length;
    var nY = deltaY / length;
    var sX = nX * step;
    var sY = nY * step;

    var steps = (int) (deltaX / sX);

    var x = from.x;
    var y = from.y;
    var coords = new double[(2 + steps) * 2];
    for (var s = 0; s <= steps; s++) {
        coords[s * 2] = x + sX * s;
        coords[s * 2 + 1] = y + sY * s;
    }
    coords[coords.length - 2] = to.x;
    coords[coords.length - 1] = to.y;
    return coords;
}

Coordinates filterByThreshold(double threshold, double[] split, double[] elevation) {
    var list = new List();
    var ele = elevation[0];
    for (var i = 1; i < elevation.length - 1; i++) {
        if (abs(elevation[i] - ele) > threshold || abs(elevation[i + 1] - ele) > threshold) {
            ele = elevation[i];
            list.add(new Coordinate(split[i * 2], split[i * 2 + 1], ele));
        }
    }
    return list;
}