0

I am working with street numbering and toponymy and I have an issue in hands:

Given some data provided by an Institution, i.e.:

  1. points as door numbering and
  2. lines as the streets

Those streets were digitized randomly and not considering the door numbering. In this case, I came up with a road hierarchy which I classified from 1 to 5, being 1 the top important order, like highways, and 5 the less important ones, like neighbour access roads.

In the next picture I can show a sample:

enter image description here

This sample was retrieved from my QGIS project, showing in darker colours the highest ranked ways, and in lighter colours the lowest ranked ones. The colours represent:

  1. Dark purple: class 3, avenues, roundabouts, plazas, aka main ways;
  2. Red: class 4, mostly streets with inferior classification to the later;
  3. Orange: class 5, roads without expression, accesses to houses and without assessment of door numbering.

As you can also see, the arrows show the digitizing direction of the lines, as were provided. Where there's no arrows it means that those are ways which where the door numbering can't be created automatically, like squares, roundabouts, plazas, etc., at least for this project. Also, the lines that receive those arrows, are the only ones acceptable to have door numbering automatically defined. The ones that don't have arrows, somehow will have to be treated manually.

Anyways, the main idea here is to develop a script which checks for intersection between different features, and if a Class 4 road 'touches' a Class 3 road, then the digitizing direction should start from that intersection, i.e., from higher to lower classes. This classes have been created in a column in the lines shapefile, by hand, and are the main principle to define the digitizing direction, which will affect the points of house numbering.

I have already drawn a script, which returns no error, but doesn't do a thing to the digitizing direction. There are also some rules regarding the digitizing direction, like if the mean bearing is East-West (or vice-versa), digitizing should be from Right to Left, or is North-South (or vice-versa), should be South-North. The reader should be able to spot this in the following code:

# Carregar a shapefile
shapefile_path = r"path_to_input"
gdf = gpd.read_file(shapefile_path)

# Order GeoDataFrame by descending class
gdf = gdf.sort_values(by='Classe', ascending=False)

# Start a dictionary to store higher ways geometries
superior_lines = {}

# Iterate over lines
for index, row in gdf.iterrows():
    # print("Analysing feature with COD_VIA:", row['COD_VIATRO'])
    geometry = row['geometry']
    classe = row['Classe']
    tipo_parte = row['TipoParte']  # Add TipoParte column to the shapefile

    if tipo_parte == 2:  # Check if TipoParte == 2, i.e., if it's a sqaure, roundabout...
        continue  # continue if TipoParte == 2
    
    if geometry.geom_type == 'LineString':
        coords = list(geometry.coords)
        total_numerador = 0
        total_denominador = 0
        
        # Find higher way
        superior_line = None
        for superior_classe in range(classe - 1, 0, -1):
            if superior_classe in superior_lines:
                superior_line = superior_lines[superior_classe]
                break
        
        # Adjust line bearing
        adjusted_line = adjust_line_orientation(geometry, superior_line, False)
        
        # Update higher way geometry (for later use)
        superior_lines[classe] = adjusted_line
        
        # arctan calculation for the bearing
        for i in range(len(coords) - 1):
            Mpartida, Ppartida = coords[i]
            Mseguinte, Pseguinte = coords[i + 1]
            numerador = Mseguinte - Mpartida
            denominador = Pseguinte - Ppartida
            # print("Numerador:", numerador, "Denominador:", denominador)
            total_numerador += numerador
            total_denominador += denominador
        
        rumo_radianos = rumoArctg(total_numerador, 0, total_denominador, 0)
        rumo_graus = m.degrees(rumo_radianos)
        rumo_arredondado = round(rumo_graus, 2)
        # print("Rumo médio da feature (graus) - COD_VIA:", row['COD_VIA'], rumo_arredondado)
        
        # store bearing to the attribute table
        gdf.loc[index, 'RUMO'] = rumo_arredondado
        
    elif geometry.geom_type == 'MultiLineString':
        total_numerador = 0
        total_denominador = 0
        total_parts = 0
        
        for line_string in geometry.geoms:
            coords = list(line_string.coords)
            # Encontrar a via superior
            superior_line = None
            for superior_classe in range(classe - 1, 0, -1):
                if superior_classe in superior_lines:
                    superior_line = superior_lines[superior_classe]
                    break
            
            # adjust bearing
            adjusted_line = adjust_line_orientation(line_string, superior_line, False)
            
            # Update higher way geometry (for later use)
            superior_lines[classe] = adjusted_line
            
            for i in range(len(coords) - 1):
                Mpartida, Ppartida = coords[i]
                Mseguinte, Pseguinte = coords[i + 1]
                numerador = Mseguinte - Mpartida
                denominador = Pseguinte - Ppartida
                # print("Numerador:", numerador, "Denominador:", denominador) 
                total_numerador += numerador
                total_denominador += denominador
            total_parts += 1
        
        if total_parts > 0:
            rumo_radianos = rumoArctg(total_numerador, 0, total_denominador, 0)
            rumo_graus = m.degrees(rumo_radianos)
            rumo_medio = rumo_graus / total_parts
            rumo_arredondado = round(rumo_medio, 2)
            # print("Rumo médio da feature (graus) - COD_VIA:", row['COD_VIA'], rumo_arredondado)
            
            # Store bearing result
            gdf.loc[index, 'RUMO'] = rumo_arredondado
        
    else:
        print("Geometry type not supported:", geometry.geom_type)

# Save GeoDataFrame back to the shapefile with the new "RUMO" (bearing)
output_tratadasComRumos = r"path_to_output"
gdf.to_file(output_processedWithBearings)

What am I missing here?

This process can be done completely in a manual fashion, using the "reverse line" tool in QGIS. But since I am working with a big set of data, this whole script is meant to do it automatically, using only a manual intervention to correct the topology, add "Classe" and "TipoParte" to the attributes table.

enter image description here

PS: as you can see, some streets have the door numbering order against the arrows direction.

3
  • @BERA the this is the existing door numbering is in some cases against the street orientation. The intention is to reverse the street orientation when it comes from main and important places, like if some streets erradiate from a plaza, a square, etc. And also a street is from an inferior category like an avenue, so all streets that intersect from an avenue, as example, should also have their digitizing orientation starting on that/those avenue(s). If the door numbering is incorrect, following these rules, it should be corrected afterwards. Thank you
    – PatsonJS
    Commented Oct 2, 2023 at 11:48
  • What do you mean by "digitizing direction"? Do you just mean the direction of the nodes of the LineString/MultiLineString? Reversing lines in PyQGIS is covered in: gis.stackexchange.com/questions/313998/… Commented Oct 3, 2023 at 1:05
  • Can you share some sample data?
    – Bera
    Commented Oct 4, 2023 at 10:28

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.