第18章-高级应用与案例
第18章：高级应用与案例
18.1 概述
本章通过多个实际案例，展示 PostGIS 在不同领域的高级应用。这些案例涵盖了智慧城市、交通物流、环境监测、商业分析等多个领域，帮助读者将所学知识应用到实际项目中。
18.2 智慧城市应用案例
18.2.1 城市设施选址分析
需求描述：为某城市规划新建一个消防站，要求：
- 距离现有消防站至少 2 公里
- 距离主要道路不超过 500 米
- 能够覆盖尽可能多的居民区
数据表设计
-- 现有消防站
CREATE TABLE fire_stations (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
address VARCHAR(200),
capacity INTEGER,
geom GEOMETRY(POINT, 4326)
);
-- 道路网络
CREATE TABLE roads (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
road_type VARCHAR(50), -- 主干道、次干道、支路
lanes INTEGER,
geom GEOMETRY(LINESTRING, 4326)
);
-- 居民区
CREATE TABLE residential_areas (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
population INTEGER,
households INTEGER,
geom GEOMETRY(POLYGON, 4326)
);
-- 城市边界
CREATE TABLE city_boundary (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
geom GEOMETRY(POLYGON, 4326)
);
-- 创建空间索引
CREATE INDEX idx_fire_stations_geom ON fire_stations USING GIST (geom);
CREATE INDEX idx_roads_geom ON roads USING GIST (geom);
CREATE INDEX idx_residential_geom ON residential_areas USING GIST (geom);
CREATE INDEX idx_city_boundary_geom ON city_boundary USING GIST (geom);
选址分析实现
-- 步骤1：找出距离现有消防站超过2公里的区域
WITH existing_coverage AS (
SELECT ST_Union(ST_Buffer(geom::geography, 2000)::geometry) AS buffer
FROM fire_stations
),
-- 步骤2：获取主要道路500米缓冲区
main_road_buffer AS (
SELECT ST_Union(ST_Buffer(geom::geography, 500)::geometry) AS buffer
FROM roads
WHERE road_type IN ('主干道', '次干道')
),
-- 步骤3：计算可选区域
available_area AS (
SELECT ST_Intersection(
ST_Difference(
(SELECT geom FROM city_boundary WHERE name = '主城区'),
(SELECT buffer FROM existing_coverage)
),
(SELECT buffer FROM main_road_buffer)
) AS geom
),
-- 步骤4：生成候选点（网格采样）
candidate_points AS (
SELECT (ST_Dump(ST_GeneratePoints(geom, 100))).geom AS geom
FROM available_area
),
-- 步骤5：评估每个候选点的覆盖能力
scored_candidates AS (
SELECT
cp.geom,
(
SELECT COALESCE(SUM(ra.population), 0)
FROM residential_areas ra
WHERE ST_DWithin(ra.geom::geography, cp.geom::geography, 3000)
) AS covered_population,
(
SELECT COUNT(*)
FROM residential_areas ra
WHERE ST_DWithin(ra.geom::geography, cp.geom::geography, 3000)
) AS covered_areas
FROM candidate_points cp
)
-- 输出最佳候选位置
SELECT
ST_X(geom) AS longitude,
ST_Y(geom) AS latitude,
covered_population,
covered_areas,
ST_AsGeoJSON(geom) AS geojson
FROM scored_candidates
ORDER BY covered_population DESC
LIMIT 10;
服务范围可视化
-- 创建视图：各消防站的服务范围
CREATE OR REPLACE VIEW fire_station_service_areas AS
SELECT
fs.id,
fs.name,
ST_Buffer(fs.geom::geography, 3000)::geometry AS service_area_3km,
ST_Buffer(fs.geom::geography, 5000)::geometry AS service_area_5km,
(
SELECT COALESCE(SUM(population), 0)
FROM residential_areas
WHERE ST_DWithin(geom::geography, fs.geom::geography, 3000)
) AS pop_within_3km,
(
SELECT COALESCE(SUM(population), 0)
FROM residential_areas
WHERE ST_DWithin(geom::geography, fs.geom::geography, 5000)
) AS pop_within_5km
FROM fire_stations fs;
-- 分析服务覆盖缺口
WITH total_coverage AS (
SELECT ST_Union(ST_Buffer(geom::geography, 5000)::geometry) AS covered
FROM fire_stations
)
SELECT
ST_Area(ST_Difference(
(SELECT geom FROM city_boundary WHERE name = '主城区'),
(SELECT covered FROM total_coverage)
)::geography) / 1000000.0 AS uncovered_area_km2;
18.2.2 城市热力图分析
POI 密度分析
-- 创建 POI 表
CREATE TABLE poi (
id SERIAL PRIMARY KEY,
name VARCHAR(200),
category VARCHAR(50),
sub_category VARCHAR(50),
address VARCHAR(300),
geom GEOMETRY(POINT, 4326)
);
-- 创建六边形网格
CREATE OR REPLACE FUNCTION generate_hexgrid(
bounds geometry,
cell_size float,
srid integer DEFAULT 4326
)
RETURNS TABLE(geom geometry) AS $$
DECLARE
x_min float;
y_min float;
x_max float;
y_max float;
hex_width float;
hex_height float;
BEGIN
x_min := ST_XMin(bounds);
y_min := ST_YMin(bounds);
x_max := ST_XMax(bounds);
y_max := ST_YMax(bounds);
hex_width := cell_size * 2;
hex_height := cell_size * sqrt(3);
RETURN QUERY
WITH rows AS (
SELECT generate_series(0, ceil((y_max - y_min) / hex_height)::int) AS row_num
),
cols AS (
SELECT generate_series(0, ceil((x_max - x_min) / (hex_width * 0.75))::int) AS col_num
),
hex_centers AS (
SELECT
x_min + col_num * hex_width * 0.75 AS cx,
y_min + row_num * hex_height +
CASE WHEN col_num % 2 = 1 THEN hex_height / 2 ELSE 0 END AS cy
FROM rows CROSS JOIN cols
)
SELECT ST_SetSRID(
ST_MakePolygon(
ST_MakeLine(ARRAY[
ST_MakePoint(cx + cell_size, cy),
ST_MakePoint(cx + cell_size/2, cy + cell_size*sqrt(3)/2),
ST_MakePoint(cx - cell_size/2, cy + cell_size*sqrt(3)/2),
ST_MakePoint(cx - cell_size, cy),
ST_MakePoint(cx - cell_size/2, cy - cell_size*sqrt(3)/2),
ST_MakePoint(cx + cell_size/2, cy - cell_size*sqrt(3)/2),
ST_MakePoint(cx + cell_size, cy)
])
),
srid
)
FROM hex_centers;
END;
$$ LANGUAGE plpgsql;
-- 生成热力图数据
WITH hexgrid AS (
SELECT geom FROM generate_hexgrid(
(SELECT ST_Envelope(ST_Collect(geom)) FROM poi),
0.005 -- 约500米
)
)
SELECT
hg.geom,
COUNT(p.id) AS poi_count,
COUNT(p.id) FILTER (WHERE p.category = '餐饮') AS restaurant_count,
COUNT(p.id) FILTER (WHERE p.category = '购物') AS shopping_count,
COUNT(p.id) FILTER (WHERE p.category = '住宿') AS hotel_count
FROM hexgrid hg
LEFT JOIN poi p ON ST_Contains(hg.geom, p.geom)
GROUP BY hg.geom
HAVING COUNT(p.id) > 0;
18.3 交通物流应用案例
18.3.1 配送路线优化
数据表设计
-- 配送中心
CREATE TABLE distribution_centers (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
capacity INTEGER,
vehicles INTEGER,
geom GEOMETRY(POINT, 4326)
);
-- 客户点
CREATE TABLE customers (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
address VARCHAR(200),
demand INTEGER, -- 需求量
time_window_start TIME,
time_window_end TIME,
priority INTEGER DEFAULT 1,
geom GEOMETRY(POINT, 4326)
);
-- 道路网络（用于路径规划）
CREATE TABLE road_network (
id SERIAL PRIMARY KEY,
source INTEGER,
target INTEGER,
length FLOAT,
speed_limit INTEGER,
road_class VARCHAR(20),
geom GEOMETRY(LINESTRING, 4326)
);
-- 配送任务
CREATE TABLE delivery_tasks (
id SERIAL PRIMARY KEY,
center_id INTEGER REFERENCES distribution_centers(id),
customer_id INTEGER REFERENCES customers(id),
status VARCHAR(20) DEFAULT 'pending',
planned_time TIMESTAMP,
actual_time TIMESTAMP,
route_geom GEOMETRY(LINESTRING, 4326)
);
-- 创建索引
CREATE INDEX idx_dc_geom ON distribution_centers USING GIST (geom);
CREATE INDEX idx_customers_geom ON customers USING GIST (geom);
CREATE INDEX idx_road_network_geom ON road_network USING GIST (geom);
配送区域划分
-- 使用 Voronoi 图划分配送区域
WITH dc_points AS (
SELECT id, name, geom FROM distribution_centers
),
voronoi AS (
SELECT
(ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom AS poly
FROM dc_points
)
SELECT
dc.id AS center_id,
dc.name AS center_name,
ST_Intersection(v.poly, cb.geom) AS service_area
FROM voronoi v
CROSS JOIN city_boundary cb
JOIN dc_points dc ON ST_Contains(v.poly, dc.geom)
WHERE ST_IsValid(cb.geom);
-- 分配客户到最近的配送中心
UPDATE customers c
SET assigned_center = (
SELECT dc.id
FROM distribution_centers dc
ORDER BY c.geom <-> dc.geom
LIMIT 1
);
-- 计算每个配送中心的负载
SELECT
dc.id,
dc.name,
COUNT(c.id) AS customer_count,
SUM(c.demand) AS total_demand,
AVG(ST_Distance(c.geom::geography, dc.geom::geography)) AS avg_distance_m
FROM distribution_centers dc
LEFT JOIN customers c ON c.assigned_center = dc.id
GROUP BY dc.id, dc.name
ORDER BY customer_count DESC;
路径规划（需要 pgRouting 扩展）
-- 安装 pgRouting
CREATE EXTENSION IF NOT EXISTS pgrouting;
-- 准备拓扑
SELECT pgr_createTopology('road_network', 0.0001, 'geom', 'id');
-- 最短路径查询
SELECT
seq,
node,
edge,
cost,
agg_cost,
geom
FROM pgr_dijkstra(
'SELECT id, source, target, length AS cost FROM road_network',
(SELECT id FROM road_network_vertices_pgr ORDER BY geom <->
(SELECT geom FROM distribution_centers WHERE name = '配送中心A') LIMIT 1),
(SELECT id FROM road_network_vertices_pgr ORDER BY geom <->
(SELECT geom FROM customers WHERE name = '客户X') LIMIT 1),
directed := false
) AS route
JOIN road_network rn ON route.edge = rn.id;
-- 多点路径（旅行商问题近似解）
WITH customer_nodes AS (
SELECT
c.id AS customer_id,
(SELECT v.id FROM road_network_vertices_pgr v
ORDER BY v.geom <-> c.geom LIMIT 1) AS node_id
FROM customers c
WHERE c.assigned_center = 1
LIMIT 10
)
SELECT * FROM pgr_dijkstraTSP(
'SELECT id, source, target, length AS cost FROM road_network',
(SELECT array_agg(node_id) FROM customer_nodes),
1 -- 起点
);
18.3.2 车辆轨迹分析
轨迹数据存储
-- 车辆轨迹表
CREATE TABLE vehicle_tracks (
id SERIAL PRIMARY KEY,
vehicle_id VARCHAR(20),
timestamp TIMESTAMP WITH TIME ZONE,
speed FLOAT,
heading FLOAT,
geom GEOMETRY(POINT, 4326)
);
-- 分区表（按日期分区提高性能）
CREATE TABLE vehicle_tracks_partitioned (
id SERIAL,
vehicle_id VARCHAR(20),
timestamp TIMESTAMP WITH TIME ZONE,
speed FLOAT,
heading FLOAT,
geom GEOMETRY(POINT, 4326)
) PARTITION BY RANGE (timestamp);
-- 创建分区
CREATE TABLE vehicle_tracks_2024_01 PARTITION OF vehicle_tracks_partitioned
FOR VALUES FROM ('2024-01-01') TO ('2024-02-01');
CREATE TABLE vehicle_tracks_2024_02 PARTITION OF vehicle_tracks_partitioned
FOR VALUES FROM ('2024-02-01') TO ('2024-03-01');
-- 创建索引
CREATE INDEX idx_tracks_vehicle_time ON vehicle_tracks (vehicle_id, timestamp);
CREATE INDEX idx_tracks_geom ON vehicle_tracks USING GIST (geom);
CREATE INDEX idx_tracks_time ON vehicle_tracks (timestamp);
轨迹分析函数
-- 生成轨迹线
CREATE OR REPLACE FUNCTION generate_trajectory(
p_vehicle_id VARCHAR,
p_start_time TIMESTAMP,
p_end_time TIMESTAMP
)
RETURNS GEOMETRY AS $$
SELECT ST_MakeLine(geom ORDER BY timestamp)
FROM vehicle_tracks
WHERE vehicle_id = p_vehicle_id
AND timestamp BETWEEN p_start_time AND p_end_time;
$$ LANGUAGE SQL;
-- 计算行驶里程
CREATE OR REPLACE FUNCTION calculate_mileage(
p_vehicle_id VARCHAR,
p_start_time TIMESTAMP,
p_end_time TIMESTAMP
)
RETURNS FLOAT AS $$
DECLARE
trajectory GEOMETRY;
BEGIN
trajectory := generate_trajectory(p_vehicle_id, p_start_time, p_end_time);
IF trajectory IS NULL THEN
RETURN 0;
END IF;
RETURN ST_Length(trajectory::geography);
END;
$$ LANGUAGE plpgsql;
-- 停留点检测
CREATE OR REPLACE FUNCTION detect_stops(
p_vehicle_id VARCHAR,
p_start_time TIMESTAMP,
p_end_time TIMESTAMP,
p_min_duration INTERVAL DEFAULT '5 minutes',
p_max_radius FLOAT DEFAULT 50 -- 米
)
RETURNS TABLE (
stop_start TIMESTAMP,
stop_end TIMESTAMP,
duration INTERVAL,
center_point GEOMETRY,
address TEXT
) AS $$
WITH ordered_points AS (
SELECT
geom,
timestamp,
LAG(geom) OVER w AS prev_geom,
LAG(timestamp) OVER w AS prev_timestamp
FROM vehicle_tracks
WHERE vehicle_id = p_vehicle_id
AND timestamp BETWEEN p_start_time AND p_end_time
WINDOW w AS (ORDER BY timestamp)
),
stop_segments AS (
SELECT
timestamp,
geom,
CASE
WHEN ST_Distance(geom::geography, prev_geom::geography) < p_max_radius
THEN 0 ELSE 1
END AS is_moving
FROM ordered_points
WHERE prev_geom IS NOT NULL
),
stop_groups AS (
SELECT
timestamp,
geom,
SUM(is_moving) OVER (ORDER BY timestamp) AS group_id
FROM stop_segments
WHERE is_moving = 0
)
SELECT
MIN(timestamp) AS stop_start,
MAX(timestamp) AS stop_end,
MAX(timestamp) - MIN(timestamp) AS duration,
ST_Centroid(ST_Collect(geom)) AS center_point,
NULL::TEXT AS address
FROM stop_groups
GROUP BY group_id
HAVING MAX(timestamp) - MIN(timestamp) >= p_min_duration;
$$ LANGUAGE SQL;
-- 使用示例
SELECT * FROM detect_stops('京A12345', '2024-01-15 08:00:00', '2024-01-15 18:00:00');
围栏检测
-- 电子围栏表
CREATE TABLE geofences (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
fence_type VARCHAR(20), -- 'enter', 'exit', 'both'
geom GEOMETRY(POLYGON, 4326)
);
-- 围栏事件记录
CREATE TABLE geofence_events (
id SERIAL PRIMARY KEY,
vehicle_id VARCHAR(20),
fence_id INTEGER REFERENCES geofences(id),
event_type VARCHAR(10), -- 'enter', 'exit'
event_time TIMESTAMP,
geom GEOMETRY(POINT, 4326)
);
-- 围栏检测触发器
CREATE OR REPLACE FUNCTION check_geofence()
RETURNS TRIGGER AS $$
DECLARE
prev_point GEOMETRY;
fence RECORD;
BEGIN
-- 获取上一个位置点
SELECT geom INTO prev_point
FROM vehicle_tracks
WHERE vehicle_id = NEW.vehicle_id
AND timestamp < NEW.timestamp
ORDER BY timestamp DESC
LIMIT 1;
IF prev_point IS NULL THEN
RETURN NEW;
END IF;
-- 检查所有围栏
FOR fence IN SELECT * FROM geofences LOOP
-- 检测进入
IF NOT ST_Contains(fence.geom, prev_point)
AND ST_Contains(fence.geom, NEW.geom) THEN
INSERT INTO geofence_events (vehicle_id, fence_id, event_type, event_time, geom)
VALUES (NEW.vehicle_id, fence.id, 'enter', NEW.timestamp, NEW.geom);
END IF;
-- 检测离开
IF ST_Contains(fence.geom, prev_point)
AND NOT ST_Contains(fence.geom, NEW.geom) THEN
INSERT INTO geofence_events (vehicle_id, fence_id, event_type, event_time, geom)
VALUES (NEW.vehicle_id, fence.id, 'exit', NEW.timestamp, NEW.geom);
END IF;
END LOOP;
RETURN NEW;
END;
$$ LANGUAGE plpgsql;
CREATE TRIGGER trigger_check_geofence
AFTER INSERT ON vehicle_tracks
FOR EACH ROW EXECUTE FUNCTION check_geofence();
18.4 环境监测应用案例
18.4.1 空气质量监测
监测站数据管理
-- 监测站表
CREATE TABLE air_monitoring_stations (
id SERIAL PRIMARY KEY,
station_code VARCHAR(20) UNIQUE,
station_name VARCHAR(100),
station_type VARCHAR(50), -- 国控、省控、市控
status VARCHAR(20) DEFAULT 'active',
geom GEOMETRY(POINT, 4326)
);
-- 监测数据表
CREATE TABLE air_quality_data (
id SERIAL PRIMARY KEY,
station_id INTEGER REFERENCES air_monitoring_stations(id),
measure_time TIMESTAMP,
pm25 FLOAT,
pm10 FLOAT,
so2 FLOAT,
no2 FLOAT,
co FLOAT,
o3 FLOAT,
aqi INTEGER,
primary_pollutant VARCHAR(20)
);
-- 创建索引
CREATE INDEX idx_aqs_geom ON air_monitoring_stations USING GIST (geom);
CREATE INDEX idx_aqd_station_time ON air_quality_data (station_id, measure_time);
CREATE INDEX idx_aqd_time ON air_quality_data (measure_time);
空气质量插值分析
-- 反距离加权插值 (IDW)
CREATE OR REPLACE FUNCTION idw_interpolation(
target_point GEOMETRY,
measure_time TIMESTAMP,
pollutant VARCHAR,
power FLOAT DEFAULT 2
)
RETURNS FLOAT AS $$
DECLARE
result FLOAT;
BEGIN
EXECUTE format('
WITH station_values AS (
SELECT
s.geom,
d.%I AS value,
ST_Distance(s.geom::geography, $1::geography) AS distance
FROM air_monitoring_stations s
JOIN air_quality_data d ON s.id = d.station_id
WHERE d.measure_time = $2
AND d.%I IS NOT NULL
),
weighted AS (
SELECT
value,
CASE
WHEN distance < 1 THEN 1e10
ELSE 1.0 / power(distance, $3)
END AS weight
FROM station_values
)
SELECT SUM(value * weight) / SUM(weight)
FROM weighted
', pollutant, pollutant)
INTO result
USING target_point, measure_time, power;
RETURN result;
END;
$$ LANGUAGE plpgsql;
-- 生成网格插值图
WITH grid AS (
SELECT (ST_PixelAsPolygons(
ST_AsRaster(
ST_MakeEnvelope(116.0, 39.5, 117.0, 40.5, 4326),
100, 100
)
)).*
),
interpolated AS (
SELECT
g.geom,
idw_interpolation(ST_Centroid(g.geom), '2024-01-15 12:00:00', 'pm25') AS pm25_value
FROM grid g
)
SELECT
geom,
pm25_value,
CASE
WHEN pm25_value <= 35 THEN '优'
WHEN pm25_value <= 75 THEN '良'
WHEN pm25_value <= 115 THEN '轻度污染'
WHEN pm25_value <= 150 THEN '中度污染'
WHEN pm25_value <= 250 THEN '重度污染'
ELSE '严重污染'
END AS quality_level
FROM interpolated
WHERE pm25_value IS NOT NULL;
18.4.2 水质监测与流域分析
流域数据管理
-- 河流表
CREATE TABLE rivers (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
river_level INTEGER, -- 1-5级河流
length_km FLOAT,
basin_name VARCHAR(100),
geom GEOMETRY(LINESTRING, 4326)
);
-- 监测断面表
CREATE TABLE water_monitoring_sections (
id SERIAL PRIMARY KEY,
section_code VARCHAR(20),
section_name VARCHAR(100),
river_id INTEGER REFERENCES rivers(id),
section_type VARCHAR(50),
geom GEOMETRY(POINT, 4326)
);
-- 水质数据表
CREATE TABLE water_quality_data (
id SERIAL PRIMARY KEY,
section_id INTEGER REFERENCES water_monitoring_sections(id),
sample_time TIMESTAMP,
ph FLOAT,
dissolved_oxygen FLOAT,
cod FLOAT,
bod5 FLOAT,
ammonia_nitrogen FLOAT,
total_phosphorus FLOAT,
water_grade VARCHAR(10)
);
-- 汇水区域表
CREATE TABLE catchments (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
area_km2 FLOAT,
main_river_id INTEGER REFERENCES rivers(id),
geom GEOMETRY(POLYGON, 4326)
);
流域分析查询
-- 上游监测断面查询
CREATE OR REPLACE FUNCTION get_upstream_sections(
p_section_id INTEGER
)
RETURNS TABLE (
section_id INTEGER,
section_name VARCHAR,
distance_km FLOAT
) AS $$
WITH RECURSIVE section_info AS (
SELECT
wms.id,
wms.section_name,
wms.river_id,
wms.geom,
ST_LineLocatePoint(r.geom, wms.geom) AS position
FROM water_monitoring_sections wms
JOIN rivers r ON wms.river_id = r.id
WHERE wms.id = p_section_id
),
upstream AS (
SELECT
wms.id,
wms.section_name,
wms.river_id,
wms.geom,
ST_LineLocatePoint(r.geom, wms.geom) AS position
FROM water_monitoring_sections wms
JOIN rivers r ON wms.river_id = r.id
JOIN section_info si ON wms.river_id = si.river_id
WHERE ST_LineLocatePoint(r.geom, wms.geom) < si.position
)
SELECT
u.id AS section_id,
u.section_name,
ST_Distance(u.geom::geography, si.geom::geography) / 1000.0 AS distance_km
FROM upstream u
CROSS JOIN section_info si
ORDER BY distance_km;
$$ LANGUAGE SQL;
-- 污染源影响分析
WITH pollution_sources AS (
SELECT id, name, discharge_amount, geom
FROM industrial_enterprises
WHERE industry_type = '化工'
),
affected_sections AS (
SELECT
ps.name AS source_name,
wms.section_name,
ST_Distance(ps.geom::geography, wms.geom::geography) AS distance_m
FROM pollution_sources ps
CROSS JOIN water_monitoring_sections wms
WHERE ST_DWithin(ps.geom::geography, wms.geom::geography, 10000)
)
SELECT * FROM affected_sections ORDER BY distance_m;
18.5 商业分析应用案例
18.5.1 商业选址分析
数据准备
-- 商圈表
CREATE TABLE business_districts (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
district_level VARCHAR(20), -- 市级、区级、社区级
avg_rent_price FLOAT,
geom GEOMETRY(POLYGON, 4326)
);
-- 竞争对手门店
CREATE TABLE competitor_stores (
id SERIAL PRIMARY KEY,
brand VARCHAR(100),
store_name VARCHAR(200),
category VARCHAR(50),
estimated_revenue FLOAT,
geom GEOMETRY(POINT, 4326)
);
-- 人口数据（网格化）
CREATE TABLE population_grid (
id SERIAL PRIMARY KEY,
total_population INTEGER,
working_population INTEGER,
avg_income FLOAT,
consumption_level VARCHAR(20),
geom GEOMETRY(POLYGON, 4326)
);
-- 交通设施
CREATE TABLE transport_facilities (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
facility_type VARCHAR(50), -- 地铁站、公交站、停车场
daily_flow INTEGER,
geom GEOMETRY(POINT, 4326)
);
选址评分模型
-- 创建选址评分函数
CREATE OR REPLACE FUNCTION calculate_site_score(
site_point GEOMETRY,
radius_m FLOAT DEFAULT 1000
)
RETURNS TABLE (
dimension VARCHAR,
score FLOAT,
weight FLOAT,
weighted_score FLOAT,
details JSONB
) AS $$
BEGIN
-- 人口维度评分
RETURN QUERY
WITH pop_stats AS (
SELECT
COALESCE(SUM(total_population), 0) AS total_pop,
COALESCE(AVG(avg_income), 0) AS avg_income
FROM population_grid
WHERE ST_DWithin(geom::geography, site_point::geography, radius_m)
)
SELECT
'人口密度'::VARCHAR,
LEAST(total_pop / 10000.0, 10.0)::FLOAT,
0.25::FLOAT,
LEAST(total_pop / 10000.0, 10.0) * 0.25,
jsonb_build_object('population', total_pop, 'avg_income', avg_income)
FROM pop_stats;
-- 交通便利性评分
RETURN QUERY
WITH transport_stats AS (
SELECT
COUNT(*) FILTER (WHERE facility_type = '地铁站') AS subway_count,
COUNT(*) FILTER (WHERE facility_type = '公交站') AS bus_count,
COALESCE(MIN(ST_Distance(geom::geography, site_point::geography)), 9999) AS nearest_subway_m
FROM transport_facilities
WHERE ST_DWithin(geom::geography, site_point::geography, radius_m)
)
SELECT
'交通便利性'::VARCHAR,
(CASE WHEN nearest_subway_m < 500 THEN 10
WHEN nearest_subway_m < 1000 THEN 8
ELSE GREATEST(10 - nearest_subway_m / 200, 0) END)::FLOAT,
0.20::FLOAT,
(CASE WHEN nearest_subway_m < 500 THEN 10
WHEN nearest_subway_m < 1000 THEN 8
ELSE GREATEST(10 - nearest_subway_m / 200, 0) END) * 0.20,
jsonb_build_object('subway_stations', subway_count, 'bus_stops', bus_count,
'nearest_subway_m', nearest_subway_m)
FROM transport_stats;
-- 竞争环境评分
RETURN QUERY
WITH competitor_stats AS (
SELECT COUNT(*) AS competitor_count
FROM competitor_stores
WHERE ST_DWithin(geom::geography, site_point::geography, radius_m)
)
SELECT
'竞争环境'::VARCHAR,
GREATEST(10 - competitor_count * 2, 0)::FLOAT,
0.30::FLOAT,
GREATEST(10 - competitor_count * 2, 0) * 0.30,
jsonb_build_object('competitors', competitor_count)
FROM competitor_stats;
-- 商业氛围评分
RETURN QUERY
SELECT
'商业氛围'::VARCHAR,
CASE
WHEN EXISTS (SELECT 1 FROM business_districts bd
WHERE ST_Contains(bd.geom, site_point)
AND bd.district_level = '市级') THEN 10.0
WHEN EXISTS (SELECT 1 FROM business_districts bd
WHERE ST_Contains(bd.geom, site_point)
AND bd.district_level = '区级') THEN 7.0
WHEN EXISTS (SELECT 1 FROM business_districts bd
WHERE ST_Contains(bd.geom, site_point)) THEN 5.0
ELSE 2.0
END::FLOAT,
0.25::FLOAT,
CASE
WHEN EXISTS (SELECT 1 FROM business_districts bd
WHERE ST_Contains(bd.geom, site_point)
AND bd.district_level = '市级') THEN 10.0 * 0.25
WHEN EXISTS (SELECT 1 FROM business_districts bd
WHERE ST_Contains(bd.geom, site_point)
AND bd.district_level = '区级') THEN 7.0 * 0.25
WHEN EXISTS (SELECT 1 FROM business_districts bd
WHERE ST_Contains(bd.geom, site_point)) THEN 5.0 * 0.25
ELSE 2.0 * 0.25
END,
jsonb_build_object('in_business_district',
EXISTS (SELECT 1 FROM business_districts WHERE ST_Contains(geom, site_point)));
END;
$$ LANGUAGE plpgsql;
-- 批量评估候选点
WITH candidate_sites AS (
SELECT
id,
name,
geom
FROM potential_sites
),
site_scores AS (
SELECT
cs.id,
cs.name,
cs.geom,
(SELECT SUM(weighted_score) FROM calculate_site_score(cs.geom)) AS total_score
FROM candidate_sites cs
)
SELECT
id,
name,
total_score,
ST_AsGeoJSON(geom) AS geojson,
RANK() OVER (ORDER BY total_score DESC) AS ranking
FROM site_scores
ORDER BY total_score DESC;
18.5.2 客户分布分析
客户聚类分析
-- 客户表
CREATE TABLE customers_business (
id SERIAL PRIMARY KEY,
customer_id VARCHAR(50),
name VARCHAR(100),
customer_type VARCHAR(20),
total_purchase FLOAT,
visit_frequency INTEGER,
last_visit_date DATE,
geom GEOMETRY(POINT, 4326)
);
-- 基于 DBSCAN 的空间聚类
-- 需要 PostGIS 3.0+
CREATE OR REPLACE FUNCTION customer_clustering(
eps FLOAT DEFAULT 0.01, -- 约1公里
min_points INTEGER DEFAULT 5
)
RETURNS TABLE (
customer_id VARCHAR,
cluster_id INTEGER,
geom GEOMETRY
) AS $$
SELECT
c.customer_id,
ST_ClusterDBSCAN(c.geom, eps, min_points) OVER () AS cluster_id,
c.geom
FROM customers_business c;
$$ LANGUAGE SQL;
-- 分析聚类结果
WITH clustered AS (
SELECT * FROM customer_clustering(0.01, 5)
),
cluster_stats AS (
SELECT
cluster_id,
COUNT(*) AS customer_count,
ST_Centroid(ST_Collect(geom)) AS center,
ST_ConvexHull(ST_Collect(geom)) AS boundary,
AVG(c.total_purchase) AS avg_purchase,
SUM(c.total_purchase) AS total_revenue
FROM clustered cl
JOIN customers_business c ON cl.customer_id = c.customer_id
WHERE cluster_id IS NOT NULL
GROUP BY cluster_id
)
SELECT
cluster_id,
customer_count,
avg_purchase,
total_revenue,
ST_Area(boundary::geography) / 1000000.0 AS area_km2,
ST_AsGeoJSON(center) AS center_geojson,
ST_AsGeoJSON(boundary) AS boundary_geojson
FROM cluster_stats
ORDER BY total_revenue DESC;
18.6 农业应用案例
18.6.1 精准农业地块管理
农业数据模型
-- 农田地块表
CREATE TABLE farmland_parcels (
id SERIAL PRIMARY KEY,
parcel_code VARCHAR(50) UNIQUE,
owner_name VARCHAR(100),
land_type VARCHAR(50), -- 耕地、园地、林地
soil_type VARCHAR(50),
area_mu FLOAT, -- 亩
geom GEOMETRY(POLYGON, 4326)
);
-- 种植记录表
CREATE TABLE planting_records (
id SERIAL PRIMARY KEY,
parcel_id INTEGER REFERENCES farmland_parcels(id),
crop_type VARCHAR(50),
variety VARCHAR(100),
planting_date DATE,
harvest_date DATE,
yield_kg FLOAT,
season VARCHAR(20)
);
-- 土壤检测数据
CREATE TABLE soil_samples (
id SERIAL PRIMARY KEY,
parcel_id INTEGER REFERENCES farmland_parcels(id),
sample_date DATE,
ph FLOAT,
organic_matter FLOAT,
nitrogen FLOAT,
phosphorus FLOAT,
potassium FLOAT,
sample_point GEOMETRY(POINT, 4326)
);
-- 遥感影像分析结果
CREATE TABLE ndvi_data (
id SERIAL PRIMARY KEY,
parcel_id INTEGER REFERENCES farmland_parcels(id),
image_date DATE,
ndvi_mean FLOAT,
ndvi_std FLOAT,
health_status VARCHAR(20)
);
农业分析功能
-- 计算地块产量对比
WITH parcel_yields AS (
SELECT
fp.id,
fp.parcel_code,
fp.area_mu,
pr.crop_type,
pr.yield_kg,
pr.yield_kg / fp.area_mu AS yield_per_mu,
EXTRACT(YEAR FROM pr.harvest_date) AS year
FROM farmland_parcels fp
JOIN planting_records pr ON fp.id = pr.parcel_id
WHERE pr.harvest_date IS NOT NULL
),
avg_yields AS (
SELECT
crop_type,
year,
AVG(yield_per_mu) AS avg_yield_per_mu
FROM parcel_yields
GROUP BY crop_type, year
)
SELECT
py.parcel_code,
py.crop_type,
py.year,
py.yield_per_mu,
ay.avg_yield_per_mu,
(py.yield_per_mu - ay.avg_yield_per_mu) / ay.avg_yield_per_mu * 100 AS deviation_pct,
CASE
WHEN py.yield_per_mu > ay.avg_yield_per_mu * 1.2 THEN '高产'
WHEN py.yield_per_mu < ay.avg_yield_per_mu * 0.8 THEN '低产'
ELSE '正常'
END AS yield_level
FROM parcel_yields py
JOIN avg_yields ay ON py.crop_type = ay.crop_type AND py.year = ay.year
ORDER BY deviation_pct DESC;
-- 土壤质量空间分布
SELECT
fp.parcel_code,
fp.geom,
ss.ph,
ss.organic_matter,
CASE
WHEN ss.ph BETWEEN 6.5 AND 7.5 AND ss.organic_matter > 20 THEN '优良'
WHEN ss.ph BETWEEN 6.0 AND 8.0 AND ss.organic_matter > 15 THEN '良好'
ELSE '一般'
END AS soil_quality
FROM farmland_parcels fp
JOIN LATERAL (
SELECT * FROM soil_samples
WHERE parcel_id = fp.id
ORDER BY sample_date DESC
LIMIT 1
) ss ON true;
-- 农作物健康监测（基于NDVI）
SELECT
fp.parcel_code,
fp.geom,
nd.image_date,
nd.ndvi_mean,
CASE
WHEN nd.ndvi_mean > 0.7 THEN '健康'
WHEN nd.ndvi_mean > 0.5 THEN '正常'
WHEN nd.ndvi_mean > 0.3 THEN '注意'
ELSE '异常'
END AS health_status,
pr.crop_type
FROM farmland_parcels fp
JOIN planting_records pr ON fp.id = pr.parcel_id
JOIN ndvi_data nd ON fp.id = nd.parcel_id
WHERE pr.harvest_date IS NULL -- 当前种植中
AND nd.image_date = (SELECT MAX(image_date) FROM ndvi_data WHERE parcel_id = fp.id);
18.7 应急管理应用案例
18.7.1 灾害影响评估
灾害数据模型
-- 灾害事件表
CREATE TABLE disaster_events (
id SERIAL PRIMARY KEY,
event_code VARCHAR(50),
event_type VARCHAR(50), -- 洪水、地震、台风等
event_level VARCHAR(20), -- 特大、重大、较大、一般
start_time TIMESTAMP,
end_time TIMESTAMP,
affected_area GEOMETRY(POLYGON, 4326),
center_point GEOMETRY(POINT, 4326)
);
-- 避难场所表
CREATE TABLE shelters (
id SERIAL PRIMARY KEY,
name VARCHAR(100),
address VARCHAR(200),
capacity INTEGER,
shelter_type VARCHAR(50),
facilities TEXT[],
status VARCHAR(20) DEFAULT 'available',
geom GEOMETRY(POINT, 4326)
);
-- 救援队伍表
CREATE TABLE rescue_teams (
id SERIAL PRIMARY KEY,
team_name VARCHAR(100),
team_type VARCHAR(50), -- 消防、医疗、搜救
personnel_count INTEGER,
equipment TEXT[],
current_location GEOMETRY(POINT, 4326)
);
-- 受灾人口统计
CREATE TABLE affected_population (
id SERIAL PRIMARY KEY,
event_id INTEGER REFERENCES disaster_events(id),
area_name VARCHAR(100),
total_population INTEGER,
evacuated INTEGER,
injured INTEGER,
missing INTEGER,
geom GEOMETRY(POLYGON, 4326)
);
应急分析功能
-- 计算灾害影响范围内的资源
CREATE OR REPLACE FUNCTION analyze_disaster_impact(
p_event_id INTEGER
)
RETURNS TABLE (
metric_name VARCHAR,
metric_value FLOAT,
details JSONB
) AS $$
DECLARE
disaster_area GEOMETRY;
BEGIN
SELECT affected_area INTO disaster_area
FROM disaster_events WHERE id = p_event_id;
-- 受影响人口
RETURN QUERY
SELECT
'受影响人口'::VARCHAR,
COALESCE(SUM(pg.total_population), 0)::FLOAT,
jsonb_build_object('unit', '人')
FROM population_grid pg
WHERE ST_Intersects(pg.geom, disaster_area);
-- 受影响面积
RETURN QUERY
SELECT
'受影响面积'::VARCHAR,
ST_Area(disaster_area::geography) / 1000000.0,
jsonb_build_object('unit', '平方公里');
-- 可用避难场所
RETURN QUERY
SELECT
'周边避难场所'::VARCHAR,
COUNT(*)::FLOAT,
jsonb_agg(jsonb_build_object(
'name', s.name,
'capacity', s.capacity,
'distance_km', ST_Distance(s.geom::geography,
(SELECT center_point FROM disaster_events WHERE id = p_event_id)::geography) / 1000.0
))
FROM shelters s
WHERE ST_DWithin(s.geom::geography, disaster_area::geography, 10000)
AND s.status = 'available';
-- 可调配救援队伍
RETURN QUERY
SELECT
'可调配救援队伍'::VARCHAR,
COUNT(*)::FLOAT,
jsonb_agg(jsonb_build_object(
'team_name', rt.team_name,
'team_type', rt.team_type,
'personnel', rt.personnel_count
))
FROM rescue_teams rt
WHERE ST_DWithin(rt.current_location::geography, disaster_area::geography, 50000);
END;
$$ LANGUAGE plpgsql;
-- 疏散路线规划
WITH affected_areas AS (
SELECT geom FROM affected_population WHERE event_id = 1
),
nearest_shelters AS (
SELECT DISTINCT ON (aa.geom)
aa.geom AS from_area,
s.id AS shelter_id,
s.name AS shelter_name,
s.geom AS shelter_geom,
ST_Distance(ST_Centroid(aa.geom)::geography, s.geom::geography) AS distance_m
FROM affected_areas aa
CROSS JOIN shelters s
WHERE s.status = 'available'
ORDER BY aa.geom, ST_Distance(ST_Centroid(aa.geom)::geography, s.geom::geography)
)
SELECT
shelter_name,
COUNT(*) AS assigned_areas,
SUM(ap.total_population) AS total_to_evacuate,
ST_MakeLine(ST_Centroid(ns.from_area), ns.shelter_geom) AS evacuation_route
FROM nearest_shelters ns
JOIN affected_population ap ON ST_Equals(ns.from_area, ap.geom)
GROUP BY ns.shelter_id, ns.shelter_name, ns.shelter_geom;
18.8 数据迁移与 ETL 案例
18.8.1 Shapefile 批量导入
-- 创建导入日志表
CREATE TABLE import_logs (
id SERIAL PRIMARY KEY,
file_name VARCHAR(200),
table_name VARCHAR(100),
record_count INTEGER,
import_time TIMESTAMP DEFAULT NOW(),
status VARCHAR(20),
error_message TEXT
);
批量导入脚本（Shell）
#!/bin/bash
# 批量导入 Shapefile 到 PostGIS
DB_HOST="localhost"
DB_PORT="5432"
DB_NAME="gis_db"
DB_USER="postgres"
SHAPEFILE_DIR="/data/shapefiles"
SRID="4326"
ENCODING="GBK"
for shp in "$SHAPEFILE_DIR"/*.shp; do
filename=$(basename "$shp" .shp)
table_name=$(echo "$filename" | tr '[:upper:]' '[:lower:]' | tr ' ' '_')
echo "Importing $filename..."
shp2pgsql -s $SRID -I -D -W "$ENCODING" "$shp" "public.$table_name" | \
psql -h $DB_HOST -p $DB_PORT -d $DB_NAME -U $DB_USER
if [ $? -eq 0 ]; then
psql -h $DB_HOST -p $DB_PORT -d $DB_NAME -U $DB_USER -c \
"INSERT INTO import_logs (file_name, table_name, record_count, status)
SELECT '$filename', '$table_name', COUNT(*), 'success' FROM $table_name;"
else
psql -h $DB_HOST -p $DB_PORT -d $DB_NAME -U $DB_USER -c \
"INSERT INTO import_logs (file_name, table_name, status, error_message)
VALUES ('$filename', '$table_name', 'failed', 'Import error');"
fi
done
18.8.2 数据质量检查
-- 创建数据质量检查函数
CREATE OR REPLACE FUNCTION check_spatial_data_quality(
p_table_name VARCHAR,
p_geom_column VARCHAR DEFAULT 'geom'
)
RETURNS TABLE (
check_name VARCHAR,
issue_count INTEGER,
sample_ids INTEGER[],
description TEXT
) AS $$
BEGIN
-- 检查空几何
RETURN QUERY EXECUTE format('
SELECT
''空几何''::VARCHAR,
COUNT(*)::INTEGER,
array_agg(id ORDER BY id LIMIT 10),
''几何列为空''::TEXT
FROM %I
WHERE %I IS NULL
', p_table_name, p_geom_column);
-- 检查无效几何
RETURN QUERY EXECUTE format('
SELECT
''无效几何''::VARCHAR,
COUNT(*)::INTEGER,
array_agg(id ORDER BY id LIMIT 10),
''几何对象不符合 OGC 规范''::TEXT
FROM %I
WHERE NOT ST_IsValid(%I)
', p_table_name, p_geom_column);
-- 检查空 SRID
RETURN QUERY EXECUTE format('
SELECT
''缺失 SRID''::VARCHAR,
COUNT(*)::INTEGER,
array_agg(id ORDER BY id LIMIT 10),
''空间参考系统未设置''::TEXT
FROM %I
WHERE ST_SRID(%I) = 0
', p_table_name, p_geom_column);
-- 检查重复几何
RETURN QUERY EXECUTE format('
WITH duplicates AS (
SELECT id, %I,
ROW_NUMBER() OVER (PARTITION BY ST_AsBinary(%I) ORDER BY id) AS rn
FROM %I
)
SELECT
''重复几何''::VARCHAR,
COUNT(*)::INTEGER,
array_agg(id ORDER BY id LIMIT 10),
''存在完全相同的几何对象''::TEXT
FROM duplicates
WHERE rn > 1
', p_geom_column, p_geom_column, p_table_name);
-- 检查自相交多边形
RETURN QUERY EXECUTE format('
SELECT
''自相交多边形''::VARCHAR,
COUNT(*)::INTEGER,
array_agg(id ORDER BY id LIMIT 10),
''多边形边界自相交''::TEXT
FROM %I
WHERE GeometryType(%I) IN (''POLYGON'', ''MULTIPOLYGON'')
AND NOT ST_IsSimple(%I)
', p_table_name, p_geom_column, p_geom_column);
END;
$$ LANGUAGE plpgsql;
-- 使用示例
SELECT * FROM check_spatial_data_quality('poi', 'geom');
-- 批量修复无效几何
UPDATE my_table
SET geom = ST_MakeValid(geom)
WHERE NOT ST_IsValid(geom);
18.9 性能优化最佳实践
18.9.1 查询优化技巧
-- 1. 使用空间索引过滤
-- 差：全表扫描
SELECT * FROM poi WHERE ST_Distance(geom::geography, target::geography) < 1000;
-- 好：使用索引
SELECT * FROM poi WHERE ST_DWithin(geom::geography, target::geography, 1000);
-- 2. 避免在大表上使用 ST_Union
-- 差：一次性合并所有几何
SELECT ST_Union(geom) FROM large_polygons;
-- 好：分批合并
WITH batched AS (
SELECT
(row_number() OVER ()) / 1000 AS batch,
geom
FROM large_polygons
),
batch_unions AS (
SELECT ST_Union(geom) AS geom
FROM batched
GROUP BY batch
)
SELECT ST_Union(geom) FROM batch_unions;
-- 3. 使用部分索引
CREATE INDEX idx_poi_restaurant_geom
ON poi USING GIST (geom)
WHERE category = '餐饮';
-- 4. 使用 CLUSTER 优化空间查询
CLUSTER poi USING idx_poi_geom;
-- 5. 并行查询配置
SET max_parallel_workers_per_gather = 4;
SET parallel_tuple_cost = 0.001;
SET parallel_setup_cost = 10;
18.9.2 数据分区策略
-- 按地理区域分区
CREATE TABLE poi_partitioned (
id SERIAL,
name VARCHAR(100),
category VARCHAR(50),
geom GEOMETRY(POINT, 4326),
region VARCHAR(20)
) PARTITION BY LIST (region);
CREATE TABLE poi_north PARTITION OF poi_partitioned FOR VALUES IN ('北京', '天津', '河北');
CREATE TABLE poi_east PARTITION OF poi_partitioned FOR VALUES IN ('上海', '江苏', '浙江');
CREATE TABLE poi_south PARTITION OF poi_partitioned FOR VALUES IN ('广东', '福建', '海南');
-- 自动分区函数
CREATE OR REPLACE FUNCTION poi_region_trigger()
RETURNS TRIGGER AS $$
BEGIN
-- 根据经度判断区域
IF ST_X(NEW.geom) > 118 THEN
NEW.region := '上海';
ELSIF ST_Y(NEW.geom) > 35 THEN
NEW.region := '北京';
ELSE
NEW.region := '广东';
END IF;
RETURN NEW;
END;
$$ LANGUAGE plpgsql;
CREATE TRIGGER set_poi_region
BEFORE INSERT ON poi_partitioned
FOR EACH ROW EXECUTE FUNCTION poi_region_trigger();
18.10 本章小结
本章通过多个实际案例展示了 PostGIS 在不同领域的高级应用：
- 智慧城市：设施选址分析、城市热力图
- 交通物流：配送路线优化、车辆轨迹分析
- 环境监测：空气质量监测、水质监测
- 商业分析：商业选址、客户聚类分析
- 精准农业：地块管理、产量分析
- 应急管理：灾害影响评估、疏散规划
- 数据迁移：ETL 流程、数据质量检查
- 性能优化：查询优化、分区策略
18.11 总结与展望
通过本教程的学习，您已经掌握了 PostGIS 的核心概念和实际应用技能：
基础篇
- PostGIS 概述与安装配置
- 空间数据类型与空间参考系统
- 空间索引与性能优化
函数篇
- 几何构造、访问、输出函数
- 空间关系与分析函数
- 空间测量函数
高级篇
- 栅格数据处理
- 拓扑处理
- 三维与曲线几何
集成篇
- 与 GeoServer、QGIS 集成
- 编程语言集成
实战篇
- 多领域应用案例
- 性能优化最佳实践
PostGIS 作为开源空间数据库的领导者，持续保持活跃的发展。未来值得关注的方向包括：
- 云原生支持：更好的容器化和 Kubernetes 部署
- 实时处理：流式空间数据处理能力
- 机器学习集成：空间数据与 AI/ML 的结合
- 更强的栅格能力：遥感数据处理优化
- 性能持续提升：并行计算和向量化优化
希望本教程能够帮助您在 GIS 开发中充分发挥 PostGIS 的强大功能！
相关资源：
浙公网安备 33010602011771号