Skip to content
Snippets Groups Projects
test_gis_calculate.py 10.1 KiB
Newer Older
  • Learn to ignore specific revisions
  • Pedro L. Magalhães's avatar
    Pedro L. Magalhães committed
    # imports
    
    # standard
    
    from math import isclose, inf
    
    # local, external
    
    import networkx as nx
    import osmnx as ox
    from shapely.geometry import LineString
    from numpy import nan
    
    # local, internal
    
    import src.topupopt.data.gis.calculate as gis_calc
    import src.topupopt.data.gis.osm as osm
    
    # *****************************************************************************
    # *****************************************************************************
    
    
    class TestGisCalculate:
        # *************************************************************************
        # *************************************************************************
    
        def validate_edge_distances(self, G: nx.MultiDiGraph, abs_tol: float = 5):
            # get the true edge lengths
            true_lengths = {
                edge_key: (
                    G.edges[edge_key][osm.KEY_OSMNX_LENGTH]
                    if osm.KEY_OSMNX_LENGTH in G.edges[edge_key]
                    else None
                )
                for edge_key in G.edges(keys=True)
            }
    
            # get the edge lengths calculated independently
            calculated_lengths = gis_calc.edge_lengths(G)
    
            # for each edge on the graph
            for edge_key in true_lengths.keys():
                # validate
                assert isclose(
                    calculated_lengths[edge_key], true_lengths[edge_key], abs_tol=abs_tol
                )
    
        # *************************************************************************
        # *************************************************************************
    
        def test_distances(self):
            # get a graph
    
            G = ox.graph_from_point(
                (55.71654, 9.11728),
                network_type="drive",
                custom_filter='["highway"~"residential|tertiary|unclassified|service"]',
                truncate_by_edge=True,
            )
    
            # count occurrences
    
            self.example_count_occurrences(G)
    
            # validate without projected coordinates
    
            self.validate_edge_distances(G=G)
    
            # project the graph
    
            projected_G = ox.project_graph(G=G)
    
            # validate with projected coordinates
    
            self.validate_edge_distances(G=projected_G)
    
            # calculate node path lenths
    
            self.example_node_path_lengths()
    
            # calculate edge path lengths
    
            self.example_edge_path_lengths()
    
        # *************************************************************************
        # *************************************************************************
    
        def example_node_path_lengths(self):
            # identify path between two nodes
    
            G = nx.MultiDiGraph()
    
            G.add_edge(0, 1, length=10)
    
            G.add_edge(1, 2, length=10)
    
            G.add_edge(1, 2, length=5)
    
            G.add_edge(2, 3, length=5)
    
            # calculate length
    
            path = [0, 1, 2, 3]
    
            this_length = gis_calc.node_path_length(network=G, path=path)
    
            assert this_length == 20
    
            # calculate all possible lengths
    
            these_lengths = gis_calc.node_path_length(
                network=G, path=path, return_minimum_length_only=False
            )
    
            true_lengths = [25, 20]
    
            assert len(these_lengths) == len(true_lengths)
    
            for _length in these_lengths:
                assert _length in true_lengths
    
            # *********************************************************************
            # *********************************************************************
    
            # no path
    
            path = [0, 3]
    
            this_length = gis_calc.node_path_length(network=G, path=path)
    
            assert this_length == inf
    
            these_lengths = gis_calc.node_path_length(
                network=G, path=path, return_minimum_length_only=False
            )
    
            assert these_lengths[0] == inf
    
            # empty path
    
            path = []
    
            this_length = gis_calc.node_path_length(network=G, path=path)
    
            assert this_length == inf
    
            # *********************************************************************
            # *********************************************************************
    
        # *************************************************************************
        # *************************************************************************
    
        def example_edge_path_lengths(self):
            # identify path between two nodes
    
            G = nx.MultiDiGraph()
    
            G.add_edge(0, 1, length=10)
    
            G.add_edge(1, 2, length=10)
    
            G.add_edge(1, 2, length=5)
    
            G.add_edge(2, 3, length=5)
    
            # calculate length
    
            path = [(0, 1, 0), (1, 2, 0), (2, 3, 0)]
    
            this_length = gis_calc.edge_path_length(network=G, path=path)
    
            assert this_length == 25
    
            # alternative path
    
            path = [(0, 1, 0), (1, 2, 1), (2, 3, 0)]
    
            this_length = gis_calc.edge_path_length(network=G, path=path)
    
            assert this_length == 20
    
            # *********************************************************************
            # *********************************************************************
    
            # no path
    
            path = [(0, 1, 0), (1, 3, 0)]
    
            this_length = gis_calc.edge_path_length(network=G, path=path)
    
            assert this_length == inf
    
            path = []
    
            this_length = gis_calc.edge_path_length(network=G, path=path)
    
            assert this_length == inf
    
            # *********************************************************************
            # *********************************************************************
    
        # *************************************************************************
        # *************************************************************************
    
        def example_count_occurrences(self, G: nx.MultiDiGraph):
            # get a gdf from the graph via osmnx
    
            gdf = ox.utils_graph.graph_to_gdfs(G, nodes=False)
    
            # define the column (street)
    
            column = "name"
    
            # count the ocurrences
    
            count_dict = gis_calc.count_ocurrences(gdf, column=column)
    
            # expected result
    
            true_count_dict = {
                nan: 249,
                "Havremarken": 8,
                "Kornmarken": 34,
                "Kløvervej": 38,
                "Kærvej": 52,
                "Plougslundvej": 24,
                "Egevænget": 16,
                "Fyrrevænget": 6,
                "Kærhusvej": 12,
                "Kløvermarken": 52,
                "Grønningen": 38,
                "Bygager": 10,
                "Fælleden": 6,
                "Flintemarken": 52,
                "Stendyssen": 8,
                "Markskellet": 54,
                "Engdraget": 20,
                "Vestervang": 36,
                "Tingstedet": 87,
                "Tuen": 10,
                "Lillevang": 96,
                "Grenevej": 24,
                "Hedegårdsvej": 16,
                "Gravhøjen": 28,
                "Lysningen": 37,
                "Ved Søen": 20,
                "Bopladsen": 10,
                "Koldingvej": 14,
                "Bakkelien": 38,
            }
    
            # test
    
            for key, value in count_dict.items():
                assert value == true_count_dict[key]
    
            # *********************************************************************
            # *********************************************************************
    
            # count only a few column entries
    
            count_dict = gis_calc.count_ocurrences(
                gdf, column=column, column_entries=["Kløvermarken", "Grenevej", "Bopladsen"]
            )
    
            # test
    
            for key, value in count_dict.items():
                assert value == true_count_dict[key]
    
            # *********************************************************************
            # *********************************************************************
    
            # count the nans
    
            count_dict = gis_calc.count_ocurrences(gdf, column=column, column_entries=[nan])
    
            # test
    
            for key, value in count_dict.items():
                assert value == true_count_dict[key]
    
            # *********************************************************************
            # *********************************************************************
    
        # *************************************************************************
        # *************************************************************************
    
        def test_compute_great_circle_distance_linestring(self):
            # Source: https://keisan.casio.com/exec/system/1224587128
    
            length_tolerance = 1  # meters
    
            # dict_A = {'x': 12.2430512, 'y': 55.6822924}
    
            # dict_B = {'x': 12.242863, 'y': 55.682155}
    
            # dict_C = {'x': 12.252233, 'y': 55.68085}
    
            list_1_points = [(12, 56)]
    
            # 1º longitude at a latitude of 56º
    
            list_2_points_a = [(12, 56), (13, 56)]
    
            # 1º latituide at a longitude of 13º
    
            list_2_points_b = [(13, 56), (13, 57)]
    
            # 1º longitude at a latitude of 56º + 1º latituide at a longitude of 13º
    
            list_3_points = [(12, 56), (13, 56), (13, 57)]
    
            # radius = 6363.478 km at sea level
            # edge = radius*(pi/180)*angle in degrees
    
            true_length_2_points_a = 62.178959 * 1e3  # 62.178959 km with r=6371.009 km
    
            true_length_2_points_b = 111.195084 * 1e3  # 111.195084 km with r=6371.009 km
    
            true_length_3_points = true_length_2_points_a + true_length_2_points_b
    
            # make sure the function fails with a single point (sanity check)
    
    Pedro L. Magalhães's avatar
    Pedro L. Magalhães committed
            try:
                line = LineString(list_1_points)
            except Exception:
    
    Pedro L. Magalhães's avatar
    Pedro L. Magalhães committed
    
            # make sure it works with 2 points
    
            line = LineString(list_2_points_a)
    
            _length = gis_calc.great_circle_distance_along_path(line)
    
            assert isclose(_length, true_length_2_points_a, abs_tol=length_tolerance)
    
            line = LineString(list_2_points_b)
    
            _length = gis_calc.great_circle_distance_along_path(line)
    
            assert isclose(_length, true_length_2_points_b, abs_tol=length_tolerance)
    
            # make sure it works with 3 points
    
            line = LineString(list_3_points)
    
            _length = gis_calc.great_circle_distance_along_path(line)
    
            assert isclose(_length, true_length_3_points, abs_tol=length_tolerance)
    
        # *************************************************************************
        # *************************************************************************
    
    
    # *****************************************************************************
    # *****************************************************************************