专栏名称: 朝阳35处
互联网金融大数据挖掘哪家强,平安前海征信帮你忙。分享数据挖掘和人工智能前沿技术,探讨其在风险控制、反欺诈等金融业务中的实际应用。每周涨点知识,多点谈资,定期举办行业沙龙。
目录
相关文章推荐
内蒙古检察  ·  媒体+检察丨检察建议助力换电柜安全升级 ·  21 小时前  
内蒙古检察  ·  媒体+检察丨检察建议助力换电柜安全升级 ·  21 小时前  
合肥市场监管  ·  合肥市发布2025年市级重点产品质量安全监管目录 ·  2 天前  
杭州网  ·  刚刚,阿里重磅宣布! ·  2 天前  
杭州网  ·  刚刚,阿里重磅宣布! ·  2 天前  
互联网坊间八卦  ·  刚刚,阿里宣布投入3800亿元建设云和AI硬 ... ·  2 天前  
51好读  ›  专栏  ›  朝阳35处

构建自己的地理信息空间数据库及与客户端简单交互

朝阳35处  · 公众号  ·  · 2018-06-07 10:38

正文

原创: 杜雨 数据小魔方

最近研究了下postgresql数据库及其空间地理信息拓展插件——postgis。

postgis作为新一代空间数据存储标准模型,将空间地理信息数据结构规范为关系型数据库可以承载的sp模式(simple features),这样,使得之前门槛颇高的gis空间数据存储模式变得通俗易懂、简单明了。

最重要的只要接触过SQL语言,就可以利用postgis的SQL语法便捷的操纵装载着空间信息的数据框(数据表),这些二维表除了被设定了一个特殊的空间地理信息字段(带有空间投影信息、经纬度信息等)之外,与主流数据管理系统所定义的各种字段并无两样。

本篇作为postgis数据库的一个前期探索篇,主要简单分享下postgresql+postgis的环境配置,及其与R语言、Python的API接口调用,以及如何通过这些接口来将shp、json空间地理信息数据源导入postgis。

1、环境配置篇:(可执行程序安装,如果你命令行比较熟练可以参照百度中的终端命令行进行安装)

关于postgis的环境配置,要先配置好postgresql环境,直接在以下主页下载安装即可:

https://www.enterprisedb.com/downloads

版本不要下载太高,建议9.6即可。

具体过程可以直接参考百度的教程:

https://www.yiibai.com/postgresql/install-postgresql.html

其中有几个细节点需要格外注意(自己踩过的坑)

1、尽量自己命名一个主目录(英文、不要带空格)
2、安装完postgresql之后会自定提示是否安装扩展插件(勾选postgis),如果这一步失败了不用担心,只是postgis没有安装成功,可以单独下载exe文件安装。
3、postgis安装(一定记得要和postgresql的主目录保持一致)
http://postgis.net/2017/07/01/postgis-2.3.3/
4、postgis安装之后会在postgresql库中新建一个带有空间数据表格式的模板库,此时使用postgresql安装环境中自带的pgAdmin4 工具打开postgresql数据库, 并可以新建一个引用空间数据表模板的测试库,这一步也有一个坑,在新建引用模板的测试库之后,一定要先按照官网给的步骤在测试库中运行以下脚本:

-- Enable PostGIS (includes raster)
CREATE EXTENSION postgis;
-- Enable Topology
CREATE EXTENSION postgis_topology;
-- Enable PostGIS Advanced 3D-- and other geoprocessing algorithms
-- sfcgal not available with all distributions
CREATE EXTENSION postgis_sfcgal;
-- fuzzy matching needed for Tiger
CREATE EXTENSION fuzzystrmatch;
-- rule based standardizer
CREATE EXTENSION address_standardizer;
-- example rule data set
CREATE EXTENSION address_standardizer_data_us;
-- Enable US Tiger Geocoder
CREATE EXTENSION postgis_tiger_geocoder;

5、上一步完成之后,即可通过postgis安装目录中的PostGIS 2.0 Shapefile and DBF Loader Exporter工具来手动导入本地的shp文件。

导入时要先建立与测试库的连接,并加载shp数据,含有中文要设置encoding = GBK。

显示导入成功即可刷新刚才的测试库,在测试库-schemas-public-tables中即可看到你新导入的控件数据集,与普通的数据库表并没有什么两样,仅仅是新增了一列叫做geom(geometry)的空间地理信息字段。


这张表整体就是我们之前在分享 R语言的sf对象和Python中的GeoDataFrame对象的技术雏形。

可以看到地理信息列在postgis中已经被编码成一组特殊数字,而在R中的sf对象中则是嵌套列表,在Python的GeoDataFrame中则是特殊的geomtry列。

  • geomtry对象一共分为7中,分别为:

  • point/mutipoint

  • string/mutistring

  • polygon/mutipolygon

  • CollectionFetures(前几种种的集合)


我们平时使用最多的地理信息多边形便是mutipolygon格式。

如果觉得pgAdmin4界面信息过于繁杂,可以安装Navicat Premium,它可以直接与postgresql数据库连接,作为一个桌面可视化管理界面。


Navicat Premium界面干净整洁,几乎没有任何冗余信息,具备常用的数据查询、管理功能,非常方便。

2、postgis与R语言通讯:

在R语言中调用postgis库表,需要依赖以下两个包(RPostgreSQL\rpostgis):

library("rpostgis")
library("RPostgreSQL")
library("sf")
library("ggplot2")
library("magrittr")

读取空间数据

conn #驱动名称
     dbname='mytest',  #要连接的库名称
     host='localhost', #本机地址
     port='5432',      #port编码
     user='postgres',  #用户名(在安装时默认生成,也可自定义,记清楚就好)
     password='*****'  #密码(自定义)
     )                 #建立连接池 
postgresqlpqExec(conn, "SET client_encoding = 'gbk'")  
#设定编码(多字节字符串)

以上两句联通了postgresql数据库表mytest,并设定了编码适应中文环境。

读入方法1

my_spdf #连接池名称
      name=c("public","bou2_4p"),    #指定schemas和表名,长度为2的向量,顺序不要乱
      geom = "geom"                  #指定表中的地理信息字段列名称
      ) %>% st_as_sf()               #导入数据默认为sp格式,转换为sf格式

读入方法2

map_data "bou2_4p")


方法一实在是太麻烦了,sf包的导入函数中封装了更加简便高效的导入函数:
直接指定连接池和测试库中空间数据表表明即可。

写入空间数据:

写入空间数据时,一般要以sp格式写入(就是之前用的最多的,maptools、rgdal包导入的默认格式),但是好在sf包中提供了一键转化sf和sp对象的函数,所以这里的写入数据格式转换非常高效。

world_map "shape/nc.shp", package="sf")) %>% 
 as("Spatial")

写入方法1

pgInsert(
     conn,name=c("public","world_map"),   #链接池名称+写入的表表明(自定义)
     data.obj = world_map                 #本地sp对象表
     )

写入方法2

sf包中也封装了直接写入postgis数据库的函数:

nc "shape/nc.shp", package="sf"))

st_write(
     nc,            #本地表名
     dsn = conn,    #指定连接池
     "world_data",  #写入后在库中名称
     layer_options = "OVERWRITE=true"
     )

写入之后在pgAdmin平台上刷新对应测试库之后即可看到新写入的表内容。


关闭链接

dbDisconnect(conn)

最后用刚才导入的表做一个简单的填充图。

ggplot() +
  geom_sf(data = map_data,aes(fill = perimeter)) +
  coord_sf(crs = 4326) +
  scale_fill_distiller(palette ='BrBG') +
  theme_void()







请到「今天看啥」查看全文