de holle romp

het creëren van een clusterrand met behulp van een K-benadering van Dichtstbijzijnde buren

een paar maanden geleden schreef ik hier op Medium een artikel over het in kaart brengen van de Britse verkeersongevallen hot spots. Ik was vooral bezorgd over het illustreren van het gebruik van het dbscan clustering algoritme op geografische data. In het artikel gebruikte ik geografische informatie gepubliceerd door de Britse regering over gemelde verkeersongevallen. Mijn doel was om een dichtheidsgebaseerd clustering proces uit te voeren om de gebieden te vinden waar verkeersongevallen het vaakst worden gemeld. Het eindresultaat was de creatie van een set van geo-hekken vertegenwoordigen deze ongevallen hot spots.

door alle punten in een gegeven cluster te verzamelen, kunt u een idee krijgen over hoe het cluster eruit zou zien op een kaart, maar u zult een belangrijk stukje informatie missen: de buitenste vorm van het cluster. In dit geval hebben we het over een gesloten veelhoek die op een kaart kan worden weergegeven als een geo-fence. Elk punt binnen de geo-fence kan worden verondersteld tot de cluster te behoren, wat deze vorm een interessant stukje informatie maakt: je kunt het gebruiken als een discriminator functie. Alle nieuw bemonsterde punten die binnen de veelhoek vallen, kunnen worden verondersteld tot de overeenkomstige cluster te behoren. Zoals ik al aangaf in het artikel, je zou dergelijke polygonen kunnen gebruiken om je rijrisico te laten gelden, door ze te gebruiken om je eigen bemonsterde GPS-locatie te classificeren.

de vraag is nu hoe je een betekenisvolle veelhoek maakt uit een wolk van punten die een specifiek cluster vormen. Mijn aanpak in het eerste artikel was enigszins naïef en weerspiegelde een oplossing die ik al in productiecode had gebruikt. Deze oplossing omvatte het plaatsen van een cirkel gecentreerd op elk van de punten van de cluster en vervolgens het samenvoegen van alle cirkels samen om een wolk-vormige veelhoek te vormen. Het resultaat is niet erg mooi, noch realistisch. Door cirkels te gebruiken als basisvorm voor het bouwen van de uiteindelijke veelhoek, zullen deze meer punten hebben dan een meer gestroomlijnde vorm, waardoor de opslagkosten toenemen en de inclusiedetectiealgoritmen trager worden uitgevoerd.

aan de andere kant heeft deze benadering het voordeel dat ze rekenkundig eenvoudig is (tenminste vanuit het perspectief van de ontwikkelaar), omdat ze gebruik maakt van de functie cascaded_union om alle cirkels samen te voegen. Een ander voordeel is dat de vorm van de veelhoek impliciet wordt gedefinieerd met behulp van alle punten in de cluster.

voor meer geavanceerde benaderingen moeten we op de een of andere manier de randpunten van de cluster identificeren, die de puntwolkvorm lijken te definiëren. Interessant, met sommige dbscan implementaties kunt u daadwerkelijk herstellen van de grenspunten als een bijproduct van het clustering proces. Helaas is deze informatie (blijkbaar) niet beschikbaar op de implementatie van SciKit Learn, dus we moeten het doen.

de eerste benadering die in me opkwam was het berekenen van de convexe romp van de verzameling punten. Dit is een goed begrepen algoritme, maar lijdt aan het probleem van het niet omgaan met concave vormen, zoals deze:

deze vorm geeft de essentie van de onderliggende punten niet correct weer. Als het als discriminator zou worden gebruikt, zouden sommige punten ten onrechte worden geclassificeerd als zich binnen het cluster wanneer ze dat niet zijn. We hebben een andere aanpak nodig.

het Concave Rompalternatief

gelukkig zijn er alternatieven voor deze stand van zaken: we kunnen een concave romp berekenen. Hier is hoe de holle romp eruit ziet wanneer toegepast op dezelfde set van punten als in de vorige afbeelding:

of misschien deze:

zoals u kunt zien, en in tegenstelling tot de convexe romp, is er geen enkele definitie van wat de holle romp van een verzameling punten is. Met het algoritme dat ik hier presenteer, wordt de keuze van hoe concaaf je wilt dat je rompen worden gemaakt door middel van een enkele parameter: k — het aantal dichtstbijzijnde buren beschouwd tijdens de rompberekening. Eens kijken hoe dit werkt.

het algoritme

