In [1]:
import pandas as pd
import geopandas as gpd
import numpy as np
import os
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = 'Microsoft YaHei'
import rasterio
import json
from shapely.geometry import Point
import contextily as ctx
In [ ]:
youbike_sta = pd.read_csv('./dataset/youbike_station.csv')
youbike_sta = gpd.GeoDataFrame(youbike_sta,
geometry=gpd.points_from_xy(youbike_sta['lon'], youbike_sta['lat']),crs='EPSG:4326').to_crs('EPSG:3826')
spot = gpd.read_file('./dataset/Scenic_spot_C_f.shp').to_crs('EPSG:3826')
spot = spot[spot['Region'] == '臺北市']
school = pd.read_csv('./dataset/112年_各級學校名錄_240828.csv')
school = school[school['縣市別'] == '臺北市']
school = gpd.GeoDataFrame(school, geometry=gpd.points_from_xy(school['經度'], school['緯度']), crs='EPSG:4326').to_crs('EPSG:3826')
FF = gpd.read_file('./dataset/TPE_FastFood.shp', encoding='big5', crs='EPSG:3826')
TPE = gpd.read_file('./dataset/Taipei_Vill.shp', encoding='big5', crs='EPSG:3826')
tif = rasterio.open('./dataset/taipei_s2.tif')
r = tif.read(3).astype(float)
nir = tif.read(7).astype(float)
NDVI = (nir - r) / (nir + r + 1e-8)
In [10]:
def classify_schools(school_gdf):
"""分類學校層級"""
# 使用學校級別欄位進行分類
school_gdf['school_type'] = 'others'
# 國小相關
elementary_mask = school_gdf['學校級別'].isin(['國民小學', '附設國民小學'])
school_gdf.loc[elementary_mask, 'school_type'] = 'elementary'
# 國中高中相關
middle_high_mask = school_gdf['學校級別'].isin(['國民中學', '附設國民中學', '高級中等學校'])
school_gdf.loc[middle_high_mask, 'school_type'] = 'middle_high'
# 大學及專科相關
university_mask = school_gdf['學校級別'].isin(['大學', '專科學校', '學院'])
school_gdf.loc[university_mask, 'school_type'] = 'university'
return school_gdf
def create_fastfood_restriction_zones(school_gdf, ff_gdf):
"""建立速食店管制區域"""
# 轉換至投影坐標系進行緩衝區計算
school_proj = school_gdf.to_crs('EPSG:3826')
ff_proj = ff_gdf.to_crs('EPSG:3826')
# 建立不同學校類型的緩衝區
elementary = school_proj[school_proj['school_type'] == 'elementary']
middle_high = school_proj[school_proj['school_type'] == 'middle_high']
all_buffers = []
# 國小500公尺緩衝區
if len(elementary) > 0:
elementary_buffer = elementary.geometry.buffer(500)
elem_gdf = gpd.GeoDataFrame({
'school_type': ['elementary'] * len(elementary_buffer),
'buffer_distance': [500] * len(elementary_buffer)
}, geometry=elementary_buffer, crs='EPSG:3826')
all_buffers.append(elem_gdf)
# 國高中300公尺緩衝區
if len(middle_high) > 0:
middle_high_buffer = middle_high.geometry.buffer(300)
mh_gdf = gpd.GeoDataFrame({
'school_type': ['middle_high'] * len(middle_high_buffer),
'buffer_distance': [300] * len(middle_high_buffer)
}, geometry=middle_high_buffer, crs='EPSG:3826')
all_buffers.append(mh_gdf)
# 合併所有管制區域
if all_buffers:
all_buffers = pd.concat(all_buffers, ignore_index=True)
else:
# 如果沒有緩衝區，建立空的 GeoDataFrame
all_buffers = gpd.GeoDataFrame({
'school_type': [],
'buffer_distance': []
}, geometry=[], crs='EPSG:3826')
# 分析違規速食店
violation_analysis = {}
for school_type in ['elementary', 'middle_high']:
buffers = all_buffers[all_buffers['school_type'] == school_type]
if len(buffers) > 0:
buffer_union = buffers.geometry.unary_union
# 確保 buffer_union 不是 None
if buffer_union is not None:
# 找出緩衝區內的速食店
violations = ff_proj[ff_proj.geometry.within(buffer_union)]
else:
violations = gpd.GeoDataFrame(columns=ff_proj.columns, crs=ff_proj.crs)
else:
violations = gpd.GeoDataFrame(columns=ff_proj.columns, crs=ff_proj.crs)
violation_analysis[school_type] = {
'count': len(violations),
'locations': violations
}
return all_buffers.to_crs('EPSG:4326'), violation_analysis
def analyze_green_accessibility(school_gdf, spot_gdf, youbike_gdf):
"""分析綠地可及性"""
school_proj = school_gdf.to_crs('EPSG:3826')
spot_proj = spot_gdf.to_crs('EPSG:3826')
youbike_proj = youbike_gdf.to_crs('EPSG:3826')
accessibility_results = []
for idx, school in school_proj.iterrows():
school_point = school.geometry
# 步行可及性分析 (100, 300, 500公尺)
walking_distances = [100, 300, 500]
walking_access = {}
for dist in walking_distances:
buffer = school_point.buffer(dist)
accessible_spots = spot_proj[spot_proj.geometry.within(buffer)]
walking_access[f'walk_{dist}m'] = len(accessible_spots)
# 騎行可及性分析 (1000, 2000公尺至YouBike + 300公尺步行至綠地)
cycling_access = {}
for bike_dist in [1000, 2000]:
bike_buffer = school_point.buffer(bike_dist)
accessible_youbikes = youbike_proj[youbike_proj.geometry.within(bike_buffer)]
# 從YouBike站點300公尺內的綠地
green_from_bike = 0
for _, bike in accessible_youbikes.iterrows():
green_buffer = bike.geometry.buffer(300)
nearby_greens = spot_proj[spot_proj.geometry.within(green_buffer)]
green_from_bike += len(nearby_greens)
cycling_access[f'bike_{bike_dist}m'] = green_from_bike
# 整合結果
result = {
'school_id': idx,
'school_name': school['學校名稱'],
'school_type': school['school_type'],
**walking_access,
**cycling_access
}
accessibility_results.append(result)
return pd.DataFrame(accessibility_results)
def extract_ndvi_for_schools(school_gdf, ndvi_array, transform):
"""提取學校周邊NDVI值"""
# 確保學校資料與NDVI使用相同的坐標系
school_proj = school_gdf.to_crs(tif.crs) # 使用tif檔案的坐標系
ndvi_results = []
for idx, school in school_proj.iterrows():
# 提取學校位置的NDVI值
try:
# 轉換座標至raster索引
row, col = rasterio.transform.rowcol(transform, school.geometry.x, school.geometry.y)
if 0 <= row < ndvi_array.shape[0] and 0 <= col < ndvi_array.shape[1]:
ndvi_value = ndvi_array[row, col]
else:
ndvi_value = np.nan
ndvi_results.append({
'school_id': idx,
'ndvi_value': ndvi_value
})
except Exception as e:
print(f"Error processing school {idx}: {e}")
ndvi_results.append({
'school_id': idx,
'ndvi_value': np.nan
})
return pd.DataFrame(ndvi_results)
def create_health_zones(accessibility_df, violation_analysis, ndvi_df):
"""建立健康區域評分系統 - 分析所有學校"""
# 合併所有分析結果
health_df = accessibility_df.copy()
health_df = health_df.merge(ndvi_df, on='school_id', how='left')
# 計算健康指標分數 (0-100)
health_df['green_access_score'] = (
health_df['walk_100m'] * 3 +
health_df['walk_300m'] * 2 +
health_df['walk_500m'] * 1 +
health_df['bike_1000m'] * 0.5 +
health_df['bike_2000m'] * 0.3
) / 10 * 100
# NDVI分數標準化 (處理NaN值)
ndvi_min = health_df['ndvi_value'].min()
ndvi_max = health_df['ndvi_value'].max()
if pd.notna(ndvi_min) and pd.notna(ndvi_max) and ndvi_max != ndvi_min:
health_df['ndvi_score'] = ((health_df['ndvi_value'] - ndvi_min) /
(ndvi_max - ndvi_min) * 100)
else:
health_df['ndvi_score'] = 50 # 預設分數
# 速食店風險分數計算 - 針對所有學校
health_df['fastfood_risk_score'] = 100 # 預設滿分
# 對有速食店管制的學校類型進行扣分
for school_type in ['elementary', 'middle_high']:
if school_type in violation_analysis:
violation_count = violation_analysis[school_type]['count']
# 對該類型學校進行風險評估
school_mask = health_df['school_type'] == school_type
if violation_count > 0:
# 根據違規速食店數量調整風險分數
risk_penalty = min(violation_count * 10, 50) # 每家違規店扣10分，最多扣50分
health_df.loc[school_mask, 'fastfood_risk_score'] -= risk_penalty
# 綜合健康指標 - 對所有學校適用
health_df['overall_health_score'] = (
health_df['green_access_score'] * 0.4 +
health_df['ndvi_score'] * 0.3 +
health_df['fastfood_risk_score'] * 0.3
)
# 分級 - 適用於所有學校
health_df['health_zone'] = pd.cut(health_df['overall_health_score'],
bins=[0, 30, 60, 100],
labels=['需改善', '一般', '優良'])
return health_df
# === 3. 執行分析 ===
# 檢查學校分類結果
school = classify_schools(school)
print("學校分類結果:")
print(school['school_type'].value_counts())
print("\n各級學校數量:")
print(school['學校級別'].value_counts())
# 分類學校
school = classify_schools(school)
# 速食店管制分析
restriction_zones, violations = create_fastfood_restriction_zones(school, FF)
# 綠地可及性分析
accessibility_results = analyze_green_accessibility(school, spot, youbike_sta)
# NDVI分析
ndvi_results = extract_ndvi_for_schools(school, NDVI, tif.transform)
# 健康區域劃分
health_zones = create_health_zones(accessibility_results, violations, ndvi_results)
print("分析執行完成！")
print(f"分析學校數量: {len(school)}")
print(f"國小數量: {len(school[school['school_type'] == 'elementary'])}")
print(f"國高中數量: {len(school[school['school_type'] == 'middle_high'])}")
# === 4. 視覺化分析結果 ===
fig, axes = plt.subplots(2, 2, figsize=(20, 16))
# 轉換所有資料到相同的坐標系統進行視覺化
TPE_viz = TPE.to_crs('EPSG:4326')
FF_viz = FF.to_crs('EPSG:4326')
school_viz = school.to_crs('EPSG:4326')
spot_viz = spot.to_crs('EPSG:4326')
# 安全的轉換 restriction_zones
if len(restriction_zones) > 0 and not restriction_zones.is_empty.all():
restriction_zones_viz = restriction_zones.copy()
# 確保幾何有效性
restriction_zones_viz = restriction_zones_viz[~restriction_zones_viz.is_empty]
restriction_zones_viz = restriction_zones_viz[restriction_zones_viz.geometry.is_valid]
if len(restriction_zones_viz) > 0:
# 檢查邊界值是否有效
bounds = restriction_zones_viz.total_bounds
if np.isfinite(bounds).all():
pass # 邊界值正常
else:
restriction_zones_viz = gpd.GeoDataFrame() # 清空無效資料
else:
restriction_zones_viz = gpd.GeoDataFrame()
else:
restriction_zones_viz = gpd.GeoDataFrame()
# 4.1 速食店管制區域圖
ax1 = axes[0, 0]
TPE_viz.plot(ax=ax1, color='lightgray', alpha=0.5)
# 安全地繪製管制區域
if len(restriction_zones_viz) > 0:
try:
restriction_zones_viz.plot(ax=ax1, color='red', alpha=0.3)
except (ValueError, AttributeError) as e:
print(f"警告：無法繪製管制區域 - {e}")
FF_viz.plot(ax=ax1, color='orange', markersize=5, label='速食店')
school_viz.plot(ax=ax1, color='blue', markersize=3, label='學校')
ax1.set_title('學校周邊速食店管制區域')
ax1.legend()
ax1.set_aspect('equal') # 手動設定比例
ax1.axis('off')
# 4.2 綠地可及性熱圖
ax2 = axes[0, 1]
try:
school_with_health = school_viz.merge(health_zones, left_index=True, right_on='school_id')
school_with_health = gpd.GeoDataFrame(school_with_health, geometry=school_with_health.geometry, crs='EPSG:4326')
TPE_viz.plot(ax=ax2, color='lightgray', alpha=0.3)
if 'green_access_score' in school_with_health.columns:
# 過濾掉無效值
valid_data = school_with_health.dropna(subset=['green_access_score'])
if len(valid_data) > 0:
valid_data.plot(ax=ax2, column='green_access_score',
cmap='Greens', markersize=20, legend=True,
legend_kwds={'shrink': 0.5})
spot_viz.plot(ax=ax2, color='green', markersize=1, alpha=0.7)
except Exception as e:
print(f"警告：繪製綠地可及性圖時發生錯誤 - {e}")
TPE_viz.plot(ax=ax2, color='lightgray', alpha=0.3)
ax2.set_title('學校綠地可及性分析')
ax2.set_aspect('equal')
ax2.axis('off')
# 4.3 NDVI分布圖
ax3 = axes[1, 0]
try:
TPE_viz.plot(ax=ax3, color='lightgray', alpha=0.3)
if 'school_with_health' in locals() and 'ndvi_score' in school_with_health.columns:
# 過濾掉無效值
valid_ndvi = school_with_health.dropna(subset=['ndvi_score'])
if len(valid_ndvi) > 0:
valid_ndvi.plot(ax=ax3, column='ndvi_score',
cmap='RdYlGn', markersize=20, legend=True,
legend_kwds={'shrink': 0.5})
except Exception as e:
print(f"警告：繪製NDVI圖時發生錯誤 - {e}")
TPE_viz.plot(ax=ax3, color='lightgray', alpha=0.3)
ax3.set_title('學校周邊綠度分析(NDVI)')
ax3.set_aspect('equal')
ax3.axis('off')
# 4.4 綜合健康區域圖
ax4 = axes[1, 1]
try:
TPE_viz.plot(ax=ax4, color='lightgray', alpha=0.3)
colors = {'需改善': 'red', '一般': 'yellow', '優良': 'green'}
if 'school_with_health' in locals() and 'health_zone' in health_zones.columns:
for zone in health_zones['health_zone'].unique():
if pd.notna(zone):
subset = school_with_health[school_with_health['health_zone'] == zone]
if len(subset) > 0:
subset.plot(ax=ax4, color=colors.get(zone, 'gray'), markersize=20,
label=f'{zone}區域', alpha=0.7)
except Exception as e:
print(f"警告：繪製健康區域圖時發生錯誤 - {e}")
TPE_viz.plot(ax=ax4, color='lightgray', alpha=0.3)
ax4.set_title('校園周邊綜合健康區域分級')
ax4.legend()
ax4.set_aspect('equal')
ax4.axis('off')
plt.tight_layout()
plt.show()
# ...existing code...
# === 5. 統計分析報告 ===
print("\n=== 台北市校園周邊健康生活分析報告 ===\n")
print("1. 速食店管制區域分析:")
print(f" - 國小500m緩衝區內違規速食店: {violations['elementary']['count']}家")
print(f" - 國高中300m緩衝區內違規速食店: {violations['middle_high']['count']}家")
print(" - 大學院校：無速食店管制限制")
print("\n2. 綠地可及性統計:")
accessibility_stats = accessibility_results.describe()
print(f" - 平均步行500m內綠地數量: {accessibility_stats.loc['mean', 'walk_500m']:.2f}")
print(f" - 平均騎行1000m+步行300m可達綠地數量: {accessibility_stats.loc['mean', 'bike_1000m']:.2f}")
print("\n3. 健康區域分布（涵蓋所有學校類型）:")
if 'health_zone' in health_zones.columns:
zone_counts = health_zones['health_zone'].value_counts()
for zone, count in zone_counts.items():
percentage = count / len(health_zones) * 100
print(f" - {zone}區域: {count}所學校 ({percentage:.1f}%)")
# 各學校類型的健康區域分布
print("\n各學校類型健康區域分布:")
school_type_mapping = {
'elementary': '國小',
'middle_high': '國高中',
'university': '大學院校',
'others': '其他'
}
for school_type in health_zones['school_type'].unique():
if pd.notna(school_type):
type_name = school_type_mapping.get(school_type, school_type)
type_data = health_zones[health_zones['school_type'] == school_type]
type_zones = type_data['health_zone'].value_counts()
print(f" {type_name}:")
for zone, count in type_zones.items():
percentage = count / len(type_data) * 100
print(f" - {zone}: {count}所 ({percentage:.1f}%)")
print("\n4. 政策建議:")
print(" - 國小及國高中：優先改善'需改善'區域的綠地建設和YouBike站點配置")
print(" - 國小及國高中：加強管制違規速食店的設立")
print(" - 大學院校：著重綠地可及性和環境品質提升")
print(" - 在綠地不足區域增設小型公園或綠化空間")
print(" - 建立校園健康環境監測機制")
# 在進階分析中增加各學校類型的詳細統計
print("\n=== 進階分析結果 ===")
# 7.1 學校類型與健康區域交叉分析（包含所有類型）
if 'health_zone' in health_zones.columns:
print("\n7.1 學校類型與健康區域交叉分析:")
cross_analysis = pd.crosstab(health_zones['school_type'], health_zones['health_zone'])
# 重命名索引以便閱讀
cross_analysis.index = cross_analysis.index.map(school_type_mapping)
print(cross_analysis)
# 7.2 綠地可及性詳細統計（所有學校類型）
print("\n7.2 綠地可及性詳細統計:")
green_stats = accessibility_results.groupby('school_type').agg({
'walk_100m': ['mean', 'std'],
'walk_300m': ['mean', 'std'],
'walk_500m': ['mean', 'std'],
'bike_1000m': ['mean', 'std'],
'bike_2000m': ['mean', 'std']
}).round(2)
# 重命名索引
green_stats.index = green_stats.index.map(school_type_mapping)
print(green_stats)
# 7.4 各學校類型健康指標平均值比較
print("\n7.4 各學校類型健康指標平均值:")
health_by_type = health_zones.groupby('school_type').agg({
'green_access_score': 'mean',
'ndvi_score': 'mean',
'fastfood_risk_score': 'mean',
'overall_health_score': 'mean'
}).round(2)
# 重命名索引
health_by_type.index = health_by_type.index.map(school_type_mapping)
print(health_by_type)
學校分類結果: school_type middle_high 160 elementary 154 others 62 Name: count, dtype: int64 各級學校數量: 學校級別 國民小學 147 高級中等學校 70 國民中學 59 附設國民中學 31 國中小補校 28 大專校院 24 附設國民小學 7 特殊教育學校 4 宗教研修學院 2 空大及大專校院附設進修學校 2 軍警大專校院 2 Name: count, dtype: int64
C:\Users\Chang Tzu Cheng\AppData\Local\Temp\ipykernel_38544\1768818812.py:65: FutureWarning: `unary_union` returned None due to all-None GeoSeries. In future, `unary_union` will return 'GEOMETRYCOLLECTION EMPTY' instead. buffer_union = buffers.geometry.unary_union C:\Users\Chang Tzu Cheng\AppData\Local\Temp\ipykernel_38544\1768818812.py:65: FutureWarning: `unary_union` returned None due to all-None GeoSeries. In future, `unary_union` will return 'GEOMETRYCOLLECTION EMPTY' instead. buffer_union = buffers.geometry.unary_union
分析執行完成！ 分析學校數量: 376 國小數量: 154 國高中數量: 160
=== 台北市校園周邊健康生活分析報告 ===
1. 速食店管制區域分析:
- 國小500m緩衝區內違規速食店: 0家
- 國高中300m緩衝區內違規速食店: 0家
- 大學院校：無速食店管制限制
2. 綠地可及性統計:
- 平均步行500m內綠地數量: 2.59
- 平均騎行1000m+步行300m可達綠地數量: 40.53
3. 健康區域分布（涵蓋所有學校類型）:
- 優良區域: 83所學校 (22.1%)
- 一般區域: 21所學校 (5.6%)
- 需改善區域: 0所學校 (0.0%)
各學校類型健康區域分布:
其他:
- 優良: 18所 (29.0%)
- 一般: 1所 (1.6%)
- 需改善: 0所 (0.0%)
國高中:
- 優良: 35所 (21.9%)
- 一般: 8所 (5.0%)
- 需改善: 0所 (0.0%)
國小:
- 優良: 30所 (19.5%)
- 一般: 12所 (7.8%)
- 需改善: 0所 (0.0%)
4. 政策建議:
- 國小及國高中：優先改善'需改善'區域的綠地建設和YouBike站點配置
- 國小及國高中：加強管制違規速食店的設立
- 大學院校：著重綠地可及性和環境品質提升
- 在綠地不足區域增設小型公園或綠化空間
- 建立校園健康環境監測機制
=== 進階分析結果 ===
7.1 學校類型與健康區域交叉分析:
health_zone 一般 優良
school_type
國小 12 30
國高中 8 35
其他 1 18
7.2 綠地可及性詳細統計:
walk_100m walk_300m walk_500m bike_1000m \
mean std mean std mean std mean std
school_type
國小 0.06 0.25 0.95 1.61 2.49 3.61 41.16 48.75
國高中 0.04 0.21 0.78 1.25 2.52 2.99 39.24 46.39
其他 0.13 0.38 1.02 1.58 3.02 3.71 42.31 47.31
bike_2000m
mean std
school_type
國小 131.44 128.71
國高中 142.72 134.75
其他 138.85 135.15
7.4 各學校類型健康指標平均值:
green_access_score ndvi_score fastfood_risk_score \
school_type
國小 645.90 53.77 100.0
國高中 666.38 48.04 100.0
其他 682.45 51.81 100.0
overall_health_score
school_type
國小 304.49
國高中 310.96
其他 318.52
In [ ]: