標籤：

構建自己的地理信息空間資料庫及與客戶端簡單交互

06-07
構建自己的地理信息空間資料庫及與客戶端簡單交互

來自專欄 R語言中文社區

最近研究了下postgresql資料庫及其空間地理信息拓展插件——postgis。

postgis作為新一代空間數據存儲標準模型，將空間地理信息數據結構規範為關係型資料庫可以承載的sp模式（simple features），這樣，使得之前門檻頗高的gis空間數據存儲模式變得通俗易懂、簡單明了。

最重要的只要接觸過SQL語言，就可以利用postgis的SQL語法便捷的操縱裝載著空間信息的數據框（數據表），這些二維表除了被設定了一個特殊的空間地理信息欄位（帶有空間投影信息、經緯度信息等）之外，與主流數據管理系統所定義的各種欄位並無兩樣。

本篇作為postgis資料庫的一個前期探索篇，主要簡單分享下postgresql+postgis的環境配置，及其與R語言、Python的API介面調用，以及如何通過這些介面來將shp、json空間地理信息數據源導入postgis。

1、環境配置篇：（可執行程序安裝，如果你命令行比較熟練可以參照百度中的終端命令行進行安裝）

關於postgis的環境配置，要先配置好postgresql環境，直接在以下主頁下載安裝即可：

enterprisedb.com/downlo

版本不要下載太高，建議9.6即可。

具體過程可以直接參考百度的教程：

yiibai.com/postgresql/i

其中有幾個細節點需要格外注意（自己踩過的坑）

1、盡量自己命名一個主目錄（英文、不要帶空格）

2、安裝完postgresql之後會自定提示是否安裝擴展插件（勾選postgis）,如果這一步失敗了不用擔心，只是postgis沒有安裝成功，可以單獨下載exe文件安裝。

3、postgis安裝（一定記得要和postgresql的主目錄保持一致）

postgis.net/2017/07/01/

4、postgis安裝之後會在postgresql庫中新建一個帶有空間數據表格式的模板庫，此時使用postgresql安裝環境中自帶的pgAdmin4 工具打開postgresql資料庫，並可以新建一個引用空間數據表模板的測試庫，這一步也有一個坑，在新建引用模板的測試庫之後，一定要先按照官網給的步驟在測試庫中運行以下腳本：

