1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

-- center 中心点
-- polygons 障碍物集合
-- radius 半径范围（米）
-- rays 射线条数
-- heading 面向方向（正东为0°） 0°~360°
-- fov 视角 0°~360°
CREATE OR REPLACE FUNCTION ISOVIST(
    IN center geometry,
    IN polygons geometry[],
    IN radius numeric DEFAULT 150,
    IN rays integer DEFAULT 36,
    IN heading integer DEFAULT 0,
    IN fov integer DEFAULT 360
    )
RETURNS geometry  AS $$
DECLARE
    arc numeric;
    angle_0 numeric;
    geomout geometry;
BEGIN
    arc := fov::numeric / rays::numeric;
    IF fov = 360 THEN
        angle_0 := 0;
    ELSE
        angle_0 := heading - 0.5 * fov;
    END IF;

    WITH
    buildings_0 AS(
        SELECT
            t.geom
        FROM unnest(polygons) as t(geom)
    ),
   buildings_crop AS(
        SELECT
            geom
        FROM buildings_0
        WHERE ST_DWithin(center::geography, geom::geography, radius)
    ),
    buildings AS(
        SELECT geom FROM buildings_crop
        UNION ALL
        SELECT ST_buffer(center::geography, radius)::geometry as geom
    ),
    rays AS(
        SELECT
            t.n as id,
            ST_SetSRID(
                ST_MakeLine(
                    center,
                    ST_Project(
                       center::geography,
                       radius + 1,
                       radians(angle_0 + t.n::numeric * arc)
                    )::geometry
                ),
             4326) AS geom
        FROM generate_series(0, rays) as t(n)
    ),
    intersections AS(
        SELECT
            r.id,
            (ST_Dump(ST_Intersection(ST_Boundary(b.geom),r.geom))).geom AS point
        FROM
            rays r
        LEFT JOIN
            buildings b
        ON
            ST_Intersects(b.geom,r.geom)
    ),
    intersections_distances AS(
        SELECT
            id,
            point as geom,
            row_number() over(partition by id order by center <-> point) as ranking
        FROM intersections
    ),
    intersection_closest AS(
        SELECT
            -1 as id,
            CASE WHEN fov = 360 THEN null::geometry ELSE center END as geom
        UNION ALL
        (SELECT
            id,
            geom
        FROM intersections_distances
        WHERE ranking = 1
        ORDER BY ID)
        UNION ALL
        SELECT
            999999 as id,
            CASE WHEN fov = 360 THEN null::geometry  ELSE center END as geom
    ),
    isovist_0 AS(
        SELECT
            ST_MakePolygon(ST_MakeLine(geom)) as geom
        FROM intersection_closest
    ),
    isovist_buildings AS(
        SELECT
            ST_CollectionExtract(ST_union(b.geom),3) as geom
        FROM
            isovist_0 i,
            buildings_crop b
        WHERE ST_Intersects(b.geom,i.geom)
    )
    SELECT
        coalesce(ST_Difference(i.geom, b.geom), i.geom) into geomout
    FROM
        isovist_0 i,
        isovist_buildings b;

    RETURN geomout;
END;
$$ language plpgsql IMMUTABLE;

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

-- 中心点 (114.42334 30.5136)
-- osm_hubei_building 障碍物表（此处用的是OSM的建筑物表）
-- 范围 500m
-- 射线条数 120
-- 面向方向 0°
-- 视角大小 360°
select ISOVIST(st_geomfromtext('point (114.42334 30.5136)', 4326), 
	array(select geom from osm_hubei_building),
	500, 120, 0, 360)

-- 优化 sql 运行时间：先获取范围内的建筑物，再进行视域分析
with buildings as (select geom from osm_hubei_building o where ST_Intersects(ST_buffer(
st_geomfromtext('point (114.42334 30.5136)', 4326)::geography, 1000)::geometry, o.geom)
)
select ISOVIST(st_geomfromtext('point (114.42334 30.5136)', 4326), 
	array(select geom from buildings),
	1000, 360, 0, 360)