Road matching data

1
import pandas as pd
2
import pyarrow as pa
3
import pyarrow.parquet as pq
4

5
roads = pd.DataFrame({
6
    "id": [5707295, 5707296, 5707297],
7
    "geometry": [
8
        "LINESTRING (4.8952 52.3702, 4.8960 52.3710)",
9
        "LINESTRING (4.8960 52.3710, 4.8975 52.3725)",
10
        "LINESTRING (4.8975 52.3725, 4.8990 52.3740)",
11
    ],
12
    "is_navigable": [True, True, False],
13
})
14

15
schema = pa.schema([("id", pa.int64()), ("geometry", pa.string()), ("is_navigable", pa.bool_())])
16
pq.write_table(pa.Table.from_pandas(roads, schema=schema, preserve_index=False), "roads.parquet")

1
import geopandas as gpd
2
import pandas as pd
3

4
roads = gpd.read_file("roads.shp").to_crs("EPSG:4326").explode(index_parts=False)
5
roads = roads[roads.geom_type == "LineString"].reset_index(drop=True)
6

7
part = roads.groupby("road_id").cumcount()
8
road_id = roads["road_id"].astype(str)
9

10
pd.DataFrame({
11
    "id": road_id.where(part == 0, road_id + "-" + part.astype(str)),
12
    "geometry": roads.geometry.to_wkt(),
13
}).to_parquet("roads.parquet", index=False)

1
from pyspark.sql import SparkSession, functions as F
2

3
spark = SparkSession.builder.appName("gem-road-input").getOrCreate()
4

5
roads = spark.read.parquet("source_roads").select(
6
    F.col("road_id").alias("id"),
7
    F.col("wkt").alias("geometry"),
8
    F.col("navigable").cast("boolean").alias("is_navigable"),
9
)
10
roads.coalesce(1).write.mode("overwrite").parquet("roads_spark")

1
az storage blob upload --account-name "<STORAGE_ACCOUNT_NAME>" \
2
                       --container-name "default" \
3
                       --name "roads.parquet" \
4
                       --file "$(ls roads_spark/part-*.parquet)" \
5
                       --auth-mode login

1
pip install pyarrow shapely
2
python validate_road_input.py roads.parquet

1
"""Check a road matching input file before uploading it to GEM.
2

3
Usage: python validate_road_input.py <file.parquet>
4
"""
5
import sys
6
from pathlib import Path
7

8
import numpy as np
9
import pyarrow as pa
10
import pyarrow.compute as pc
11
import pyarrow.parquet as pq
12
import shapely
13

14
REQUIRED_COLUMNS = ["id", "geometry"]
15
EXPECTED_TYPES = {
16
    "id": ([pa.int32(), pa.int64(), pa.string(), pa.large_string()], "'integer' or 'string'"),
17
    "geometry": ([pa.string(), pa.large_string()], "'string'"),
18
    "is_navigable": ([pa.bool_()], "'boolean'"),
19
}
20
TYPE_NAMES = {"int8": "tinyint", "int16": "smallint", "int32": "int", "int64": "bigint",
21
              "bool": "boolean", "large_string": "string", "large_binary": "binary"}
22
LINE_TYPES = [shapely.GeometryType.LINESTRING, shapely.GeometryType.LINEARRING]
23

24

25
def quoted(names):
26
    return ", ".join(f"'{name}'" for name in names)
27

28

29
def schema_errors(schema, name):
30
    errors = []
31
    missing = [column for column in REQUIRED_COLUMNS if column not in schema.names]
32
    if missing:
33
        errors.append(f"Input file '{name}' is missing required column(s): {quoted(missing)}. "
34
                      f"Required columns: {quoted(REQUIRED_COLUMNS)}.")
35
    for column, (allowed, expected) in EXPECTED_TYPES.items():
36
        if column in schema.names and schema.field(column).type not in allowed:
37
            actual = str(schema.field(column).type)
38
            errors.append(f"Column '{column}' in input file '{name}' has unsupported type "
39
                          f"'{TYPE_NAMES.get(actual, actual)}'; expected {expected}.")
40
    return errors
41

42

43
def data_findings(path, name):
44
    errors, warnings = [], []