-- Enable PostGIS (includes raster)CREATE EXTENSION postgis;-- Enable TopologyCREATE EXTENSION postgis_topology;-- Enable PostGIS Advanced 3D-- and other geoprocessing algorithms-- sfcgal not available with all distributionsCREATE EXTENSION postgis_sfcgal;-- fuzzy matching needed for TigerCREATE EXTENSION fuzzystrmatch;-- rule based standardizerCREATE EXTENSION address_standardizer;-- example rule data setCREATE EXTENSION address_standardizer_data_us;-- Enable US Tiger GeocoderCREATE 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
postgis）：library("rpostgis")library("RPostgreSQL")library("sf")library("ggplot2")library("magrittr")讀取空間數據conn <- dbConnect( PostgreSQL(), #驅動名稱 dbname=mytest, #要連接的庫名稱 host=localhost, #本機地址 port=5432, #port編碼 user=postgres, #用戶名（在安裝時默認生成，也可自定義，記清楚就好） password=***** #密碼（自定義） ) #建立連接池 postgresqlpqExec(conn, "SET client_encoding = gbk") #設定編碼(多位元組字元串)以上兩句聯通了postgresql資料庫表mytest，並設定了編碼適應中文環境。讀入方法1my_spdf <- pgGetGeom( conn, #連接池名稱 name=c("public","bou2_4p"), #指定schemas和表名，長度為2的向量，順序不要亂 geom = "geom" #指定表中的地理信息欄位列名稱 ) %>% st_as_sf() #導入數據默認為sp格式，轉換為sf格式讀入方法2map_data <- st_read(conn,"bou2_4p")讀入方法3conn <- "PG:dbname=mytest host=localhost port=5432 user=postgres password=******"mylayer <-readOGR(dsn= conn, layer = "bou2_4p",encoding = "utf-8",use_iconv = TRUE,stringsAsFactors = TRUE) #讀入方法3

方法一實在是太麻煩了，sf包的導入函數中封裝了更加簡便高效的導入函數：直接指定連接池和測試庫中空間數據表表明即可。寫入空間數據：寫入空間數據時，一般要以sp格式寫入（就是之前用的最多的，maptools、rgdal包導入的默認格式），但是好在sf包中提供了一鍵轉化sf和sp對象的函數，所以這裡的寫入數據格式轉換非常高效。world_map <- st_read(system.file("shape/nc.shp", package="sf")) %>% as("Spatial")寫入方法1pgInsert( conn,name=c("public","world_map"), #鏈接池名稱+寫入的表表明（自定義） data.obj = world_map #本地sp對象表 )寫入方法2sf包中也封裝了直接寫入postgis資料庫的函數：nc <- st_read(system.file("shape/nc.shp", package="sf"))st_write( nc, #本地表名 dsn = conn, #指定連接池 "world_data", #寫入後在庫中名稱 layer_options = "OVERWRITE=true" )寫入方法3conn <- "PG:dbname=mytest host=localhost port=5432 user=postgres password=708965"writeOGR( obj = mylayer, dsn = conn, driver = "PostgreSQL", layer = "mytest_map", encoding="gbk") 寫入之後在pgAdmin平台上刷新對應測試庫之後即可看到新寫入的表內容。

關閉鏈接

dbDisconnect(conn)

最後用剛才導入的表做一個簡單的填充圖。

ggplot() +

geom_sf(data = map_data,aes(fill = perimeter)) +

coord_sf(crs = 4326) +

scale_fill_distiller(palette =BrBG) +

theme_void()

3、postgis與Python交互：

python中的資料庫交互操作路徑比較多，這裡僅使用geopandas提供的寫入介面以及sqlalchemy寫出介面，探索出可行的讀寫代碼即可。

import geopandas as gpd

from sqlalchemy import create_engine

from geoalchemy2 import Geometry,WKTElement

建立連接

engine = create_engine(

"postgresql://postgres:******@localhost:5432/mytest",

encoding = "gbk"

)

create_engine函數創建客戶端與postgresql資料庫的鏈接，連接參數是一個包含資料庫驅動、用戶名、用戶密碼、本機環境：埠及測試庫名稱的長字元串。最後指定編碼（這裡涉及到中文字元串）

df = gpd.read_postgis("select * from bou2_4p",engine,crs = 4326)

使用geopandas包中提供的postgis介面函數，導入engine連接池mytest庫中的bou2_4p表所有數據。

china_map = gpd.GeoDataFrame.from_file("D:/R/rstudy/CHN_adm/bou2_4p.shp", encoding = gb18030)

china_map[geometry] = china_map[geometry].apply(lambda x: WKTElement(x.wkt,4326))

因為GeoDataFrame中默認存的simple features 對象與postgis庫中定義的默認simple features對象存在差異，這裡需要使用geoalchemy2包提供的Geometry,geoalchemy2介面將GeoDataFrame轉換為postgis庫可識別的模式。

最後利用pandas封裝的sqlalchemy寫出函數，將剛才規範過之後的表china_map寫入postgis庫中。

china_map.to_sql(

name = "world", #表名

con = engine, #連接池

if_exists= replace, #若已存在同名表，則替換

dtype = {

geometry:

Geometry(geometry_type =POLYGON,srid = 4326)}

#聲明地理信息欄位模式

)

本篇分享僅就postgis的安裝，與R和Python之間的簡單交互做一下梳理，之後如果有新的成果再做分享，期待各位小夥伴兒利用postgis結合R語言中的ggplot2+sf+shiny以及Python中的GeoDataFrame做出更多具備商業價值的東西，感興趣的也可以一起分享交流心得。

寫在最後！！！

最近陸陸續續整理了一些關於sf模型空間地理信息數據可視化得新的，特別是R中得sf包、Python中得geodataframe、basemap包等資料，這戲資料在中文互聯網上極其稀缺，陸陸續續整理了好幾天，逛了N多個Stack Overflow才打通這個流程，一定要珍惜哦，如果覺著這還不夠過癮，最近正在錄製的課程《R語言商務圖表與數據可視化》已經更新到第九章了，足足四章的地理信息可視化模型、原理、應用一定會讓你收穫滿滿，趕快來瞧瞧吧~

edu.hellobi.com/course/

閱讀原文


推薦閱讀：

※把英雄從網格領回街區
※珠峰有多高？（指從山腳下到頂峰的垂直距離）?
※經緯度和公里網之間，怎麼互相轉換？
※在武漢大學資環院就讀是一種什麼樣的體驗？
※觀念的水位

TAG:GIS地理信息系統 | 資料庫 | 地理信息 |