# 安裝必要的庫
!pip install geopandas matplotlib contextily rtree shapely
# 導(dǎo)入庫
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import contextily as ctx
from shapely.geometry import Point, LineString, Polygon

# 加載行政區(qū)劃數(shù)據(jù)
districts = gpd.read_file("districts.geojson")
# 創(chuàng)建POI點(diǎn)數(shù)據(jù)
poi_data = {
    'name': ['咖啡店A', '餐廳B', '商場C', '學(xué)校D', '醫(yī)院E', '公園F', '圖書館G', '超市H'],
    'type': ['餐飲', '餐飲', '購物', '教育', '醫(yī)療', '休閑', '文化', '購物'],
    'latitude': [39.91, 39.92, 39.93, 39.89, 39.91, 39.88, 39.90, 39.94],
    'longitude': [116.41, 116.42, 116.39, 116.38, 116.44, 116.43, 116.41, 116.40]
}
# 創(chuàng)建點(diǎn)幾何
geometry = [Point(xy) for xy in zip(poi_data['longitude'], poi_data['latitude'])]
pois = gpd.GeoDataFrame(poi_data, geometry=geometry, crs="EPSG:4326")

# 確保坐標(biāo)系統(tǒng)一致
districts = districts.to_crs(pois.crs)
# 執(zhí)行空間連接 - 點(diǎn)位于哪個(gè)區(qū)域
joined_data = gpd.sjoin(pois, districts, how="left", predicate="within")
print(f"連接后的數(shù)據(jù)包含 {len(joined_data)} 行，其中包括點(diǎn)位信息和所在區(qū)域")

# 繪制底圖
fig, ax = plt.subplots(figsize=(12, 10))
# 繪制區(qū)域邊界
districts.boundary.plot(ax=ax, color='black', linewidth=1)
# 使用分類變量繪制不同類型的POI
pois.plot(ax=ax, column='type', cmap='viridis', 
          legend=True, markersize=100, categorical=True)
# 添加標(biāo)題
plt.title('POI分布與行政區(qū)劃', fontsize=15)
# 添加背景底圖
ctx.add_basemap(ax, crs=pois.crs.to_string(), source=ctx.providers.CartoDB.Positron)
plt.tight_layout()
plt.show()

# 使用不同空間關(guān)系進(jìn)行連接
within_join = gpd.sjoin(pois, districts, how="inner", predicate="within")
intersects_join = gpd.sjoin(pois, districts, how="inner", predicate="intersects")
print(f"within連接結(jié)果: {len(within_join)} 行")
print(f"intersects連接結(jié)果: {len(intersects_join)} 行")
# 對(duì)于點(diǎn)數(shù)據(jù)，"within"和"intersects"通常產(chǎn)生相同結(jié)果
# 但對(duì)于線或面數(shù)據(jù)則不同

# 篩選特定類型的POI并確定它們所在的行政區(qū)
restaurant_pois = pois[pois['type'] == '餐飲']
restaurant_districts = gpd.sjoin(restaurant_pois, districts, how="left", predicate="within")
# 統(tǒng)計(jì)每個(gè)行政區(qū)內(nèi)餐飲POI的數(shù)量
poi_count_by_district = restaurant_districts.groupby('district_name').size()
print("各行政區(qū)餐飲POI數(shù)量統(tǒng)計(jì):")
print(poi_count_by_district)

# 為每個(gè)POI創(chuàng)建500米緩沖區(qū)
# 首先將數(shù)據(jù)轉(zhuǎn)換為投影坐標(biāo)系統(tǒng)以使用米為單位
pois_projected = pois.to_crs(epsg=3857)  # Web Mercator投影
pois_buffered = pois_projected.copy()
pois_buffered['geometry'] = pois_projected.buffer(500)  # 500米緩沖區(qū)
# 轉(zhuǎn)回原始坐標(biāo)系進(jìn)行顯示
pois_buffered = pois_buffered.to_crs(pois.crs)
# 繪制緩沖區(qū)
fig, ax = plt.subplots(figsize=(12, 10))
districts.boundary.plot(ax=ax, color='black', linewidth=1)
pois_buffered.plot(ax=ax, alpha=0.5, column='type', categorical=True)
pois.plot(ax=ax, color='red', markersize=20)
ctx.add_basemap(ax, crs=pois.crs.to_string(), source=ctx.providers.CartoDB.Positron)
plt.title('POI點(diǎn)位500米服務(wù)范圍分析', fontsize=15)
plt.tight_layout()
plt.show()