45
    ids = pq.read_table(path, columns=["id"]).column("id").combine_chunks()
46
    duplicated = pc.sum(pc.greater(pc.value_counts(ids).field("counts"), 1)).as_py() or 0
47
    if duplicated:
48
        errors.append(f"Column 'id' in input file '{name}' contains {duplicated} duplicated value(s). "
49
                      "Each road must have a unique 'id'.")
50
    if ids.null_count:
51
        warnings.append(f"Column 'id' has {ids.null_count} empty value(s). Give every road an ID.")
52

53
    invalid = empty = outside = 0
54
    for batch in pq.ParquetFile(path).iter_batches(batch_size=100_000, columns=["geometry"]):
55
        geometries = shapely.from_wkt(batch.column(0).to_numpy(zero_copy_only=False), on_invalid="ignore")
56
        lines = np.isin(shapely.get_type_id(geometries), LINE_TYPES)
57
        invalid += int((~lines).sum())
58
        empty += int((lines & shapely.is_empty(geometries)).sum())
59
        bounds = shapely.bounds(geometries[lines])
60
        outside += int(((np.abs(bounds[:, [0, 2]]) > 180) | (np.abs(bounds[:, [1, 3]]) > 90)).any(axis=1).sum())
61
    if invalid:
62
        errors.append(f"Column 'geometry' in input file '{name}' contains {invalid} value(s) "
63
                      "that are not valid WKT LINESTRING geometries.")
64
    if empty:
65
        warnings.append(f"Column 'geometry' has {empty} empty LINESTRING value(s), which can't be matched.")
66
    if outside:
67
        warnings.append(f"Column 'geometry' has {outside} value(s) with coordinates outside the "
68
                        "longitude/latitude range. Use WGS84 (EPSG:4326) with longitude first.")
69
    return errors, warnings
70

71

72
def main(file_name):
73
    path = Path(file_name)
74
    if not path.is_file():
75
        print(f"ERROR: {path} does not exist.")
76
        return 1
77
    try:
78
        schema = pq.read_schema(path)
79
    except pa.ArrowInvalid:
80
        print(f"ERROR: Input file '{path.name}' is not a valid Parquet file.")
81
        return 1
82

83
    errors, warnings = schema_errors(schema, path.name), []
84
    if not errors:
85
        errors, warnings = data_findings(path, path.name)
86
    if "is_navigable" not in schema.names:
87
        warnings.append("Column 'is_navigable' is missing, so GEM treats all roads as navigable.")
88

89
    for message in errors:
90
        print(f"ERROR: {message}")
91
    for message in warnings:
92
        print(f"WARNING: {message}")
93
    if not errors:
94
        print(f"OK: {path.name} is ready for road matching.")
95
    return 1 if errors else 0
96

97

98
if __name__ == "__main__":
99
    if len(sys.argv) != 2:
100
        sys.exit(__doc__)
101
    sys.exit(main(sys.argv[1]))

1
ERROR: Column 'id' in input file 'roads.parquet' contains 2 duplicated value(s). Each road must have a unique 'id'.
2
WARNING: Column 'geometry' has 15 value(s) with coordinates outside the longitude/latitude range. Use WGS84 (EPSG:4326) with longitude first.

1
Your road 5707296   0%              40%                          100%
2
                    |===============|============================|
3
                    |---------------|                              lr_id 0.0-40.0#5f1a9c3b-…
4
                                    |----------------------------| lr_id 40.0-100.0#9d2e71a3-…

1
import pandas as pd
2

3
results = pd.read_parquet("roads.results.parquet")
4
for column in ("lr_id", "lr_gers"):
5
    parts = results[column].str.extract(r"^(?P<start>[\d.]+)-(?P<end>[\d.]+)#(?P<ref>.+)$")
6
    results[f"{column}_start"] = parts["start"].astype(float)
7
    results[f"{column}_end"] = parts["end"].astype(float)
8

9
print(results[["id", "gers", "confidence", "lr_id_start", "lr_id_end"]])

1
az storage blob download-batch --account-name "<STORAGE_ACCOUNT_NAME>" \
2
                               --source "default" \
3
                               --pattern "<name>-visualization-data/*" \
4
                               --destination "/path/to/destination" \
5
                               --auth-mode login