sklearn - plot k means python

Detectando clusters geográficos

Question

Tengo un cuadro de datos R que contiene longitud, latitud que se extiende sobre todo el mapa de EE. UU. Cuando X número de entradas están todas dentro de una pequeña región geográfica de unos pocos grados de longitud y unos pocos grados de latitud, quiero poder detectar esto y luego hacer que mi programa devuelva las coordenadas para el cuadro delimitador geográfico. ¿Hay algún paquete Python o R CRAN que ya haga esto? Si no es así, ¿cómo voy a averiguar esta información?

Answer 1

Algunas ideas:

• Ad-hoc y aproximado: El "histograma 2-D". Cree contenedores arbitrarios "rectangulares", del ancho de grado de su elección, asigne una ID a cada contenedor. Colocar un punto en un contenedor significa "asociar el punto con el ID del contenedor". Después de cada agregado a un recipiente, pregúntele cuántos puntos tiene. Desventaja: no "ve" correctamente un conjunto de puntos que se extienden sobre un límite de bin; y: los contenedores de "ancho longitudinal constante" en realidad son (espacialmente) más pequeños a medida que te mueves hacia el norte.
• Utilice la biblioteca "Shapely" para Python. Siga su ejemplo de stock para "puntos de búfer" y realice una unión en cascada de los búferes. Busque globos sobre un área determinada, o que "contengan" un cierto número de puntos originales. Tenga en cuenta que Shapely no es intrínsecamente "geo-savy", por lo que tendrá que agregar correcciones si las necesita.
• Utilice una base de datos verdadera con procesamiento espacial. MySQL, Oracle, Postgres (con PostGIS), MSSQL (creo) tienen tipos de datos de "Geometría" y "Geografía", y puede hacer consultas espaciales en ellos (desde sus scripts de Python).

Cada uno de estos tiene diferentes costos en dólares y tiempo (en la curva de aprendizaje) ... y diferentes grados de precisión geoespacial. Usted tiene que elegir lo que se adapte a su presupuesto y / o requisitos.

Answer 2

Hago esto de forma regular creando primero una matriz de distancia y luego ejecutando el agrupamiento en clústeres. Aquí está mi código.

library(geosphere)
library(cluster)
clusteramounts <- 10
distance.matrix <- (distm(points.to.group[,c("lon","lat")]))
clustersx <- as.hclust(agnes(distance.matrix, diss = T))
points.to.group$group <- cutree(clustersx, k=clusteramounts)

No estoy seguro de si soluciona completamente su problema. Es posible que desee realizar pruebas con diferentes k, y también hacer una segunda ejecución de agrupación de algunos de los primeros clústeres en caso de que sean demasiado grandes, como si tuviera un punto en Minnesota y mil en California. Cuando tenga el grupo points.to.group $, puede obtener los cuadros delimitadores si encuentra el máximo y mínimo lat lon por grupo.

Si quieres que X tenga 20, y tienes 18 puntos en Nueva York y 22 en Dallas, debes decidir si quieres una pequeña y una realmente grande (20 puntos cada una), si es mejor tener la caja de Dallas incluye 22 puntos, o si quieres dividir los 22 puntos en Dallas en dos grupos. La agrupación basada en la distancia puede ser buena en algunos de estos casos. Pero, por supuesto, depende de por qué quieres agrupar los puntos.

/ Chris

Answer 3

Pude combinar la respuesta de Joran con el comentario de Dan H. Este es un ejemplo de salida:

El código de Python emite funciones para R: map () y rect (). Este mapa de ejemplo de EE. UU. Se creó con:

map(''state'', plot = TRUE, fill = FALSE, col = palette())

y luego puede aplicar los rect () en consecuencia desde el intérprete de la GUI R (ver más abajo).

import math
from collections import defaultdict

to_rad = math.pi / 180.0   # convert lat or lng to radians
fname = "site.tsv"        # file format: LAT/tLONG
threshhold_dist=50         # adjust to your needs
threshhold_locations=15    # minimum # of locations needed in a cluster

def dist(lat1,lng1,lat2,lng2):
    global to_rad
    earth_radius_km = 6371

    dLat = (lat2-lat1) * to_rad
    dLon = (lng2-lng1) * to_rad
    lat1_rad = lat1 * to_rad
    lat2_rad = lat2 * to_rad

    a = math.sin(dLat/2) * math.sin(dLat/2) + math.sin(dLon/2) * math.sin(dLon/2) * math.cos(lat1_rad) * math.cos(lat2_rad)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)); 
    dist = earth_radius_km * c
    return dist