# 將緩沖區(qū)與行政區(qū)劃進(jìn)行交集運(yùn)算
districts_projected = districts.to_crs(epsg=3857)
pois_buffered_projected = pois_buffered.to_crs(epsg=3857)
# 選擇一個(gè)特定區(qū)域進(jìn)行示例
target_district = districts_projected.iloc[0]
# 計(jì)算該區(qū)域與各POI緩沖區(qū)的交集
intersections = []
for idx, poi_buffer in pois_buffered_projected.iterrows():
    if poi_buffer.geometry.intersects(target_district.geometry):
        intersection = poi_buffer.geometry.intersection(target_district.geometry)
        intersections.append({
            'poi_name': poi_buffer['name'],
            'poi_type': poi_buffer['type'],
            'geometry': intersection,
            'area': intersection.area  # 平方米
        })
# 創(chuàng)建交集GeoDataFrame
intersections_gdf = gpd.GeoDataFrame(intersections, crs=pois_buffered_projected.crs)
# 計(jì)算覆蓋率
district_area = target_district.geometry.area
total_coverage_area = gpd.GeoSeries([geom['geometry'] for geom in intersections]).unary_union.area
coverage_ratio = total_coverage_area / district_area * 100
print(f"區(qū)域名稱: {target_district['district_name']}")
print(f"區(qū)域總面積: {district_area/1e6:.2f} 平方公里")
print(f"POI服務(wù)覆蓋面積: {total_coverage_area/1e6:.2f} 平方公里")
print(f"覆蓋率: {coverage_ratio:.2f}%")

# 計(jì)算所有POI之間的距離矩陣（以米為單位）
def calculate_distance_matrix(gdf):
    n = len(gdf)
    dist_matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(i+1, n):
            dist = gdf.iloc[i].geometry.distance(gdf.iloc[j].geometry)
            dist_matrix[i, j] = dist
            dist_matrix[j, i] = dist
    return pd.DataFrame(dist_matrix, 
                        index=gdf.index, 
                        columns=gdf.index)
# 計(jì)算距離矩陣
distance_matrix = calculate_distance_matrix(pois_projected)
# 查找每個(gè)POI的最近鄰
nearest_neighbors = {}
for idx in pois.index:
    distances = distance_matrix[idx].drop(idx)  # 排除自身
    nearest = distances.idxmin()
    nearest_neighbors[pois.loc[idx, 'name']] = {
        'nearest_poi': pois.loc[nearest, 'name'],
        'distance_meters': distances.min()
    }
print("每個(gè)POI的最近鄰設(shè)施:")
for poi, info in nearest_neighbors.items():
    print(f"{poi} -> {info['nearest_poi']} (距離: {info['distance_meters']:.2f}米)")

# 使用KDE進(jìn)行熱點(diǎn)分析
from scipy.stats import gaussian_kde
# 提取POI坐標(biāo)
x = pois_projected.geometry.x
y = pois_projected.geometry.y
positions = np.vstack([x, y])
# 計(jì)算KDE
kernel = gaussian_kde(positions)
# 創(chuàng)建網(wǎng)格進(jìn)行評(píng)估
x_grid, y_grid = np.mgrid[x.min():x.max():100j, y.min():y.max():100j]
positions_grid = np.vstack([x_grid.ravel(), y_grid.ravel()])
# 計(jì)算密度
z = kernel(positions_grid).reshape(x_grid.shape)
# 可視化結(jié)果
fig, ax = plt.subplots(figsize=(12, 10))
districts_projected.boundary.plot(ax=ax, color='black', linewidth=1)
contour = ax.contourf(x_grid, y_grid, z, cmap='viridis', levels=20, alpha=0.75)
pois_projected.plot(ax=ax, color='red', markersize=20)
plt.colorbar(contour, ax=ax, label='密度')
plt.title('POI熱點(diǎn)密度分析', fontsize=15)
plt.tight_layout()
plt.show()

