如何在postgis中创建一个以米为单位的圆圈？

Question

我想问如何用radius=4km创建一个圈子。我尝试了ST_Buffer 函数，但它创建了一个更大的圆圈。 （我通过将多边形插入新的 kml 文件来查看创建的圆。）

这就是我正在尝试的。

INSERT INTO camera(geom_circle) VALUES(geometry(ST_Buffer(georgaphy(ST_GeomFromText('POINT(21.304116745663165 38.68607570952619)')), 4000)))

圆的中心是一个 lon lat 点，但我不知道它的 SRID，因为我是从 kml 文件中导入它的。 我需要SRID 来转换几何图形等吗？

Answer 1

KML 文件始终是纬度/经度并使用 SRID=4326。如果您使用 geography，则隐含此 SRID。地理是在纬度/经度数据上混合 4 公里公制度量的好方法……你试过这个太好了！

试试这个语句来修复强制转换，并使用参数化点构造函数：

SELECT ST_Buffer(ST_MakePoint(21.304116745663165, 38.68607570952619)::geography, 4000);

如果您需要将其转换回几何，请在末尾添加 ::geometry 转换。

---

准确性更新

上一个答案在内部将几何图形（通常）重新投影到该点适合的 UTM 区域（请参阅ST_Buffer）。如果该点位于两个 UTM 边界的边缘，这可能会导致轻微失真。大多数人不会关心这些错误的大小，但通常会是几米。但是，如果您需要亚毫米精度，请考虑构建动态azimuthal equidistant projection。这需要PostGIS 2.3的ST_Transform，改编自another answer：

CREATE OR REPLACE FUNCTION geodesic_buffer(geom geometry, dist double precision,
                                           num_seg_quarter_circle integer)
  RETURNS geometry AS $$
  SELECT ST_Transform(
    ST_Buffer(ST_Point(0, 0), $2, $3),
      ('+proj=aeqd +x_0=0 +y_0=0 +lat_0='
       || ST_Y(ST_Centroid($1))::text || ' +lon_0=' || ST_X(ST_Centroid($1))::text),
      ST_SRID($1))
  $$ LANGUAGE sql IMMUTABLE STRICT COST 100;
CREATE OR REPLACE FUNCTION geodesic_buffer(geom geometry, dist double precision)
  RETURNS geometry AS 'SELECT geodesic_buffer($1, $2, 8)'
  LANGUAGE sql IMMUTABLE STRICT COST 100;
-- Optional warppers for geography type
CREATE OR REPLACE FUNCTION geodesic_buffer(geog geography, dist double precision)
  RETURNS geography AS 'SELECT geodesic_buffer($1::geometry, $2)::geography'
LANGUAGE sql IMMUTABLE STRICT COST 100;
CREATE OR REPLACE FUNCTION geodesic_buffer(geog geography, dist double precision,
                                           num_seg_quarter_circle integer)
  RETURNS geography AS 'SELECT geodesic_buffer($1::geometry, $2, $3)::geography'
  LANGUAGE sql IMMUTABLE STRICT COST 100;

运行其中一个功能的简单示例是：

SELECT geodesic_buffer(ST_MakePoint(21.304116745663165, 38.68607570952619)::geography, 4000);

为了比较每个缓冲点的距离，这里是每个geodesic 的长度（旋转椭圆体上的最短路径，即 WGS84）。首先这个函数：

SELECT count(*), min(buff_dist), avg(buff_dist), max(buff_dist)
FROM (
  SELECT ST_Distance((ST_DumpPoints(geodesic_buffer(poi, dist)::geometry)).geom, poi) AS buff_dist
  FROM (SELECT ST_MakePoint(21.304116745663165, 38.68607570952619)::geography AS poi, 4000 AS dist) AS f
) AS f;

 count |      min       |       avg       |      max
-------+----------------+-----------------+----------------
    33 | 3999.999999953 | 3999.9999999743 | 4000.000000001

将此与 ST_Buffer（答案的第一部分）进行比较，这表明它偏离了大约 1.56 m：

SELECT count(*), min(buff_dist), avg(buff_dist), max(buff_dist)
FROM (
  SELECT ST_Distance((ST_DumpPoints(ST_Buffer(poi, dist)::geometry)).geom, poi) AS buff_dist
  FROM (SELECT ST_MakePoint(21.304116745663165, 38.68607570952619)::geography AS poi, 4000 AS dist) AS f
) AS f;

 count |      min       |       avg        |      max
-------+----------------+------------------+----------------
    33 | 4001.560675049 | 4001.56585986067 | 4001.571105793