[NUUG kart] Forenkle kystlinjer o.a. med JOSM "Delete unnecessary nodes"

[Steinar Hamre]

On Mon, Jul 07, 2008 at 05:29:38PM +0200, Helge Hafting wrote:
> Kystlinjer inneholder ofte massevis av korte linjestykker som kommer 
> rett etter hverandre. Disse kan erstattes av en enkelt rett linje,
> noe som sparer massevis av punkter. Kartet er like bra, men
> mindre data å prosessere.

Dette er en av flere sideeffekter av den ikke-helt-optimale måten
disse kystlinjene har blitt laget og konvertert på.

> JOSM har en knapp "delete unnecessary nodes" som gjør akkurat dette.

Hmm... Hvor, og hvilken versjon? Hvilke begrensninger har den?

> Er det noen grunn til å ikke bruke denne funksjonen?


Jeg startet for en stund tilbake med noe kode for å gjøre noe med dette.
(tanken var å kjøre det før import av kystlinjene i Nordland,
men så kom noen meg i forkjøpet...)

Parametrene (derror < 0.017 and lerror < 0.000009)
er tillat feilmargin i vinkel (radianer) og avstand (grader)...
0.017 tilsvarer en grad sving og 0.000009 tilsvarer en meter.

Hvis det er interesse, så er det en lett sak  å få den til å generere
differ til OSM... (enklere enn å laste inn hele Norge i JOSM)


import math
from cElementTree import ElementTree

et = ElementTree()

tree = et.parse("test.osm")

allnodes = {}
delete = []

for node in tree.getiterator('node'):
    allnodes[node.get('id')] = (float(node.get('lat')), float(node.get('lon')))

for way in tree.getiterator('way'):
    tags={}
    for tag in way.getiterator("tag"):
        tags[tag.get('k')] = tag.get('v')
    if tags.get('created_by') != 'almien_coastlines': break
    if tags.get('natural') != 'coastline': break
    source = tags.get('source')
    if not source or len(source) < 3 or source[:3] != 'PGS': break
    waynodes=[ (nd.get('ref'),) + allnodes[nd.get('ref')]
               for nd in way.getiterator("nd") ]
    sid, slat, slon = waynodes[0]
    mid, mlat, mlon = waynodes[1]
    print way.get('id')
    for id,lat,lon in waynodes[2:]:
        mdir = math.atan2(mlat-slat, mlon-slon)
        dir = math.atan2(lat-slat, lon-slon)
        length = math.hypot(slat-lat, slon-lon)
        derror = abs(dir-mdir)
        lerror = length * abs(dir-mdir)
        print mid, derror*180/math.pi, lerror*1e7/90, length*1e7/90
        # error less than 1 meter and 1 degree in heading
        if (derror < 0.017 and lerror < 0.000009):
            delete.append(mid)
            print "delete"
            mid = id
        else:
            sid, slat, slon = mid, mlat, mlon
            mid, mlat, mlon = id, lat, lon
            
print len(allnodes),len(delete)