het algoritme dat ik hier presenteer is meer dan tien jaar geleden beschreven door Adriano Moreira en Maribel Yasmina Santos van de Universiteit van Minho, Portugal . Uit het abstract:

dit document beschrijft een algoritme om de envelop van een verzameling punten in een vlak te berekenen, die convex op niet-convexe rompen genereert die het gebied vertegenwoordigen dat door de gegeven punten wordt bezet. Het voorgestelde algoritme is gebaseerd op een k-dichtstbijzijnde buren aanpak, waarbij de waarde van k, de enige algoritme parameter, wordt gebruikt om de “gladheid” van de uiteindelijke oplossing te controleren.

omdat ik dit algoritme zal toepassen op geografische informatie, moesten enkele wijzigingen worden aangebracht, namelijk bij het berekenen van hoeken en afstanden . Maar deze veranderen op geen enkele manier de essentie van het algoritme, dat in grote lijnen kan worden beschreven door de volgende stappen:

1. zoek het punt met de laagste y (breedtegraad) coördinaat en maak het de huidige.
2. Zoek de k-punten die het dichtst bij het huidige punt liggen.
3. kies uit de k-dichtstbijzijnde punten die overeenkomt met de grootste bocht naar rechts van de vorige hoek. Hier zullen we het concept van lager gebruiken en beginnen met een hoek van 270 graden (pal West).
4. Controleer of door het toevoegen van het nieuwe punt aan de groeiende lijnstring, het zichzelf niet snijdt. Als dit het geval is, selecteer dan een ander punt uit de k-dichtstbijzijnde of herstart met een grotere waarde van k.
5. Maak het nieuwe punt het huidige punt en verwijder het uit de lijst.
6. Na k iteraties het eerste punt weer aan de lijst toevoegen.
7. lus naar nummer 2.

het algoritme lijkt vrij eenvoudig, maar er zijn een aantal details die aandacht moeten krijgen, vooral omdat we te maken hebben met geografische coördinaten. Afstanden en hoeken worden op een andere manier gemeten.

de Code

hier publiceer ik een aangepaste versie van de code van het vorige artikel. U vindt nog steeds dezelfde clustering code en dezelfde cloud-vormige cluster generator. De bijgewerkte versie bevat nu een pakket met de naam geomath.hullswaar u de klasse ConcaveHull kunt vinden. Ga als volgt te werk om uw holle rompen te maken:

in de bovenstaande code is points een array van dimensies (N, 2), waarin de rijen de waargenomen punten bevatten en de kolommen de geografische coördinaten (lengtegraad, breedtegraad) bevatten. De resulterende array heeft precies dezelfde structuur, maar bevat alleen de punten die in de veelhoekige vorm van de cluster horen. Een soort filter.

omdat we arrays zullen behandelen is het alleen maar natuurlijk om NumPy in de strijd te brengen. Alle berekeningen werden naar behoren vectorized, waar mogelijk, en inspanningen werden gedaan om de prestaties te verbeteren bij het toevoegen en verwijderen van items uit arrays (spoiler: ze worden niet verplaatst helemaal). Een van de ontbrekende verbeteringen is code parallellisatie. Maar dat kan wachten.

ik organiseerde de code rond het algoritme zoals weergegeven in het papier, hoewel enkele optimalisaties werden gemaakt tijdens de vertaling. Het algoritme is opgebouwd rond een aantal subroutines die duidelijk geïdentificeerd worden door het papier, dus laten we die nu uit de weg ruimen. Voor uw leescomfort gebruik ik dezelfde namen als in de krant.

CleanList – het opschonen van de lijst met punten wordt uitgevoerd in de klasse constructor:

zoals u kunt zien, de lijst met punten is geïmplementeerd als een NumPy array om prestatie-redenen. Het schoonmaken van de lijst wordt uitgevoerd op lijn 10, waar alleen de unieke punten worden bewaard. De data set array is georganiseerd met waarnemingen in rijen en geografische coördinaten in de twee kolommen. Merk op dat ik ook een Booleaanse array aan het maken ben op Regel 13 die zal worden gebruikt om te indexeren in de hoofddataset array, waardoor het karwei van het verwijderen van items wordt vergemakkelijkt en, af en toe, ze weer wordt toegevoegd. Ik heb deze techniek genaamd een “masker” op de NumPy documentatie gezien, en het is zeer krachtig. Wat de priemgetallen betreft, Ik zal ze later bespreken.

FindMinYPoint – dit vereist een kleine functie:

deze functie wordt aangeroepen met de dataset array als argument en retourneert de index van het punt met de laagste breedtegraad. Merk op dat rijen gecodeerd zijn met Lengtegraad in de eerste kolom en breedtegraad in de tweede.

RemovePoint
AddPoint-dit zijn geen-brainers, vanwege het gebruik van de indices array. Deze array wordt gebruikt om de actieve indices op te slaan in de hoofddataset array, dus het verwijderen van items uit de dataset is in een handomdraai.

hoewel het algoritme zoals beschreven in de paper vraagt om het toevoegen van een punt aan de array waaruit de romp bestaat, wordt dit in feite geïmplementeerd als:

Later zal de variabele test_hull worden toegewezen aan hull wanneer de regelstring wordt geacht niet te snijden. Maar ik loop voor op het spel. Het verwijderen van een punt uit de dataset array is zo eenvoudig als:

self.indices = False

het toevoegen van het terug is gewoon een kwestie van het omdraaien van die array waarde op dezelfde index terug naar waar. Maar al dit gemak komt met de prijs van het hebben van onze tabbladen op de indexen te houden. Meer hierover later.

NearestPoints-dingen beginnen hier interessant te worden omdat we niet te maken hebben met vlakke coördinaten, dus uit met Pythagoras en in met Haversine:

merk op dat de tweede en derde parameters arrays zijn in het datasetformaat: Lengtegraad in de eerste kolom en breedtegraad in de tweede. Zoals je kunt zien, geeft deze functie een array van afstanden in meters terug tussen het punt in het tweede argument en de punten in het derde argument. Zodra we deze hebben, kunnen we de dichtstbijzijnde buren op de gemakkelijke manier pakken. Maar daar is een gespecialiseerde functie voor en het verdient enige uitleg:

de functie begint met het maken van een array met de basisindices. Dit zijn de indices van de punten die niet uit de reeks gegevens zijn verwijderd. Bijvoorbeeld, als op een tien punt cluster we begonnen met het verwijderen van het eerste punt, de basisindices array zou zijn . Vervolgens berekenen we de afstanden en sorteren we de resulterende array indices. De eerste k worden geëxtraheerd en vervolgens gebruikt als een masker om de basisindices op te halen. Het is een beetje gekronkeld, maar het werkt. Zoals je kunt zien, retourneert de functie geen array van coördinaten, maar een array van indices in de data set array.

SortByAngle-er zijn meer problemen hier omdat we niet berekenen eenvoudige hoeken, we berekenen lagers. Deze worden gemeten als nul graden naar het noorden, waarbij de hoeken met de klok mee toenemen. Hier is de kern van de code die de lagers berekent:

de functie retourneert een array van lagers gemeten van het punt waarvan de index in het eerste argument is, tot de punten waarvan de indices in het derde argument zijn. Sorteren is eenvoudig:

op dit punt, de kandidaten array bevat de indexen van de k-dichtstbijzijnde punten gesorteerd op aflopende volgorde van lager.

IntersectQ-in plaats van mijn eigen lijnkruisfuncties te rollen, wendde ik me tot Shapely voor hulp. In feite, tijdens het bouwen van de veelhoek, behandelen we in wezen een lijnstring, waarbij we segmenten toevoegen die niet snijden met de vorige. Testen voor dit is eenvoudig: we pakken de romp array in aanbouw, converteren het naar een vormige lijn string object, en testen of het eenvoudig is (niet zelf-kruisend) of niet.

In een notendop wordt een vormige lijnstring complex als het zichzelf kruist, zodat het is_simple predicaat onwaar wordt. Gemakkelijk.

Pointpolygon — dit bleek de lastigste te implementeren. Laat me uitleggen door te kijken naar de code die de laatste romp polygoon validatie uitvoert (controleert of alle punten van het cluster zijn opgenomen in de polygoon):

Shapely ‘ s functies om te testen op kruising en inclusie zou genoeg moeten zijn om te controleren of de uiteindelijke romp polygoon overlapt alle punten van de cluster, maar ze waren niet. Waarom? Gevormd is coördinaat-agnostisch, dus het zal Geografische coördinaten uitgedrukt in breedtegraden en lengtegraden precies op dezelfde manier behandelen als coördinaten op een Cartesiaans vlak. Maar de wereld gedraagt zich anders als je op een bol leeft, en hoeken (of lagers) zijn niet constant langs een geodeet. Het voorbeeld over referentie van een geodetische lijn die Bagdad met Osaka verbindt, illustreert dit perfect. Het gebeurt zo dat onder sommige omstandigheden, het algoritme kan een punt op basis van het lager criterium bevatten, maar later, met behulp van Shapely ‘ s vlakke algoritmen, worden geacht ooit zo iets buiten de veelhoek. Dat is wat de kleine afstand correctie Daar doet.

