I am working with street numbering and toponymy and I have an issue in hands:
Given some data provided by an Institution, i.e.:
- points as door numbering and
- 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:
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:
- Dark purple: class 3, avenues, roundabouts, plazas, aka main ways;
- Red: class 4, mostly streets with inferior classification to the later;
- 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.
PS: as you can see, some streets have the door numbering order against the arrows direction.
