User:DressyPear4/WaterwayCulvert

Waterway culvert

O código tem como objetivo fazer divisão de qualquer linha com tag waterway=* e acrescentar tunnel=culvert + layer=-1 onde há interseção com highway. Nesse caso ele intera sobre o dataset, caso a quantidade de dados for muito grande pode demorar um pouco para finalizar.

Como funciona?

• Selecionar uma linha com tag waterway=*
• Executar o script
• Caso na interseção a linha com tag waterway=* já tiver tunnel=culvert + layer=-1, a divisão nesse local não será feita.
• Caso na interseção a linha com tag highway=* já tiver bridge=yes + layer=1, a divisão nesse local não será feita.

Demonstração

Imagem.gif, clique para visualizar.

Código

import math
from javax.swing import UIManager
from org.openstreetmap.josm.data.osm import Node, Way
from org.openstreetmap.josm.data.coor import LatLon
from org.openstreetmap.josm.gui import MainApplication, Notification
from org.openstreetmap.josm.data import UndoRedoHandler
from org.openstreetmap.josm.command import SequenceCommand, AddCommand, ChangeNodesCommand, DeleteCommand

def main():
    layer = MainApplication.getLayerManager().getActiveLayer()
    if layer is None or not hasattr(layer, "data"):
        Notification("Nenhuma camada de dados ativa encontrada")\
            .setIcon(UIManager.getIcon("OptionPane.errorIcon")).show()
        return

    ds = layer.data

    waterways = [w for w in ds.getSelectedWays() if w.get("waterway")]
    if not waterways:
        Notification("Nenhum waterway selecionado")\
            .setIcon(UIManager.getIcon("OptionPane.warningIcon")).show()
        return

    highways = [w for w in ds.getWays() if w.get("highway") and w.get("bridge") is None]

    def calculate_bearing(lat1, lon1, lat2, lon2):
        phi1 = math.radians(lat1)
        phi2 = math.radians(lat2)
        dlon = math.radians(lon2 - lon1)
        x = math.sin(dlon) * math.cos(phi2)
        y = math.cos(phi1)*math.sin(phi2) - math.sin(phi1)*math.cos(phi2)*math.cos(dlon)
        bearing = math.atan2(x, y)
        return (math.degrees(bearing) + 360) % 360

    def offset_point(lat, lon, distance, bearing):
        R = 6378137.0
        lat1 = math.radians(lat)
        lon1 = math.radians(lon)
        bearing_rad = math.radians(bearing)
        lat2 = math.asin(math.sin(lat1)*math.cos(distance/R) + math.cos(lat1)*math.sin(distance/R)*math.cos(bearing_rad))
        lon2 = lon1 + math.atan2(math.sin(bearing_rad)*math.sin(distance/R)*math.cos(lat1), math.cos(distance/R)-math.sin(lat1)*math.sin(lat2))
        return math.degrees(lat2), math.degrees(lon2)

    def seg_intersection(p1, p2, q1, q2):
        avg_lat = math.radians((p1[0]+p2[0]+q1[0]+q2[0])/4.0)
        def to_xy(lat, lon):
            x = math.radians(lon)*math.cos(avg_lat)
            y = math.radians(lat)
            return x, y
        p1x, p1y = to_xy(p1[0], p1[1])
        p2x, p2y = to_xy(p2[0], p2[1])
        q1x, q1y = to_xy(q1[0], q1[1])
        q2x, q2y = to_xy(q2[0], q2[1])
        denom = (p2x-p1x)*(q2y-q1y) - (p2y-p1y)*(q2x-q1x)
        if abs(denom) < 1e-15:
            return None
        t = ((q1x-p1x)*(q2y-q1y) - (q1y-p1y)*(q2x-q1x))/denom
        u = ((q1x-p1x)*(p2y-p1y) - (q1y-p1y)*(p2x-p1x))/denom
        if t < 0 or t > 1 or u < 0 or u > 1:
            return None
        ix = p1x + t*(p2x-p1x)
        iy = p1y + t*(p2y-p1y)
        lat_int = math.degrees(iy)
        lon_int = math.degrees(ix / math.cos(avg_lat))
        return (lat_int, lon_int), t

    commands = []
    total_culverts = 0

    for waterway in waterways:
        if waterway.get("tunnel") == "culvert":
            continue
        w_nodes = list(waterway.getNodes())

        new_nodes_info = []

        for i in range(len(w_nodes) - 1):
            p1 = (w_nodes[i].getCoor().lat(), w_nodes[i].getCoor().lon())
            p2 = (w_nodes[i+1].getCoor().lat(), w_nodes[i+1].getCoor().lon())

            for highway in highways:
                hw_nodes = list(highway.getNodes())
                for j in range(len(hw_nodes) - 1):
                    q1 = (hw_nodes[j].getCoor().lat(), hw_nodes[j].getCoor().lon())
                    q2 = (hw_nodes[j+1].getCoor().lat(), hw_nodes[j+1].getCoor().lon())

                    shared_nodes = set([w_nodes[i], w_nodes[i+1]]) & set([hw_nodes[j], hw_nodes[j+1]])
                    if shared_nodes:
                        continue

                    inter_data = seg_intersection(p1, p2, q1, q2)
                    if inter_data is None:
                        continue

                    inter, t = inter_data

                    bearing_to_p1 = calculate_bearing(inter[0], inter[1], p1[0], p1[1])
                    bearing_to_p2 = calculate_bearing(inter[0], inter[1], p2[0], p2[1])
                    lat_off1, lon_off1 = offset_point(inter[0], inter[1], 5, bearing_to_p1)
                    lat_off2, lon_off2 = offset_point(inter[0], inter[1], 5, bearing_to_p2)

                    node_before = Node(LatLon(lat_off1, lon_off1))
                    node_after = Node(LatLon(lat_off2, lon_off2))

                    new_nodes_info.append((i, t, node_before, node_after))

        if not new_nodes_info:
            continue

        # Ordenar inserções por índice do segmento + parâmetro t na linha
        new_nodes_info.sort(key=lambda x: (x[0], x[1]))

        # Inserir nós do fim para o início para não bagunçar os índices
        for i, t, node_before, node_after in reversed(new_nodes_info):
            commands.append(AddCommand(ds, node_before))
            commands.append(AddCommand(ds, node_after))
            w_nodes.insert(i+1, node_after)
            w_nodes.insert(i+1, node_before)

        offset_nodes = []
        for _, _, node_before, node_after in new_nodes_info:
            offset_nodes.append(node_before)
            offset_nodes.append(node_after)

        offset_indices = sorted([w_nodes.index(n) for n in offset_nodes])

        # Preparar segmentos sem sobreposição e com sequência correta
        segments = []
        prev = 0
        idx = 0
        while idx < len(offset_indices):
            start = prev
            mid = offset_indices[idx]
            idx += 1
            if idx < len(offset_indices):
                end = offset_indices[idx]
                idx += 1
            else:
                end = len(w_nodes) - 1

            segments.append( (start, mid) )
            segments.append( (mid, end) )
            prev = end

        if prev < len(w_nodes)-1:
            segments.append( (prev, len(w_nodes)-1) )

        # Encontrar maior segmento (maior número de nós)
        largest_segment = max(range(len(segments)), key=lambda i: segments[i][1]-segments[i][0])

        new_ways = []
        for i, (start, end) in enumerate(segments):
            if start > end:
                segment_nodes = list(reversed(w_nodes[end:start+1]))
            else:
                segment_nodes = w_nodes[start:end+1]

            if i == largest_segment:
                commands.append(ChangeNodesCommand(waterway, segment_nodes))
            else:
                way_seg = Way()
                way_seg.setNodes(segment_nodes)
                for k, v in waterway.getKeys().items():
                    way_seg.put(k, v)
                if i % 2 == 1:
                    way_seg.put("tunnel", "culvert")
                    way_seg.put("layer", "-1")
                new_ways.append(way_seg)

        for way in new_ways:
            commands.append(AddCommand(ds, way))

        total_culverts += len(new_nodes_info)

    if total_culverts > 0:
        UndoRedoHandler.getInstance().add(SequenceCommand("Dividir waterways com culvert", commands))
        Notification(u"Processo concluído: {} culverts criados".format(total_culverts))\
            .setIcon(UIManager.getIcon("OptionPane.informationIcon")).show()
    else:
        Notification(u"Nenhuma interseção com highway encontrada."
                     u"\nSe na interseção tiver tags com 'culvert' 'bridge',"
                     u"\nou ambas compartilharem um nó em comum,"
                     u"\na divisao não e feita.")\
            .setIcon(UIManager.getIcon("OptionPane.warningIcon")).show()

main()