het kostte me een tijdje om dit uit te zoeken. Mijn debughulp was QGIS, een geweldig stuk vrije software. Bij elke stap van de verdachte berekeningen zou ik de uitvoer van de gegevens in wkt-formaat in een CSV-bestand te lezen in als een laag. Een echte redder in nood!

ten slotte, als de veelhoek niet alle punten van het cluster dekt, is de enige optie om K te verhogen en het opnieuw te proberen. Hier voegde ik een beetje van mijn eigen intuïtie.

priemgetal k

het artikel suggereert dat de waarde van k met één wordt verhoogd en dat het algoritme opnieuw vanaf nul wordt uitgevoerd. Mijn vroege tests met deze optie waren niet erg bevredigend: de looptijd zou traag zijn op problematische clusters. Dit was te wijten aan de langzame toename van k, dus besloot ik om een ander verhogingsschema te gebruiken: een tabel met priemgetallen. Het algoritme begint al met k = 3, dus het was een gemakkelijke uitbreiding om het te laten evolueren op een lijst van priemgetallen. Dit is wat je ziet gebeuren in de recursieve aanroep:

ik heb iets met priemgetallen, Weet je …

Blow Up

de holle polygonen van de romp die door dit algoritme worden gegenereerd, moeten nog wat verder worden verwerkt, omdat ze alleen punten in de romp zullen onderscheiden, maar niet dichtbij. De oplossing is om wat padding toe te voegen aan deze skinny clusters. Hier gebruik ik exact dezelfde techniek als voorheen, en hier is hoe het eruit ziet:

hier heb ik Shapely ‘ s bufferfunctie gebruikt om de truc te doen.

de functie accepteert een gevormde veelhoek en geeft een opgeblazen versie van zichzelf terug. De tweede parameter is de straal in meters van de toegevoegde opvulling.

het uitvoeren van de Code

begin met het trekken van de code van de GitHub repository naar je lokale machine. Het bestand dat u wilt uitvoeren is ShowHotSpots.py in de hoofdmap. Bij de eerste uitvoering, de code zal lezen in het Verenigd Koninkrijk verkeersongevallengegevens van 2013 tot 2016 en cluster het. De resultaten worden vervolgens gecached als een CSV-bestand voor de volgende runs.

u krijgt dan twee kaarten te zien: de eerste wordt gegenereerd met behulp van de cloudvormige clusters, terwijl de tweede het hier besproken algoritme voor concave clustering gebruikt. Terwijl de polygoongeneratiecode wordt uitgevoerd, kunt u zien dat er een paar fouten worden gemeld. Om te begrijpen waarom het algoritme er niet in slaagt een holle romp aan te maken, schrijft de code de clusters naar CSV-bestanden naar de data/out/failed/ map. Zoals gewoonlijk kunt u QGIS gebruiken om deze bestanden als lagen te importeren.

in wezen faalt dit algoritme wanneer het niet genoeg punten vindt om de vorm” rond te gaan ” zonder zelf-kruising. Dit betekent dat u klaar moet zijn om deze clusters weg te gooien of om er een andere behandeling op toe te passen (convexe romp of samengevoegde luchtbellen).

Concave

het is een wrap. In dit artikel presenteerde ik een methode van post-processing dbscan gegenereerde geografische clusters in concave vormen. Deze methode kan een beter passende buiten veelhoek voor de clusters in vergelijking met andere alternatieven verstrekken.

Bedankt voor het lezen en veel plezier met het sleutelen aan de code!

Kryszkiewicz M., Lasek P. (2010) TI-DBSCAN: Clustering with DBSCAN by Means of the Triangle Inequality. In: Szczuka M., Kryszkiewicz M., Ramanna S., Jensen R., Hu Q. (eds) Rough Sets and Current Trends in Computing. RSCTC 2010. Lecture Notes in Computer Science, vol 6086. Springer, Berlin, Heidelberg

sikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011

Moreira, A. and Santos, M. Y., 2007, Concave Hull: A K-nearest neighbors approach for the computation of the region occupied by a set of points

Calculate distance, bearing and more between Latitude / Longitude points

GitHub repository