Обнаружение географических кластеров

У меня есть R data.frame, содержащий долготу, широту, которая охватывает всю карту США. Когда X количество записей все в пределах небольшого географического региона, скажем, несколько градусов долготы & amp; в нескольких градусах широты, я хочу иметь возможность обнаружить это, а затем заставить мою программу вернуть координаты географической ограничительной рамки. Есть ли пакет Python или CRAN, который уже делает это? Если нет, то как мне узнать эту информацию?

 luke14free14 апр. 2012 г., 23:23
Вы пробовали кластеризацию k-средних с расстоянием Haversine?
 Ben Bolker11 апр. 2012 г., 17:30
Для R см.cran.r-project.org/web/views/Spatial.html (поиск "кластера")
 Josh O'Brien11 апр. 2012 г., 19:00
Вот этоinformative thread from R-sig-geo, Это начинается с обсужденияSaTScanЭто бесплатное (но не с открытым исходным кодом) программное обеспечение для решения вопросов, подобных вашему, и исследует наличие аналогичного программного обеспечения в декабре 2010 года.

Ответы на вопрос(5)

сначала создав матрицу расстояний, а затем запустив на ней кластеризацию. Вот мой код

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)

Я не уверен, что это полностью решит вашу проблему. Возможно, вы захотите протестировать с другим k, а также, возможно, выполнить второй прогон кластеризации некоторых из первых кластеров, если они слишком велики, например, если у вас есть одна точка в Миннесоте и тысяча в Калифорнии. Когда у вас есть группа points.to.group $, вы можете получить ограничивающие рамки, найдя max и min lat lon для каждой группы.

Если вы хотите, чтобы X было 20, и у вас есть 18 очков в Нью-Йорке и 22 в Далласе, вы должны решить, хотите ли вы одну маленькую и одну действительно большую коробку (20 очков каждая), лучше ли иметь коробку Далласа включите 22 балла, или если вы хотите разделить 22 балла в Далласе на две группы. Кластеризация на основе расстояния может быть хорошей в некоторых из этих случаев. Но это, конечно, зависит от того, почему вы хотите сгруппировать очки.

/Крис

функция cluster_points вернуть ограничивающую рамку кластера через свойство .bounds стройной геометрии, например, так:

clusterlist.append(cluster, (poly.buffer(-b)).bounds)
Решение Вопроса

вывода: cluster map

Код Python испускает функции для R: map () и rect (). Эта карта примера США была создана с:

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

и затем вы можете соответственно применять функции rect () из with в интерпретаторе R GUI (см. ниже).

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

Вот пример файла 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

С моим набором данных вывод моего скрипта на python показан на карте США. Я изменил цвета для ясности.

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

Дополнение 2013-05-01 для Yacob

Эти 2 строки дают вам общую цель ...

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

Если вы хотите сузить часть карты, вы можете использоватьylim а также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"))

Вам захочется использовать «мир» карта...

map("world", plot=T )

Прошло много времени с тех пор, как я использовал этот код Python, который я разместил ниже, поэтому я постараюсь помочь вам.

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.

Вот полный пример. Файл TSV находится на pastebin.com. Я также включил изображение, сгенерированное из R, которое содержит выходные данные всех команд 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()

pyclusters.py / R image result

 jftuga13 мая 2013 г., 19:29
@Yacob: да, это правильно.
 jftuga02 мая 2013 г., 17:18
@Yacob: Надеюсь, это поможет!
 07 мая 2013 г., 15:16
Если я вас понимаю, я должен запустить ваш скрипт на python как ./pyclusters.py для команды linux, а затем использовать вывод для построения в R?
 25 февр. 2015 г., 09:21
Какую библиотеку вы используете в R для построения карты?
 30 апр. 2013 г., 17:37
Привет, Jftuga, я хотел бы использовать твой скрипт на python для географической кластеризации точек на основе широты и долготы по всему миру. Подскажите, пожалуйста, как мне это сделать.

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)
 11 апр. 2012 г., 17:24
вот почему я оставил его в своей собственной функции ... Я не уверен, как рассчитать широту и долготу ...
 11 апр. 2012 г., 23:32
спасибо за это пригодится, если мне нужно возиться с широтой
 11 апр. 2012 г., 17:12
Как правило, использование лат и долг в качестве декартовых координат, равных расстоянию, является действительно плохой идеей (как вы делаете с вычислением «гипотенузы» выше).
 11 апр. 2012 г., 21:31
Если вам нужно расстояние между парами широта / долгота, возможно, это лучший ресурс в Интернете для алгоритмов и объяснений:movable-type.co.uk/scripts/latlong.html, Много разных формул, в зависимости от вашей точности нужен бюджет. Для расстояния 100 км или около того (около градуса или около того) используется «Равноугольное приближение». это достаточно хорошо для многих целей. Вы должны настроить долготу с помощью косинуса средней широты, например: R * sqrt ((lat1-lat2) ** 2 + ((lon1-lon2) * cos ((lat1 + lat2) / 2)) ** 2), где R - радиус Земли (в той же единице, в которой вы хотите получить результат).

Ad-hoc & approximate: The "2-D histogram". Create arbitrary "rectangular" bins, of the degree width of your choice, assign each bin an ID. Placing a point in a bin means "associate the point with the ID of the bin". Upon each add to a bin, ask the bin how many points it has. Downside: doesn't correctly "see" a cluster of points that stradle a bin boundary; and: bins of "constant longitudinal width" actually are (spatially) smaller as you move north. Use the "Shapely" library for Python. Follow it's stock example for "buffering points", and do a cascaded union of the buffers. Look for globs over a certain area, or that "contain" a certain number of original points. Note that Shapely is not intrinsically "geo-savy", so you'll have to add corrections if you need them. Use a true DB with spatial processing. MySQL, Oracle, Postgres (with PostGIS), MSSQL all (I think) have "Geometry" and "Geography" datatypes, and you can do spatial queries on them (from your Python scripts).

Каждый из них имеет разные затраты в долларах и времени (в кривой обучения) ... и разные степени геопространственной точности. Вы должны выбрать то, что соответствует вашему бюджету и / или требованиям.

Ваш ответ на вопрос