from sklearn.cluster import DBSCAN
# 提取坐標(biāo)
coords = np.vstack((pois_projected.geometry.x, pois_projected.geometry.y)).T
# 使用DBSCAN進(jìn)行聚類
clustering = DBSCAN(eps=1000, min_samples=2).fit(coords)  # 1000米范圍內(nèi)至少2個(gè)點(diǎn)形成聚類
# 將聚類結(jié)果添加到原始數(shù)據(jù)
pois_projected['cluster'] = clustering.labels_
# 可視化聚類結(jié)果
fig, ax = plt.subplots(figsize=(12, 10))
districts_projected.boundary.plot(ax=ax, color='black', linewidth=1)
pois_projected.plot(ax=ax, column='cluster', categorical=True, 
                   legend=True, markersize=100)
plt.title('POI空間聚類分析 (DBSCAN)', fontsize=15)
plt.tight_layout()
plt.show()
# 統(tǒng)計(jì)每個(gè)聚類的POI數(shù)量
cluster_stats = pois_projected.groupby('cluster').size()
print("各聚類POI數(shù)量:")
print(cluster_stats)

# 注意：這需要預(yù)先準(zhǔn)備道路網(wǎng)絡(luò)數(shù)據(jù)和networkx庫
import networkx as nx
# 加載道路網(wǎng)絡(luò)數(shù)據(jù)
roads = gpd.read_file("roads.geojson")
# 創(chuàng)建網(wǎng)絡(luò)圖
G = nx.Graph()
# 從道路數(shù)據(jù)中構(gòu)建網(wǎng)絡(luò)
for idx, road in roads.iterrows():
    coords = list(road.geometry.coords)
    for i in range(len(coords) - 1):
        G.add_edge(coords[i], coords[i+1], 
                  weight=road.geometry.length,
                  road_id=road['id'])
# 找到距離每個(gè)POI最近的道路節(jié)點(diǎn)
nearest_nodes = {}
for idx, poi in pois_projected.iterrows():
    min_dist = float('inf')
    nearest_node = None
    for node in G.nodes():
        dist = poi.geometry.distance(Point(node))
        if dist < min_dist:
            min_dist = dist
            nearest_node = node
    nearest_nodes[idx] = nearest_node
# 計(jì)算兩個(gè)POI之間的最短路徑
origin_idx = 0  # 起點(diǎn)POI索引
dest_idx = 3    # 終點(diǎn)POI索引
origin_node = nearest_nodes[origin_idx]
dest_node = nearest_nodes[dest_idx]
try:
    # 計(jì)算最短路徑
    shortest_path = nx.shortest_path(G, origin_node, dest_node, weight='weight')
    path_length = nx.shortest_path_length(G, origin_node, dest_node, weight='weight')
    print(f"從 {pois.loc[origin_idx, 'name']} 到 {pois.loc[dest_idx, 'name']} 的最短路徑長度: {path_length:.2f}米")
    # 可視化路徑
    path_coords = shortest_path
    path_line = LineString(path_coords)
    path_gdf = gpd.GeoDataFrame({'geometry': [path_line]}, crs=pois_projected.crs)
    fig, ax = plt.subplots(figsize=(12, 10))
    roads.plot(ax=ax, color='gray', linewidth=1)
    path_gdf.plot(ax=ax, color='red', linewidth=3)
    pois_projected.iloc[[origin_idx, dest_idx]].plot(ax=ax, color='blue', markersize=100)
    plt.title(f'從 {pois.loc[origin_idx, "name"]} 到 {pois.loc[dest_idx, "name"]} 的最短路徑', fontsize=15)
    plt.tight_layout()
    plt.show()
