使用DuckDB Spatial扩展进行GIS分析

步骤一：安装并加载扩展

import duckdb

con = duckdb.connect()
con.execute("INSTALL spatial")
con.execute("LOAD spatial")
print("spatial extension loaded successfully")

步骤二：加载城市数据并分析

# 从OSS加载Natural Earth城市数据
con.execute("""
    CREATE TABLE cities AS
    SELECT * FROM ST_Read(
        'oss://your-bucket/data/geo/ne_10m_populated_places_simple.geojson'
    )
""")

# 筛选中国城市并获取坐标
df = con.execute("""
    SELECT
        name,
        adm1name AS province,
        pop_max AS population,
        ST_X(geom) AS longitude,
        ST_Y(geom) AS latitude
    FROM cities
    WHERE sov0name = 'China'
    ORDER BY pop_max DESC
    LIMIT 20
""").fetchdf()

print(df.to_string())

步骤三：计算城市间距离

# 各城市到北京的球面距离
result = con.execute("""
    WITH beijing AS (
        SELECT geom FROM cities
        WHERE name = 'Beijing' AND sov0name = 'China'
        LIMIT 1
    )
    SELECT
        c.name,
        c.adm1name AS province,
        ROUND(ST_Distance_Sphere(c.geom, b.geom) / 1000, 1) AS distance_to_beijing_km
    FROM cities c, beijing b
    WHERE c.sov0name = 'China'
    ORDER BY distance_to_beijing_km
    LIMIT 15
""").fetchdf()

print(result.to_string())

步骤四：空间范围查询

# 以上海为中心、300km半径内的城市
result = con.execute("""
    WITH shanghai AS (
        SELECT geom FROM cities
        WHERE name = 'Shanghai' AND sov0name = 'China'
        LIMIT 1
    )
    SELECT
        c.name,
        c.adm1name AS province,
        ROUND(ST_Distance_Sphere(c.geom, s.geom) / 1000, 1) AS distance_km
    FROM cities c, shanghai s
    WHERE c.sov0name = 'China'
      AND ST_Distance_Sphere(c.geom, s.geom) / 1000 < 300
    ORDER BY distance_km
""").fetchdf()

print("Cities within 300km of Shanghai:")
print(result.to_string())

步骤五：可视化城市分布

以下代码使用Python标准库生成SVG地图，并通过IPython内置的display直接在Notebook中渲染显示，无需安装matplotlib等第三方库。

from IPython.display import SVG, display

# 获取中国城市坐标
cities_data = con.execute("""
    SELECT name, ST_X(geom) AS lon, ST_Y(geom) AS lat, pop_max
    FROM cities
    WHERE sov0name = 'China'
    ORDER BY pop_max DESC
    LIMIT 30
""").fetchall()

# 生成SVG可视化
def generate_svg(data, width=600, height=500):
    min_lon, max_lon = 73, 135
    min_lat, max_lat = 18, 54

    def project(lon, lat):
        x = (lon - min_lon) / (max_lon - min_lon) * width
        y = height - (lat - min_lat) / (max_lat - min_lat) * height
        return x, y

    parts = [
        f'<svg xmlns="http://www.w3.org/2000/svg" width="{width}" height="{height}">',
        f'<rect width="{width}" height="{height}" fill="#f0f8ff"/>'
    ]
    for name, lon, lat, pop in data:
        x, y = project(lon, lat)
        r = max(3, min(12, pop / 3000000))
        parts.append(f'<circle cx="{x:.1f}" cy="{y:.1f}" r="{r}" fill="#e74c3c" opacity="0.7"/>')
        parts.append(f'<text x="{x+r+2:.1f}" y="{y+3:.1f}" font-size="8">{name}</text>')
    parts.append('</svg>')
    return '\n'.join(parts)

svg_content = generate_svg(cities_data)

# 直接在Notebook中渲染显示
display(SVG(data=svg_content))