Phát hiện các cụm địa lý

Question

Tôi có một khung dữ liệu R. chứa kinh độ, vĩ độ trải rộng trên toàn bộ bản đồ Hoa Kỳ. Khi số lượng mục X là tất cả trong một vùng địa lý nhỏ có nghĩa là kinh độ một vài độ & một vài độ vĩ độ, tôi muốn có thể phát hiện điều này và sau đó có chương trình của tôi rồi trả về tọa độ cho hộp giới hạn địa lý. Có một gói Python hay R CRAN đã thực hiện điều này? Nếu không, làm thế nào tôi sẽ đi về xác định thông tin này?

Answer 1

Tôi đã có thể kết hợp câu trả lời của Joran cùng với nhận xét của Dan H. Đây là ví dụ về ouput:

Mã python phát ra các hàm cho R: map() và rect(). Bản đồ ví dụ ở Hoa Kỳ này đã được tạo với:

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

và sau đó bạn có thể áp dụng rect() cho phù hợp với trình thông dịch R GUI (xem bên dưới).

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

Dưới đây là một ví dụ tập tin 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 

Với bộ dữ liệu của tôi, đầu ra của kịch bản python của tôi, hiển thị trên bản đồ Hoa Kỳ. Tôi đã thay đổi màu sắc cho sự rõ ràng.

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

---

Addition trên 2013/05/01 cho Yacob

---

Những 2 dây chuyền cung cấp cho bạn trên tất cả các mục tiêu ...

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

Nếu bạn muốn thu hẹp ở trên một phần của bản đồ, bạn có thể sử dụng ylim và 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")) 

Bạn sẽ muốn sử dụng bản đồ 'thế giới' ...

map("world", plot=T) 

Đã lâu rồi kể từ khi tôi đã sử dụng mã python này tôi đã đăng bên dưới để tôi sẽ cố hết sức để giúp bạn.

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. 

Dưới đây là ví dụ hoàn chỉnh. Tệp TSV nằm trên pastebin.com. Tôi cũng đã bao gồm một hình ảnh được tạo ra từ R có chứa đầu ra của tất cả các lệnh 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

có lẽ cái gì đó như

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

Một vài ý tưởng:

• Ad-hoc & gần đúng: Các "2-D histogram". Tạo các thùng "hình chữ nhật" tùy ý, theo chiều rộng độ của sự lựa chọn của bạn, gán mỗi bin một ID. Đặt một điểm vào thùng có nghĩa là "liên kết điểm với ID của thùng". Sau mỗi lần thêm vào thùng, hãy hỏi bin có bao nhiêu điểm. Nhược điểm: không chính xác "nhìn thấy" một cụm điểm mà đặt một ranh giới bin; và: các thùng của "chiều rộng theo chiều dọc không đổi" thực sự là (không gian) nhỏ hơn khi bạn di chuyển về phía bắc.
• Sử dụng thư viện "Shapely" cho Python. Thực hiện theo ví dụ cổ phiếu của nó cho "điểm đệm", và làm một liên minh cascaded của bộ đệm. Hãy tìm những quả bóng trên một khu vực nhất định, hoặc "chứa" một số điểm ban đầu nhất định. Lưu ý rằng Shapely không thực chất là "geo-savy", vì vậy bạn sẽ phải thêm sửa đổi nếu bạn cần chúng.
• Sử dụng DB thực với xử lý không gian. MySQL, Oracle, Postgres (với PostGIS), MSSQL tất cả (tôi nghĩ) có các kiểu dữ liệu "Hình học" và "Địa lý", và bạn có thể thực hiện các truy vấn không gian trên chúng (từ các kịch bản Python của bạn).

Mỗi phòng trong số này có chi phí khác nhau về đô la và thời gian (trong đường cong học tập) ... và các mức độ khác nhau về độ chính xác không gian địa lý. Bạn phải chọn những gì phù hợp với ngân sách và/hoặc yêu cầu của bạn.

Answer 4

tôi đang làm điều này một cách thường xuyên bằng cách đầu tiên tạo ra một ma trận khoảng cách và sau đó chạy clustering trên nó. Đây là mã của tôi.

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) 

Tôi không chắc liệu nó có giải quyết được hoàn toàn vấn đề của bạn hay không. Bạn có thể muốn thử nghiệm với các k khác nhau, và cũng có thể thực hiện một nhóm thứ hai của một số cụm đầu tiên trong trường hợp chúng quá lớn, như thể bạn có một điểm ở Minnesota và một nghìn ở California. Khi bạn có điểm.to.group $ group, bạn có thể nhận các hộp giới hạn bằng cách tìm vĩ độ tối đa và tối cho mỗi nhóm.

Nếu bạn muốn X là 20, và bạn có 18 điểm ở New York và 22 ở Dallas, bạn phải quyết định xem bạn muốn một hộp nhỏ và một hộp thực sự lớn (20 điểm), nếu tốt hơn có hộp Dallas bao gồm 22 điểm, hoặc nếu bạn muốn chia 22 điểm ở Dallas thành hai nhóm. Clustering dựa trên khoảng cách có thể được tốt trong một số các trường hợp này. Nhưng tất nhiên phụ thuộc vào lý do tại sao bạn muốn nhóm các điểm.

/Chris

Answer 5

nếu bạn sử dụng cân đối, bạn có thể mở rộng của tôi cluster_points function để trả lại khung giới hạn của cụm thông qua .bounds tài sản của hình học cân đối, ví dụ như thế này:

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