PostGIS 计算长度和面积

参考官网：http://postgis.net/docs/manual-3.0/using_postgis_dbmanagement.html

相关函数参考：http://postgis.net/docs/manual-3.0/reference.html#Spatial_Relationships

 

The basis for the PostGIS geometry type is a plane. 
The shortest path between two points on the plane is a straight line. 
That means calculations on geometries (areas, distances, lengths, intersections, etc) 
can be calculated using cartesian mathematics and straight line vectors.

The basis for the PostGIS geographic type is a sphere. 
The shortest path between two points on the sphere is a great circle arc. 
That means that calculations on geographies (areas, distances, lengths, intersections, etc) 
must be calculated on the sphere, using more complicated mathematics. 
For more accurate measurements,
 the calculations must take the actual spheroidal shape of the world into account.

Because the underlying mathematics is much more complicated, 
there are fewer functions defined for the geography type than for the geometry type.
 Over time, 
as new algorithms are added, the capabilities of the geography type will expand.

 

Geography Basics

CREATE TABLE ptgeogwgs(gid serial PRIMARY KEY, geog geography(POINT) );

CREATE TABLE ptzgeogwgs84(gid serial PRIMARY KEY, geog geography(POINTZ,4326) );

CREATE TABLE global_points (
    id SERIAL PRIMARY KEY,
    name VARCHAR(64),
    location GEOGRAPHY(POINT,4326)
  );

查看表的地理字段类型

select srid from geometry_columns where f_table_name='poi';

select srid from geography_columns where f_table_name='poi';

我平时导入Pg的表都是 geometry类型，srid = 4326。

如果需要计算长度和面积，需要转换地理坐标系为投影坐标系。如何选择投影坐标系，需要根据数据的范围和大小决定。

参考：https://www.jianshu.com/p/be5049ad8884

查看POSTGIS里面的包含的坐标系

SELECT * FROM spatial_ref_sys where srid = 4326;

SELECT * FROM spatial_ref_sys where srid = 4527;

 

1. 计算长度

-- Distance calculation using GEOGRAPHY (122.2km)
SELECT ST_Distance('LINESTRING(-122.33 47.606, 0.0 51.5)'::geography, 'POINT(-21.96 64.15)'::geography);

-- Distance calculation using GEOMETRY (13.3 "degrees")
SELECT ST_Distance('LINESTRING(-122.33 47.606, 0.0 51.5)'::geometry, 'POINT(-21.96 64.15)'::geometry);



// 单位是米  59712.49434085
select st_distance(ST_GeomFromText('POINT(113.258 23.369)', 4326)::geography,ST_GeomFromText('POINT(113.789 23.1454)', 4326)::geography);

// 单位是米 59714.426548488445
select st_distance(ST_Transform(ST_GeomFromText('POINT(113.258 23.369)', 4326),4547),ST_Transform(ST_GeomFromText('POINT(113.789 23.1454)', 4326),4547));


// 单位是度  0.5761579297380248
select st_distance(ST_GeomFromText('POINT(113.258 23.369)', 4326),ST_GeomFromText('POINT(113.789 23.1454)', 4326));

 

2.查询 某一个点所在的区域数据  ，用st_intersects 函数

select mc,dm,level from xzqh_pg where st_intersects(geom,'srid=4326;point(113.258 23.365)');

# 空间连接
select count(*),t.mc 
    from pcs as t
    join jq as tt
    on ST_Contains(t.geom,tt.geom)
    group by t.mc
    order by count(*) desc;

ST_Contains()函数 换成 ST_Intersect()也可以。

3. 查看某一个点范围内的数据

// 0.0008度范围内的poi
select id ,geom ,name from poi where st_dwithin(geom,'srid=4326;point(113.250 23.173)',0.0008);

// 投影坐标系，传入的单位是米
select id ,geom ,name from poi where st_dwithin(ST_Transform(geom,3857),'srid=3857;point(12608698 2646774)',50);

 

4. 计算面积

// 单位是平方公里  3857投影 计算出来的面积明显变大了  8700多
select st_area(ST_Transform(geom,3857))/1000000 as area,mc from xzqh_pg where level ='3';


// 单位是平方公里  4527 国家2000投影坐标系 面积7357平方公里
select st_area(ST_Transform(geom,4527))/1000000 as area,mc from xzqh_pg where level ='3';


// 广州地区 使用使用4547 ，国家2000投影坐标系中的114E带。
select st_area(ST_Transform(geom,4547))/1000000 as area,mc from xzqh_pg where level ='3';

**********
4547 投影坐标系的参数说明：
CGCS2000 / 3-degree Gauss-Kruger CM 114E
WKT
PROJCRS["CGCS2000 / 3-degree Gauss-Kruger CM 114E",
    BASEGEOGCRS["China Geodetic Coordinate System 2000",
        DATUM["China 2000",
            ELLIPSOID["CGCS2000",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4490]],
    CONVERSION["3-degree Gauss-Kruger CM 114E",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",114,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["China - 112.5°E to 115.5°E onshore"],
        BBOX[21.52,112.5,45.45,115.5]],
    ID["EPSG",4547]]
Proj4
+proj=tmerc +lat_0=0 +lon_0=114 +k=1 +x_0=500000 +y_0=0 +ellps=GRS80 +units=m +no_defs
范围
112.50, 21.52, 115.50, 45.45