def bounding_box(src, neighbors):
    neighbors.append(src)
    # nw = NorthWest se=SouthEast
    nw_lat = -360
    nw_lng = 360
    se_lat = 360
    se_lng = -360

    for (y,x) in neighbors:
        if y > nw_lat: nw_lat = y
        if x > se_lng: se_lng = x

        if y < se_lat: se_lat = y
        if x < nw_lng: nw_lng = x

    # add some padding
    pad = 0.5
    nw_lat += pad
    nw_lng -= pad
    se_lat -= pad
    se_lng += pad

    # sutiable for r''s map() function
    return (se_lat,nw_lat,nw_lng,se_lng)

def sitesDist(site1,site2): 
    #just a helper to shorted list comprehension below 
    return dist(site1[0],site1[1], site2[0], site2[1])

def load_site_data():
    global fname
    sites = defaultdict(tuple)

    data = open(fname,encoding="latin-1")
    data.readline() # skip header
    for line in data:
        line = line[:-1]
        slots = line.split("/t")
        lat = float(slots[0])
        lng = float(slots[1])
        lat_rad = lat * math.pi / 180.0
        lng_rad = lng * math.pi / 180.0
        sites[(lat,lng)] = (lat,lng) #(lat_rad,lng_rad)
    return sites

def main():
    sites_dict = {}
    sites = load_site_data()
    for site in sites: 
        #for each site put it in a dictionary with its value being an array of neighbors 
        sites_dict[site] = [x for x in sites if x != site and sitesDist(site,x) < threshhold_dist] 

    results = {}
    for site in sites: 
        j = len(sites_dict[site])
        if j >= threshhold_locations:
            coord = bounding_box( site, sites_dict[site] )
            results[coord] = coord

    for bbox in results:
        yx="ylim=c(%s,%s), xlim=c(%s,%s)" % (results[bbox]) #(se_lat,nw_lat,nw_lng,se_lng)
        print(''map("county", plot=T, fill=T, col=palette(), %s)'' % yx)
        rect=''rect(%s,%s, %s,%s, col=c("red"))'' % (results[bbox][2], results[bbox][0], results[bbox][3], results[bbox][2])
        print(rect)
        print("")

main()

Aquí hay un ejemplo de archivo TSV (site.tsv)

LAT     LONG
36.3312 -94.1334
36.6828 -121.791
37.2307 -121.96
37.3857 -122.026
37.3857 -122.026
37.3857 -122.026
37.3895 -97.644
37.3992 -122.139
37.3992 -122.139
37.402  -122.078
37.402  -122.078
37.402  -122.078
37.402  -122.078
37.402  -122.078
37.48   -122.144
37.48   -122.144
37.55   126.967

Con mi conjunto de datos, la salida de mi script de python, que se muestra en el mapa de EE.UU. Cambié los colores para mayor claridad.

rect(-74.989,39.7667, -73.0419,41.5209, col=c("red"))
rect(-123.005,36.8144, -121.392,38.3672, col=c("green"))
rect(-78.2422,38.2474, -76.3,39.9282, col=c("blue"))

Adición el 2013-05-01 para Yacob

Estas 2 líneas te dan la meta general ...

map("county", plot=T )
rect(-122.644,36.7307, -121.46,37.98, col=c("red"))

Si desea reducir una parte de un mapa, puede usar ylim y xlim

map("county", plot=T, ylim=c(36.7307,37.98), xlim=c(-122.644,-121.46))
# or for more coloring, but choose one or the other map("country") commands
map("county", plot=T, fill=T, col=palette(), ylim=c(36.7307,37.98), xlim=c(-122.644,-121.46))
rect(-122.644,36.7307, -121.46,37.98, col=c("red"))

Usted querrá usar el mapa ''mundo'' ...

map("world", plot=T )

Ha pasado mucho tiempo desde que usé este código de Python que he publicado a continuación, así que haré todo lo posible para ayudarte.

threshhold_dist is the size of the bounding box, ie: the geographical area
theshhold_location is the number of lat/lng points needed with in
    the bounding box in order for it to be considered a cluster.

Aquí hay un ejemplo completo. El archivo TSV se encuentra en pastebin.com. También he incluido una imagen generada desde R que contiene la salida de todos los comandos rect ().

# pyclusters.py
# May-02-2013
# -John Taylor

# latlng.tsv is located at http://pastebin.com/cyvEdx3V
# use the "RAW Paste Data" to preserve the tab characters

import math
from collections import defaultdict

# See also: http://www.geomidpoint.com/example.html
# See also: http://www.movable-type.co.uk/scripts/latlong.html

to_rad = math.pi / 180.0  # convert lat or lng to radians
fname = "latlng.tsv"      # file format: LAT/tLONG
threshhold_dist=20        # adjust to your needs
threshhold_locations=20   # minimum # of locations needed in a cluster
earth_radius_km = 6371

