PostGIS中使用ST_Area计算面积的单位问题

owlion

已于 2023-10-14 14:48:34 修改

阅读量1.4k

点赞数 1

分类专栏： PostGIS 数据库开发文章标签：数据库 postgresql 空间计算

于 2023-10-14 14:30:06 首次发布

本文链接：https://blog.csdn.net/owlion/article/details/133823009

版权

数据库开发同时被 2 个专栏收录

3 篇文章 0 订阅

订阅专栏

PostGIS

2 篇文章 0 订阅

订阅专栏

1. 问题描述

Postgres 数据库中存储了一张 Polygon 类型的空间数据表，投影类型为 WGS 4326。想计算每个多边形的面积。但使用PostGIS提供的 ST_Area 函数返回结果非常小，明显不对。

SELECT
	ST_Area(
		ST_SetSRID(
			ST_GeomFromText('POLYGON((
				106.06615188500007 29.25636267300007,
				106.06616592500006 29.25771628100007,
				106.06770954400008 29.25770396000007,
				106.06769548400007 29.25635035300007,
				106.06615188500007 29.25636267300007))'
			),
			4326
		)
	);

# 返回结果
-----------------------
 2.089613822508709e-06
(1 row)

2. 问题原因

2.1 先查询ST_Area官方文档，说明如下：

float ST_Area(geometry g1);
float ST_Area(geography geog, boolean use_spheroid = true);
Description: Returns the area of a polygonal geometry. For geometry types a 2D Cartesian (planar) area is computed, with units specified by the SRID. For geography types by default area is determined on a spheroid with units in square meters. To compute the area using the faster but less accurate spherical model use ST_Area(geog,false).

返回多边形几何的面积。对于几何类型，计算二维笛卡尔(平面)面积，其单位由SRID指定。对于地理类型，默认情况下，区域是在以平方米为单位的球体上确定的。要使用更快但精度较低的球面模型计算面积，请使用ST_Area(geog,false)。
2.2 查询确认数据的空间参考( Spatial Reference ) 类型所使用的单位如下
```
SELECT srtext FROM spatial_ref_sys WHERE srid = (
	SELECT st_srid(geom) FROM my_data_table LIMIT 1
)


# 返回结果
srtext
-----------------------------------------------------------
GEOGCS["WGS 84",
	DATUM[
		"WGS_1984",
		SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
		AUTHORITY["EPSG","6326"]
	],
	PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
	UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
	AUTHORITY["EPSG","4326"]
]
```
- spatial_ref_sys 是PotGIS安装时默认自带的空间参考系数据表。
- 从返回结果看，EPSG：4326 只提供了地理坐标系 GEOGCS （Geographic Coordinate System）的参数信息，是在椭球体上按经纬度的度数进行描述的。并没有投影到平面坐标上。
- 从返回结果可见，EPSG：4326 空间参考系规定使用的单位是 度(UNIT: degree)，不是米，所以面积计算结果有问题。

3. 解决方案

3.1 在计算面积前，先将空间数据的参考系从 EPSG:4326 转换成其它使用 米(UNIT: metre) 为单位的空间参考系，然后再进行面积计算。
3.2 为了减少面积计算的误差，采用中国的 CGCS2000 投影来进行转换。2000坐标系在EPSG体系里的编号是 EPSG：4527，定义如下：

SELECT srtext FROM spatial_ref_sys srs WHERE srid = 4527
	
# 返回结果
srtext
-----------------------------------------------------------

PROJCS[
	"CGCS2000 / 3-degree Gauss-Kruger zone 39",
	GEOGCS["China Geodetic Coordinate System 2000",
		DATUM[
			"China_2000",
			SPHEROID["CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],
			AUTHORITY["EPSG","1043"]
		],
		PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
		UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
		AUTHORITY["EPSG","4490"]
	],
	PROJECTION["Transverse_Mercator"],
	PARAMETER["latitude_of_origin",0],
	PARAMETER["central_meridian",117],
	PARAMETER["scale_factor",1],
	PARAMETER["false_easting",39500000],
	PARAMETER["false_northing",0],
	UNIT["metre",1,AUTHORITY["EPSG","9001"]],
	AUTHORITY["EPSG","4527"]
]

从返回结果看，EPSG：4527 除了提供GEOCS地理坐标系的定义外，还提供了 地图投影(Projectioin) 的相关信息，是方便计算面积的平面坐标。且使用米作为计量单位。
PostGIS 里进行空间坐标系转换的函数是 ST_Transfor(geometry g1, int srid)，将几何对象g1的空间参考系转换成 srid 指定的参考系

4. 计算结果

SELECT
	ST_Area(
		ST_Transform(
			ST_SetSRID(
				ST_GeomFromText('POLYGON((
					106.06615188500007 29.25636267300007,
					106.06616592500006 29.25771628100007,
					106.06770954400008 29.25770396000007,
					106.06769548400007 29.25635035300007,
					106.06615188500007 29.25636267300007))'
				),
				4326
			),	4527
		)
	);


# 返回结果
-----------------------
23149.050374493825
(1 row)