这是PostgreSQL 和 PostGIS的PostGIS部分,主要讲PostgreSQL的空间相关功能的使用
postgres 账户下,修改 ~/.profile
,加入如下:
export POSTGIS_ENABLE_OUTDB_RASTERS=1
export POSTGIS_GDAL_ENABLED_DRIVERS=ENABLE_ALL
退出postgres账户,重新进入,输入如下命令检查设置是否成功:
printenv | grep POSTGIS
一般而言,Windows在安装时会提示增加这些环境变量,在安装时点“是”这些配置就已经完成了。
如果你点了否,开始菜单搜索“编辑系统环境变量”进行相应的修改。
使用如下SQL语句,在数据库中创建空间扩展
CREATE EXTENSION IF NOT EXISTS plpgsql;
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_raster; -- OPTIONAL
CREATE EXTENSION postgis_topology; -- OPTIONAL
Geometry和Geograpy是PostGIS空间类型的两大类数据类型,异同点为:
- Geometry是定义在平面上的空间图形,两点之间的路径是直线,在笛卡尔坐标系中,点之间、点线的距离,空间关系,如交叉等很容易算出;而Geography是 (类)球面 的上的坐标,两点之间最短路径是大圆弧。因而距离、面积、空间关系等都需要在球面(类球面)上进行,其复杂度远超前者。
- Geometry支持的函数Geography多很多。当然,随着时间经过,后者功能也在快速增加补全。
- Geometry类型的同样参数的计算会比geography快很多
- Geometry与Geography可以通过地图学的投影相互转换,只要是基于经纬度的空间参考,只要能在表
spatial_ref_sys
中找到,就可以实现转换。 - 无论你使用哪种空间参考,Geography的空间函数
ST_Distance
,ST_Length
,ST_Perimeter
,ST_Area
返回的结果和ST_DWithin
都以米为基础单位。 - 标准OCG所声明的类型,包括点、多点、线、多线、多边形等等,除了曲线,Geography都支持
- 如果你的数据在一个小的区域内,从性能和功能的角度考虑,你应该找一个合适的投影并使用Geometry类型
- 如果你的数据是全球的或者横跨一个大洲,使用Geography可以免去投影方面的麻烦,将数据存储为经纬度格式,使用Geography相关的函数来轻松计算距离等
- 如果你完全不懂投影也不想学,准备使用Geography并接受他的一些限制
- Geography在默认在
类球面
上计算,Geography函数都有一个选项可以关闭这个选项以使用球面
来简化和加快计算速度。 - 日界线和极地没有特殊考虑。
- 大圆弧计算最短那一条圆弧。
- 计算欧洲、俄罗斯面积或者插入大的地理区域时运行非常缓慢。因为多边形的点太多了。参照
ST_Subdivide
函数的文档来分而治之。 - 你甚至可以自己定义其他类球面星球的参考系统
pgAdmin的导入十分不好用,这里使用命令行工具shp2pgsql。
导入之前,需要获取获取数据的参考等信息,主要是EPSG码
命令行查看:
ogrinfo -al -so SHAPEFILE.shp | grep EPSG | tail -n 1 | grep -Eo '[0-9]{4,8}'
Python程序提取
# Save this file as get_epsg_code.py
import sys
import os
from osgeo import ogr
def get_epsg(shp_path):
driver = ogr.GetDriverByName("ESRI Shapefile")
shp = driver.Open(shp_path, 0)
layer = shp.GetLayer()
sr = layer.GetSpatialRef()
if sr.GetAuthorityName(None) == "EPSG":
return sr.GetAuthorityCode(None)
else:
return '0'
if __name__ == "__main__":
shp_path = sys.argv[-1]
shp_path = os.path.realpath(shp_path)
print(get_epsg(shp_path))
以下面的命令为模板:
shp2pgsql -s EPSG_ID -d SHPFILE.shp public.TABLENAME | psql -d DBNAME -U ROLENAME -c
栅格参考查看(除了主命令从矢量的 ogrinfo
改成了 gdalinfo
之外,栅格会出现很多个 EPSG 信息,但是最后一个才是我们需要的栅格投影信息):
gdalinfo GEO.TIFF | grep EPSG | tail -n 1 | grep -Eo '[0-9]{4,8}'