def coord2cart(lat,lng):
    x = math.cos(lat) * math.cos(lng)
    y = math.cos(lat) * math.sin(lng)
    z = math.sin(lat)
    return (x,y,z)

def cart2corrd(x,y,z):
    lon = math.atan2(y,x)
    hyp = math.sqrt(x*x + y*y)
    lat = math.atan2(z,hyp)
    return(lat,lng)

def dist(lat1,lng1,lat2,lng2):
    global to_rad, earth_radius_km

    dLat = (lat2-lat1) * to_rad
    dLon = (lng2-lng1) * to_rad
    lat1_rad = lat1 * to_rad
    lat2_rad = lat2 * to_rad

    a = math.sin(dLat/2) * math.sin(dLat/2) + math.sin(dLon/2) * math.sin(dLon/2) * math.cos(lat1_rad) * math.cos(lat2_rad)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)); 
    dist = earth_radius_km * c
    return dist

def bounding_box(src, neighbors):
    neighbors.append(src)
    # nw = NorthWest se=SouthEast
    nw_lat = -360
    nw_lng = 360
    se_lat = 360
    se_lng = -360

    for (y,x) in neighbors:
        if y > nw_lat: nw_lat = y
        if x > se_lng: se_lng = x

        if y < se_lat: se_lat = y
        if x < nw_lng: nw_lng = x

    # add some padding
    pad = 0.5
    nw_lat += pad
    nw_lng -= pad
    se_lat -= pad
    se_lng += pad

    #print("answer:")
    #print("nw lat,lng : %s %s" % (nw_lat,nw_lng))
    #print("se lat,lng : %s %s" % (se_lat,se_lng))

    # sutiable for r''s map() function
    return (se_lat,nw_lat,nw_lng,se_lng)

def sitesDist(site1,site2): 
    # just a helper to shorted list comprehensioin below 
    return dist(site1[0],site1[1], site2[0], site2[1])

def load_site_data():
    global fname
    sites = defaultdict(tuple)

    data = open(fname,encoding="latin-1")
    data.readline() # skip header
    for line in data:
        line = line[:-1]
        slots = line.split("/t")
        lat = float(slots[0])
        lng = float(slots[1])
        lat_rad = lat * math.pi / 180.0
        lng_rad = lng * math.pi / 180.0
        sites[(lat,lng)] = (lat,lng) #(lat_rad,lng_rad)
    return sites

def main():
    color_list = ( "red", "blue", "green", "yellow", "orange", "brown", "pink", "purple" )
    color_idx = 0
    sites_dict = {}
    sites = load_site_data()
    for site in sites: 
        #for each site put it in a dictionarry with its value being an array of neighbors 
        sites_dict[site] = [x for x in sites if x != site and sitesDist(site,x) < threshhold_dist] 

    print("")
    print(''map("state", plot=T)'') # or use: county instead of state
    print("")


    results = {}
    for site in sites: 
        j = len(sites_dict[site])
        if j >= threshhold_locations:
            coord = bounding_box( site, sites_dict[site] )
            results[coord] = coord

    for bbox in results:
        yx="ylim=c(%s,%s), xlim=c(%s,%s)" % (results[bbox]) #(se_lat,nw_lat,nw_lng,se_lng)

        # important!
        # if you want an individual map for each cluster, uncomment this line
        #print(''map("county", plot=T, fill=T, col=palette(), %s)'' % yx)
        if len(color_list) == color_idx:
            color_idx = 0
        rect=''rect(%s,%s, %s,%s, col=c("%s"))'' % (results[bbox][2], results[bbox][0], results[bbox][3], results[bbox][1], color_list[color_idx])
        color_idx += 1
        print(rect)
    print("")


main()

Answer 4

si usa bien, podría ampliar mi función de cluster_points para devolver el cuadro delimitador del clúster a través de la propiedad .bounds de la geometría bien formada, por ejemplo, como este:

clusterlist.append(cluster, (poly.buffer(-b)).bounds)

Answer 5

tal vez algo como

def dist(lat1,lon1,lat2,lon2):
    #just return normal x,y dist
    return sqrt((lat1-lat2)**2+(lon1-lon2)**2)

def sitesDist(site1,site2):
    #just a helper to shorted list comprehensioin below
    return dist(site1.lat,site1.lon,site2.lat,site2.lon)
sites_dict = {}
threshhold_dist=5 #example dist
for site in sites:
    #for each site put it in a dictionarry with its value being an array of neighbors
    sites_dict[site] = [x for x in sites if x != site and sitesDist(site,x) < threshhold_dist]
print "/n".join(sites_dict)