Choumingzhao/postgis_quick_notes.md

## postgis_quick_notes.md

      
    Raw
  

              postgis_quick_notes.md
            
          
    PostGIS中文速查速用手册

这是PostgreSQL 和 PostGIS的PostGIS部分，主要讲PostgreSQL的空间相关功能的使用
空间相关配置

开启栅格支持

postgres 账户下，修改 ~/.profile，加入如下：
export POSTGIS_ENABLE_OUTDB_RASTERS=1
export POSTGIS_GDAL_ENABLED_DRIVERS=ENABLE_ALL
退出postgres账户，重新进入，输入如下命令检查设置是否成功：
printenv | grep POSTGIS
Windows相关设置

一般而言，Windows在安装时会提示增加这些环境变量，在安装时点“是”这些配置就已经完成了。
如果你点了否，开始菜单搜索“编辑系统环境变量”进行相应的修改。
创建PostGIS扩展

使用如下SQL语句，在数据库中创建空间扩展
CREATE EXTENSION IF NOT EXISTS plpgsql;
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_raster; -- OPTIONAL
CREATE EXTENSION postgis_topology; -- OPTIONAL
PostGIS空间概念理解

Geometry 与 Geography概念辨析

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？


如果你的数据在一个小的区域内，从性能和功能的角度考虑，你应该找一个合适的投影并使用Geometry类型
如果你的数据是全球的或者横跨一个大洲，使用Geography可以免去投影方面的麻烦，将数据存储为经纬度格式，使用Geography相关的函数来轻松计算距离等
如果你完全不懂投影也不想学，准备使用Geography并接受他的一些限制

Geography常见问题


Geography在默认在 类球面 上计算，Geography函数都有一个选项可以关闭这个选项以使用球面来简化和加快计算速度。
日界线和极地没有特殊考虑。
大圆弧计算最短那一条圆弧。
计算欧洲、俄罗斯面积或者插入大的地理区域时运行非常缓慢。因为多边形的点太多了。参照ST_Subdivide函数的文档来分而治之。
你甚至可以自己定义其他类球面星球的参考系统

数据导入

矢量数据（如shp）

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
栅格数据（如geotiff）

栅格数据空间参考查询

栅格参考查看（除了主命令从矢量的 ogrinfo 改成了 gdalinfo 之外，栅格会出现很多个 EPSG 信息，但是最后一个才是我们需要的栅格投影信息）：
gdalinfo GEO.TIFF | grep EPSG | tail -n 1 | grep -Eo '[0-9]{4,8}'