Rilevamento di cluster geografici

Question

Ho un data.frame R. contenente longitudine, latitudine che si estende su tutta la mappa USA. Quando X numero di voci sono tutte all'interno di una piccola area geografica di alcuni gradi di longitudine & di alcuni gradi di latitudine, voglio essere in grado di rilevare questo e quindi avere il mio programma quindi restituire le coordinate per il riquadro di delimitazione geografica. Esiste un pacchetto Python o R CRAN che lo fa già? In caso contrario, come farei per accertare queste informazioni?

Answer 1

Sono riuscito a combinare la risposta di Joran con il commento di Dan H. Questo è un esempio di uscita:

Il codice python emette funzioni per R: map() e rect(). Questo esempio mappa Stati Uniti d'America è stato creato con:

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

e poi si può applicare il rect() 's di conseguenza dal con nell'interprete R GUI (vedi sotto).

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() 

Ecco un esempio di file 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 il mio set di dati, l'uscita del mio script python, indicata sulla mappa USA. Ho cambiato i colori per chiarezza.

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")) 

---

aggiunta su 2013/05/01 per Yacob

---

Queste 2 linee ti danno la su tutto obiettivo ...

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

Se si vuole restringere in su una porzione di una mappa, è possibile utilizzare ylim e 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")) 

Si vuole utilizzare la mappa 'mondo' ...

map("world", plot=T) 

E 'passato molto tempo da quando ho usato questo codice Python che ho postato qui di seguito, quindi farò del mio meglio per aiutarti.

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. 

Ecco un esempio completo. Il file TSV si trova su pastebin.com. Ho anche incluso un'immagine generata da R che contiene l'output di tutti i comandi 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 2

forse qualcosa di simile

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) 

Answer 3

Alcune idee:

• Ad-hoc & approssimativa: il "2-D istogramma". Crea bins "rettangolari" arbitrari, della larghezza dei gradi a tua scelta, assegna a ogni bin un ID. Posizionare un punto in un raccoglitore significa "associare il punto all'ID del raccoglitore". Ad ogni aggiunta a un contenitore, chiedere al contenitore quanti punti ha. Lato negativo: non "vede" correttamente un gruppo di punti che si trovano a cavallo di un limite del bin; e: i contenitori di "larghezza longitudinale costante" sono in realtà (spazialmente) più piccoli mentre ci si sposta verso nord.
• Utilizzare la libreria "Shapely" per Python. Seguite il suo esempio standard per "buffering points" e fate un unione a cascata dei buffer. Cerca globi su una certa area o che "contengono" un certo numero di punti originali. Nota che Shapely non è intrinsecamente "geo-savy", quindi dovrai aggiungere correzioni se ne hai bisogno.
• Utilizzare un DB vero con elaborazione spaziale. MySQL, Oracle, Postgres (con PostGIS), MSSQL tutti (credo) hanno i tipi di dati "Geometry" e "Geography" e puoi eseguire query spaziali su di essi (dai tuoi script Python).

Ognuno di questi ha costi diversi in dollari e tempo (nella curva di apprendimento) ... e diversi gradi di accuratezza geospaziale. Devi scegliere quello che si adatta al tuo budget e/o requisiti.

Answer 4

Sto facendo questo su base regolare in primo luogo la creazione di una matrice di distanza e poi eseguire il clustering su di esso. Ecco il mio codice.

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) 

Non sono sicuro che risolva completamente il problema. Potresti voler provare con diversi k, e forse anche fare una seconda serie di cluster di alcuni dei primi cluster nel caso fossero troppo grandi, come se avessi un punto nel Minnesota e un migliaio in California. Quando si ha il gruppo points.to.group $, è possibile ottenere i riquadri di delimitazione individuando il massimo e il minimo di latitudine per gruppo.

Se vuoi che X abbia 20 anni e tu abbia 18 punti a New York e 22 a Dallas, devi decidere se vuoi una scatola piccola e una veramente grande (20 punti ciascuna), se è meglio avere avere la scatola di Dallas include 22 punti, o se si desidera dividere i 22 punti a Dallas in due gruppi. Il raggruppamento basato sulla distanza può essere buono in alcuni di questi casi. Ma ovviamente dipende dal motivo per cui vuoi raggruppare i punti.

/Chris

Answer 5

se si utilizza formosa, si potrebbe estendere la mia cluster_points function per tornare alla casella di delimitazione del cluster tramite la proprietà .bounds della geometria formosa, per esempio come questo:

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