PostGIS 空间数据学习建议
背景
我写过一些GIS的案例和文档，但是为了测试方便，文中大量使用了内置的几何point类型，并非GIS类型。
包括一些空间函数、空间数据的构建使用可能不是特别规范。
收到来自PostgreSQL社区GIS领域朋友的建议，为了防止给学习GIS的同学带来误导，请参考如下建议。
感谢这位朋友。
相关文章
相关建议
select DropGeometryColumn ('cb','geom');
drop table if exists cb;
create table cb(
objectid int not null, --共享单车编号
c1 int, -- 0表示未使用，其它表示已使用
c2 int, -- 共享单车属于哪家运营公司
constraint pk_cb_objectid primary key (objectid)
);
GIS必须明确指出地图单位是什么
民用GPS小数6位精度(6位大约为10米级精度,8位已经是亚米级精度)已经很高了,再多没有意义
c3 point – 共享单车当前位置, point在文章里只能算是自定义类型(实际上是PostgreSQL内置的几何点类型),这会给参考您文章的学习GIS的同学造成困扰
where c1=0 and c2=100 and c3 <@ circle ‘((23,3175),1000)’ order by c3 <-> point(23,3175) limit 1000;
搜索某个点附近1000距离内，属于某个公司的，没有使用的共享单车。
这样的查询条件在测试没问题,但是别人看了会造成困扰，因为没有地图单位,这个1000距离不知道是什么东西
新版本的postgis文档中已经没有ST_GeomFromText (‘POINT(121.403833486783 31.1425794813889)’, 4326)这样的函数了,虽然还支持,未来可能会删除
请使用ST_SetSRID(ST_GeomFromText (‘POINT(121.403833486783 31.1425794813889)’),4326)或ST_GeomFromText (‘SRID=4326;POINT(121.403833486783 31.1425794813889)’)
发现好几篇您写的关于gis的文章都有类似的问题,希望关于gis方面的文章严格按postgis标准
创建空间字段,根据空间参考不同,地图单位可能为米、度或其它，统称地图单位
select AddGeometryColumn ('cb','geom',4326,'POINT',2); -- 共享单车当前位置,GPS采用4326 ,类型为点,二维坐标
创建空间索引
create index gidx_cb_geom on cb using gist(geom);
或
create index gidx_cb_geomgraphy on cb using gist((geom::geography));
geometry和geography类型的区别
The geography type provides native support for spatial features represented on “geographic” coordinates (sometimes called “geodetic” coordinates, or “lat/lon”, or “lon/lat”). Geographic coordinates are spherical coordinates expressed in angular units (degrees).
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, and the mathematics becomes very complicated indeed.
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.
One restriction is that it only supports WGS 84 long lat (SRID:4326). It uses a new data type called geography. None of the GEOS functions support this new type. As a workaround one can convert back and forth between geometry and geography types.
The new geography type uses the PostgreSQL 8.3+ typmod definition format so that a table with a geography field can be added in a single step. All the standard OGC formats except for curves are supported.
经度范围：73°33′E至135°05′E
纬度范围：3°51′N至53°33′N
do $$
declare vStart bigint;
declare vEnd bigint;
declare MAXVALE bigint;
declare INTERVAL bigint;
begin
MAXVALE := 20000000;
INTERVAL := 1000; vStart := 1 ;vEnd := INTERVAL;
loop
-- 20家公司比较符合市场现状,更能反应实际情况
insert into cb
select id,(random()*1)::integer, (random()*(20-1)+1)::integer,
ST_SetSRID(ST_Point(
round((random()*(135.085831-73.406586)+73.406586)::numeric,6),
round((random()*(53.880950-3.408477)+3.408477)::numeric,6)
),4326)
from generate_series(vStart,vEnd) as id;
raise notice '%', vEnd;
vStart := vEnd + 1; vEnd := vEnd + INTERVAL;
if( vEnd > MAXVALE ) then
return;
end if;
end loop;
end$$;
ix,iy为gps经度和纬度(单位为度)
idistance为搜索距离(单位为米)
drop function if exists spatialQuery(ix float,iy float,idistance float);
create or replace function spatialQuery(ix float,iy float,idistance float)
returns table(oobjectid integer,oc1 integer,oc2 integer,odistance float,ogeom geometry)
as $$
declare
vrecord record;
vcurrentpoint geometry;
vspheroid spheroid;
begin
vspheroid := 'SPHEROID["WGS84",6378137,298.257223563]' ; --WGS84椭球体参数定义
vcurrentpoint := ST_SetSRID(ST_Point(ix,iy),4326); --
--查找圆心为vcurrentpoint,半径idistance米范围内未使用的共享单车,并按距离排序,只返回1千行
return query ( with cte as(
select * from cb
where ST_DWithin(geom::geography ,vcurrentpoint::geography,idistance,true)
) select objectid,c1,c2,ST_DistanceSpheroid(geom,vcurrentpoint,vspheroid),geom
from cte where c1=0 order by ST_DistanceSpheroid(geom,vcurrentpoint,vspheroid) limit 1000 );
end;
$$ language plpgsql;
select * from spatialQuery(102,24,5000);
查询计划
explain (analyze,verbose,costs,buffers,timing)
with cte as(
select * from cb
where ST_DWithin(geom::geography ,ST_SetSRID(ST_Point(102,24),4326)::geography,5000,true)
) select objectid,c1,c2,ST_DistanceSpheroid(geom,ST_SetSRID(ST_Point(102,24),4326),'SPHEROID["WGS84",6378137,298.257223563]'),geom
from cte where c1=0 order by ST_DistanceSpheroid(geom,ST_SetSRID(ST_Point(102,24),4326),'SPHEROID["WGS84",6378137,298.257223563]') limit 1000;
距离计算的建议：
别用平面坐标系，直接使用球坐标，用st_distancesphere函数计算球坐标表面的距离。
建议学习GIS的同学，多看GIS文档
http://postgis.net/docs/