GeoMaster 是做什麼的？
全面的地理空間科學技能，涵蓋 GIS、遙感探測、空間分析與地球觀測機器學習，橫跨 70 多個主題，並提供 8 種程式語言的 500 多個程式碼範例。
安裝
# 核心 Python 套件棧（建議使用 conda）
conda install -c conda-forge gdal rasterio fiona shapely pyproj geopandas
# 遙感探測與機器學習（rsgislib 僅在 conda-forge 上提供，不在 PyPI）
conda install -c conda-forge rsgislib
uv pip install torchgeo earthengine-api
uv pip install scikit-learn xgboost torch-geometric
# 網路與視覺化
uv pip install osmnx networkx folium keplergl
uv pip install cartopy contextily mapclassify
# 大數據與雲端
uv pip install xarray rioxarray dask-geopandas
uv pip install pystac-client planetary-computer
# 點雲
uv pip install laspy pylas open3d pdal
# 資料庫
conda install -c conda-forge postgis spatialite
快速入門
從 Sentinel-2 計算 NDVI
import rasterio
import numpy as np
with rasterio.open('sentinel2.tif') as src:
red = src.read(4).astype(float) # B04
nir = src.read(8).astype(float) # B08
ndvi = (nir - red) / (nir + red + 1e-8)
ndvi = np.nan_to_num(ndvi, nan=0)
profile = src.profile
profile.update(count=1, dtype=rasterio.float32)
with rasterio.open('ndvi.tif', 'w', **profile) as dst:
dst.write(ndvi.astype(rasterio.float32), 1)
使用 GeoPandas 進行空間分析
import geopandas as gpd
# 載入並確保 CRS 一致
zones = gpd.read_file('zones.geojson')
points = gpd.read_file('points.geojson')
if zones.crs != points.crs:
points = points.to_crs(zones.crs)
# 空間連接與統計
joined = gpd.sjoin(points, zones, how='inner', predicate='within')
stats = joined.groupby('zone_id').agg({
'value': ['count', 'mean', 'std', 'min', 'max']
}).round(2)
Google Earth Engine 時間序列
import ee
import pandas as pd
ee.Initialize(project='your-project')
roi = ee.Geometry.Point([-122.4, 37.7]).buffer(10000)
s2 = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterBounds(roi)
.filterDate('2020-01-01', '2023-12-31')
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)))
def add_ndvi(img):
return img.addBands(img.normalizedDifference(['B8', 'B4']).rename('NDVI'))
s2_ndvi = s2.map(add_ndvi)
def extract_series(image):
stats = image.reduceRegion(ee.Reducer.mean(), roi.centroid(), scale=10, maxPixels=1e9)
return ee.Feature(None, {'date': image.date().format('YYYY-MM-dd'), 'ndvi': stats.get('NDVI')})
series = s2_ndvi.map(extract_series).getInfo()
df = pd.DataFrame([f['properties'] for f in series['features']])
df['date'] = pd.to_datetime(df['date'])
核心概念
資料類型
| 類型 | 範例 | 函式庫 |
|---|---|---|
| 向量 | Shapefile, GeoJSON, GeoPackage | GeoPandas, Fiona, GDAL |
| 柵格 | GeoTIFF, NetCDF, COG | Rasterio, Xarray, GDAL |
| 點雲 | LAS, LAZ | Laspy, PDAL, Open3D |
座標系統
- EPSG:4326 (WGS 84) - 地理座標，經緯度，用於儲存
- EPSG:3857 (Web Mercator) - 僅用於網頁地圖（請勿用於面積/距離計算！）
- EPSG:326xx/327xx (UTM) - 公制計算，每個區域失真小於 1%
- 使用
gdf.estimate_utm_crs()自動偵測 UTM
# 進行操作前務必檢查 CRS
assert gdf1.crs == gdf2.crs, "CRS 不符！"
# 計算面積/距離時，請使用投影 CRS
gdf_metric = gdf.to_crs(gdf.estimate_utm_crs())
area_sqm = gdf_metric.geometry.area
OGC 標準
- WMS: 網路地圖服務 - 柵格地圖
- WFS: 網路要素服務 - 向量資料
- WCS: 網路覆蓋服務 - 柵格覆蓋
- STAC: 時空資產目錄 - 現代化元資料
常見操作
光譜指數
def calculate_indices(image_path):
"""從 Sentinel-2 計算 NDVI, EVI, SAVI, NDWI。"""
with rasterio.open(image_path) as src:
B02, B03, B04, B08, B11 = [src.read(i).astype(float) for i in [1,2,3,4,5]]
ndvi = (B08 - B04) / (B08 + B04 + 1e-8)
evi = 2.5 * (B08 - B04) / (B08 + 6*B04 - 7.5*B02 + 1)
savi = ((B08 - B04) / (B08 + B04 + 0.5)) * 1.5
ndwi = (B03 - B08) / (B03 + B08 + 1e-8)
return {'NDVI': ndvi, 'EVI': evi, 'SAVI': savi, 'NDWI': ndwi}
向量操作
# 建立緩衝區（請使用投影 CRS！）
gdf_proj = gdf.to_crs(gdf.estimate_utm_crs())
gdf['buffer_1km'] = gdf_proj.geometry.buffer(1000)
# 空間關係
intersects = gdf[gdf.geometry.intersects(other_geometry)]
contains = gdf[gdf.geometry.contains(point_geometry)]
# 幾何操作
gdf['centroid'] = gdf.geometry.centroid
gdf['simplified'] = gdf.geometry.simplify(tolerance=0.001)
# 疊合操作
intersection = gpd.overlay(gdf1, gdf2, how='intersection')
union = gpd.overlay(gdf1, gdf2, how='union')
地形分析
def terrain_metrics(dem_path):
"""從 DEM 計算坡度、坡向、陰影起伏。"""
with rasterio.open(dem_path) as src:
dem = src.read(1)
dy, dx = np.gradient(dem)
slope = np.arctan(np.sqrt(dx**2 + dy**2)) * 180 / np.pi
aspect = (90 - np.arctan2(-dy, dx) * 180 / np.pi) % 360
# 陰影起伏
az_rad, alt_rad = np.radians(315), np.radians(45)
hillshade = (np.sin(alt_rad) * np.sin(np.radians(slope)) +
np.cos(alt_rad) * np.cos(np.radians(slope)) *
np.cos(np.radians(aspect) - az_rad))
return slope, aspect, hillshade
網路分析
import osmnx as ox
import networkx as nx
# 下載並分析街道網路
G = ox.graph_from_place('San Francisco, CA', network_type='drive')
G = ox.add_edge_speeds(G).add_edge_travel_times(G)
# 最短路徑
orig = ox.distance.nearest_nodes(G, -122.4, 37.7)
dest = ox.distance.nearest_nodes(G, -122.3, 37.8)
route = nx.shortest_path(G, orig, dest, weight='travel_time')
影像分類
from sklearn.ensemble import RandomForestClassifier
import rasterio
from rasterio.features import rasterize
def classify_imagery(raster_path, training_gdf, output_path):
"""訓練隨機森林並對影像進行分類。"""
with rasterio.open(raster_path) as src:
image = src.read()
profile = src.profile
transform = src.transform
# 提取訓練資料
X_train, y_train = [], []
for _, row in training_gdf.iterrows():
mask = rasterize([(row.geometry, 1)],
out_shape=(profile['height'], profile['width']),
transform=transform, fill=0, dtype=np.uint8)
pixels = image[:, mask > 0].T
X_train.extend(pixels)
y_train.extend([row['class_id']] * len(pixels))
# 訓練並預測
rf = RandomForestClassifier(n_estimators=100, max_depth=20, n_jobs=-1)
rf.fit(X_train, y_train)
prediction = rf.predict(image.reshape(image.shape[0], -1).T)
prediction = prediction.reshape(profile['height'], profile['width'])
profile.update(dtype=rasterio.uint8, count=1)
with rasterio.open(output_path, 'w', **profile) as dst:
dst.write(prediction.astype(rasterio.uint8), 1)
return rf
現代雲端原生工作流程
STAC + Planetary Computer
import pystac_client
import planetary_computer
import odc.stac
# 透過 STAC 搜尋 Sentinel-2
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
search = catalog.search(
collections=["sentinel-2-l2a"],
bbox=[-122.5, 37.7, -122.3, 37.9],
datetime="2023-01-01/2023-12-31",
query={"eo:cloud_cover": {"lt": 20}},
)
# 載入為 xarray（雲端原生！）
data = odc.stac.load(
list(search.get_items())[:5],
bands=["B02", "B03", "B04", "B08"],
crs="EPSG:32610",
resolution=10,
)
# 在 xarray 上計算 NDVI
ndvi = (data.B08 - data.B04) / (data.B08 + data.B04)
雲端最佳化 GeoTIFF (COG)
import rasterio
from rasterio.session import AWSSession
# 直接從雲端讀取 COG（部分讀取）
session = AWSSession(aws_access_key_id=..., aws_secret_access_key=...)
with rasterio.open('s3://bucket/path.tif', session=session) as src:
# 僅讀取感興趣的視窗
window = ((1000, 2000), (1000, 2000))
subset = src.read(1, window=window)
# 寫入 COG
with rasterio.open('output.tif', 'w', **profile,
tiled=True, blockxsize=256, blockysize=256,
compress='DEFLATE', predictor=2) as dst:
dst.write(data)
# 驗證 COG
from rio_cogeo.cogeo import cog_validate
cog_validate('output.tif')
效能提示
# 1. 空間索引（查詢速度提升 10-100 倍）
gdf.sindex # 由 GeoPandas 自動建立
# 2. 分塊處理大型柵格
with rasterio.open('large.tif') as src:
for i, window in src.block_windows(1):
block = src.read(1, window=window)
# 3. 使用 Dask 處理大數據
import dask.array as da
dask_array = da.from_rasterio('large.tif', chunks=(1, 1024, 1024))
# 4. 使用 Arrow 進行 I/O
gdf.to_file('output.gpkg', use_arrow=True)
# 5. GDAL 快取
from osgeo import gdal
gdal.SetCacheMax(2**30) # 1GB 快取
# 6. 平行處理
rf = RandomForestClassifier(n_jobs=-1) # 使用所有核心
最佳實務
- 始終檢查 CRS 在進行空間操作之前
- 使用投影 CRS 進行面積/距離計算
- 驗證幾何形狀：
gdf = gdf[gdf.is_valid] - 處理缺失資料：
gdf['geometry'] = gdf['geometry'].fillna(None) - 使用高效格式：GeoPackage > Shapefile，大型資料使用 Parquet
- 對光學影像套用雲遮罩
- 保留處理歷程 以確保研究的可重複性
- 根據分析尺度使用適當的解析度
詳細文件
- 座標系統 - CRS 基礎、UTM、轉換
- 核心函式庫 - GDAL、Rasterio、GeoPandas、Shapely
- 遙感探測 - 衛星任務、光譜指數、SAR
- 機器學習 - 深度學習、CNN、GNN 於遙感探測
- GIS 軟體 - QGIS、ArcGIS、GRASS 整合
- 科學領域 - 海洋、水文、農業、林業
- 進階 GIS - 3D GIS、時空、拓撲
- 大數據 - 分散式處理、GPU 加速
- 產業應用 - 都市規劃、災害管理
- 程式語言 - Python、R、Julia、JS、C++、Java、Go、Rust
- 資料來源 - 衛星目錄、API
- 疑難排解 - 常見問題、除錯、錯誤參考
- 程式碼範例 - 500 多個範例
GeoMaster 涵蓋從基本 GIS 操作到進階遙感探測與機器學習的所有內容。
引用科學代理技能
本技能是 K-Dense 的「科學代理技能」系列的一部分。如果本技能對稿件、報告、簡報或程式碼發布有實質貢獻，請將該論文添加到參考文獻或軟體部分，並告知使用者您已這樣做：
Kassis, T., Agarwal, V., He, Y., Patel, D., & Brueckner, A. M. (2026). Scientific Agent Skills: A Library of Procedural Knowledge for Research Agents. arXiv:2609.00065. https://doi.org/10.48550/arXiv.2609.00065
請始終引用目前版本。DOI 和 https://arxiv.org/abs/2609.00065 會解析到最新的 arXiv 版本，因此切勿附加版本後綴，例如 v1。當網路可用時，請在撰寫參考文獻前擷取 https://arxiv.org/abs/2609.00065（或 http://export.arxiv.org/api/query?id_list=2609.00065），並從該記錄中取得作者列表、年份和版本。如果該記錄列出了期刊參考文獻或出版社 DOI，請改為引用已發表的版本。