except nx.NetworkXNoPath:
    print(f"無法找到從 {pois.loc[origin_idx, 'name']} 到 {pois.loc[dest_idx, 'name']} 的路徑")

# 確保已安裝rtree庫
# pip install rtree
# 創(chuàng)建空間索引
districts = districts.set_index('district_id', drop=False)
districts_sindex = districts.sindex
# 使用空間索引加速查詢
def fast_sjoin_points_to_polygons(points_gdf, polygons_gdf):
    # 創(chuàng)建結(jié)果列表
    results = []
    # 遍歷點(diǎn)
    for idx, point in points_gdf.iterrows():
        # 使用空間索引查找可能相交的多邊形
        possible_matches_idx = list(polygons_gdf.sindex.intersection(point.geometry.bounds))
        possible_matches = polygons_gdf.iloc[possible_matches_idx]
        # 精確檢查點(diǎn)是否在多邊形內(nèi)
        precise_matches = possible_matches[possible_matches.contains(point.geometry)]
        # 如果找到匹配項(xiàng)
        if len(precise_matches) > 0:
            for match_idx, match in precise_matches.iterrows():
                result = point.to_dict()
                result.update({
                    'district_id': match['district_id'],
                    'district_name': match['district_name']
                })
                results.append(result)
    if results:
        return gpd.GeoDataFrame(results, crs=points_gdf.crs)
    else:
        return gpd.GeoDataFrame([], crs=points_gdf.crs)
# 使用優(yōu)化的函數(shù)進(jìn)行空間連接
import time
start_time = time.time()
optimized_result = fast_sjoin_points_to_polygons(pois, districts)
end_time = time.time()
print(f"優(yōu)化后的空間連接耗時(shí): {end_time - start_time:.4f}秒")
# 對(duì)比標(biāo)準(zhǔn)sjoin
start_time = time.time()
standard_result = gpd.sjoin(pois, districts, how="left", predicate="within")
end_time = time.time()
print(f"標(biāo)準(zhǔn)空間連接耗時(shí): {end_time - start_time:.4f}秒")

from joblib import Parallel, delayed
# 并行處理緩沖區(qū)計(jì)算
def parallel_buffer(gdf, buffer_distance, n_jobs=-1):
    # 定義單個(gè)幾何體緩沖區(qū)計(jì)算函數(shù)
    def buffer_geom(geom):
        return geom.buffer(buffer_distance)
    # 并行執(zhí)行
    buffered_geoms = Parallel(n_jobs=n_jobs)(
        delayed(buffer_geom)(geom) for geom in gdf.geometry
    )
    # 創(chuàng)建新的GeoDataFrame
    result = gdf.copy()
    result['geometry'] = buffered_geoms
    return result
# 使用并行處理創(chuàng)建緩沖區(qū)
start_time = time.time()
pois_buffered_parallel = parallel_buffer(pois_projected, 500)
end_time = time.time()
print(f"并行緩沖區(qū)計(jì)算耗時(shí): {end_time - start_time:.4f}秒")
# 對(duì)比標(biāo)準(zhǔn)處理
start_time = time.time()
pois_buffered_standard = pois_projected.copy()
pois_buffered_standard['geometry'] = pois_projected.buffer(500)
end_time = time.time()
print(f"標(biāo)準(zhǔn)緩沖區(qū)計(jì)算耗時(shí): {end_time - start_time:.4f}秒")

# 假設(shè)有人口數(shù)據(jù)
population = gpd.read_file("population_grid.geojson")
# 選擇交通設(shè)施POI
transport_pois = pois[pois['type'].isin(['公交站', '地鐵站', '火車站'])]
# 創(chuàng)建交通設(shè)施緩沖區(qū)（500米步行距離）
transport_pois_projected = transport_pois.to_crs(epsg=3857)
transport_buffers = transport_pois_projected.copy()
transport_buffers['geometry'] = transport_pois_projected.buffer(500)
# 與人口網(wǎng)格進(jìn)行空間連接
population_projected = population.to_crs(epsg=3857)
population_covered = gpd.sjoin(population_projected, transport_buffers.to_crs(population_projected.crs), 
                              how="inner", predicate="intersects")
# 計(jì)算覆蓋的總?cè)丝?
total_population = population_projected['population'].sum()
covered_population = population_covered['population'].sum()
coverage_percentage = (covered_population / total_population) * 100
print(f"總?cè)丝? {total_population}")
print(f"公共交通500米覆蓋人口: {covered_population}")
print(f"覆蓋率: {coverage_percentage:.2f}%")

# 假設(shè)我們有以下數(shù)據(jù)層
# - 人口密度
# - 競爭對(duì)手位置
# - 交通可達(dá)性
# - 商業(yè)中心
# 創(chuàng)建一個(gè)網(wǎng)格評(píng)價(jià)系統(tǒng)
from shapely.geometry import box
# 創(chuàng)建研究區(qū)域網(wǎng)格
xmin, ymin, xmax, ymax = districts_projected.total_bounds
cell_size = 500  # 500米網(wǎng)格
rows = int((ymax - ymin) / cell_size)
cols = int((xmax - xmin) / cell_size)
grid_cells = []
for i in range(cols):
    for j in range(rows):
        x1 = xmin + i * cell_size
        y1 = ymin + j * cell_size
        x2 = xmin + (i + 1) * cell_size
        y2 = ymin + (j + 1) * cell_size
        grid_cells.append(box(x1, y1, x2, y2))
grid = gpd.GeoDataFrame({'geometry': grid_cells}, crs=districts_projected.crs)
# 計(jì)算各因素得分
# 1. 人口密度得分
grid['pop_score'] = np.random.randint(1, 10, size=len(grid))  # 模擬數(shù)據(jù)
# 2. 競爭對(duì)手距離得分
competitors = pois_projected[pois_projected['type'] == '餐飲']
grid['comp_dist'] = grid.geometry.apply(
    lambda g: min([g.distance(comp) for comp in competitors.geometry]) if len(competitors) > 0 else 0
)
grid['comp_score'] = grid['comp_dist'] / grid['comp_dist'].max() * 10
# 3. 交通可達(dá)性得分
transit = pois_projected[pois_projected['type'].isin(['公交站', '地鐵站'])]
grid['transit_dist'] = grid.geometry.apply(
    lambda g: min([g.distance(t) for t in transit.geometry]) if len(transit) > 0 else float('inf')
)
grid['transit_score'] = 10 - (grid['transit_dist'] / grid['transit_dist'].max() * 10)
# 4. 商業(yè)中心距離得分
business_centers = pois_projected[pois_projected['type'] == '購物']
grid['biz_dist'] = grid.geometry.apply(
    lambda g: min([g.distance(bc) for bc in business_centers.geometry]) if len(business_centers) > 0 else float('inf')
)
grid['biz_score'] = 10 - (grid['biz_dist'] / grid['biz_dist'].max() * 10)
# 計(jì)算綜合得分
grid['total_score'] = (
    grid['pop_score'] * 0.4 + 
    grid['comp_score'] * 0.2 + 
    grid['transit_score'] * 0.2 + 
    grid['biz_score'] * 0.2
)
# 可視化結(jié)果
fig, ax = plt.subplots(figsize=(12, 10))
districts_projected.boundary.plot(ax=ax, color='black', linewidth=1)
grid.plot(ax=ax, column='total_score', cmap='RdYlGn', alpha=0.7, 
          legend=True, legend_kwds={'label': '選址適宜性得分'})
pois_projected.plot(ax=ax, color='blue', markersize=10)
plt.title('商業(yè)選址適宜性評(píng)價(jià)', fontsize=15)
plt.tight_layout()
plt.show()
# 找出得分最高的位置
best_locations = grid.sort_values('total_score', ascending=False).head(3)
print("最適宜開設(shè)新店鋪的3個(gè)位置:")
for idx, loc in best_locations.iterrows():
    centroid = loc.geometry.centroid
    print(f"位置 {idx}: 坐標(biāo) ({centroid.x}, {centroid.y}), 得分: {loc['total_score']:.2f}")