批量將本地gis數據導入postgis資料庫

來自專欄 R語言數據分析與可視化11 人贊了文章

以前在處理gis數據的時候，都是直接導入本地shp素材、本地geojson素材，本地topojson素材，自從接觸postgis數據之後，深感使用規範的存儲系統來統一管理gis數據的好處，特別是數據量大了之後，優勢便更加明顯，你可以選擇將很多需要做空間計算的步驟轉移到Postgis資料庫內進行計算，要知道Postgis提供的空間計算能力與R和Python這種應用導向的工具相比，優勢要大得多。

在批量導入素材之前，我們可以先看下R語言目前提供的各種導入介面在I/O性能上相比有何異同。

#install.packages("geojsonio")#devtools::install_github("ropensci/geojsonio")library("geojsonio")library("rgdal")library("sf")library("maptools")

使用maptools包中的readShapePoly函數進行導入（已快被遺棄了，推薦使用sf和rgdal包）

system.time(china_map <- readShapePoly("D:/R/rstudy/CHN_adm/bou2_4p.shp"))用戶 系統 流逝 0.23 0.00 0.23 Warning message:use rgdal::readOGR or sf::st_read china_map@dataggplot2::fortify(china_map)

geojsonio包導入：

system.time(geojson1 <- geojson_read( "D:/R/rstudy/CHN_adm/bou2_4p.shp", method = "local", parse = TRUE, what = "sp", encoding="utf-8", use_iconv=TRUE ))用戶 系統 流逝 0.69 0.03 0.71

使用rgdal包：

system.time(map_data <- readOGR( "D:/R/rstudy/CHN_adm/bou2_4p.shp", encoding="utf-8", use_iconv=TRUE ))OGR data source with driver: ESRI Shapefile Source: "D:R
studyCHN_admou2_4p.shp", layer: "bou2_4p"with 925 featuresIt has 7 fieldsInteger64 fields read as strings: BOU2_4M_ BOU2_4M_ID 用戶 系統 流逝 0.66 0.09 0.75

使用sf包導入：

system.time(nepal_shp <- read_sf( "D:/R/rstudy/CHN_adm/bou2_4p.shp", options = "ENCODING=gbk" ))用戶 系統 流逝 0.05 0.00 0.05

可以看到在同一個shp文件單項導入的情況下，純粹從時間上來看：

sf > maptools > rgdal > geojsonio

這裡值得一提的是，geojsonio包是封裝的rgdal服務，性能上自然略遜rgdal一籌，以上四個包中，除sf包是基於simple features標準的模型之外，其他基本都是基於sp模型的。sf模型的性能由此可見一斑。

當然，以上sf包、rgdal包和sf包都是兼容性很好地包，可以支持非常廣泛的數據源，以下分別是在json標準下的兩種素材上進行測試。

geojson

system.time(geojson <- geojson_read( "D:/R/mapdata/State/china.geojson", method = "local", parse = TRUE, encoding="utf-8", use_iconv=TRUE, what = "sp" ))用戶 系統 流逝 0.80 0.02 0.81 system.time(map_data <- readOGR( "D:/R/mapdata/State/china.geojson", encoding="utf-8", use_iconv=TRUE, stringsAsFactors = FALSE ))OGR data source with driver: GeoJSON Source: "D:RmapdataStatechina.geojson", layer: "china"with 34 featuresIt has 2 fields用戶 系統 流逝 0.77 0.00 0.76 system.time(nepal_shp <- read_sf( "D:/R/mapdata/State/china.geojson" ))用戶 系統 流逝 0.03 0.00 0.03

topojson

system.time(map_data <- readOGR( "D:/R/mapdata/china.topojson", use_iconv=TRUE, encoding = "utf-8", stringsAsFactors = FALSE ))OGR data source with driver: GeoJSON Source: "D:Rmapdatachina.topojson", layer: "china"with 34 featuresIt has 2 fields用戶 系統 流逝 0.52 0.01 0.59 system.time(geojson <- topojson_read( "D:/R/mapdata/china.topojson", encoding="utf-8", use_iconv=TRUE ))OGR data source with driver: GeoJSON Source: "D:Rmapdatachina.topojson", layer: "china"with 34 featuresIt has 2 fields用戶 系統 流逝 0.59 0.00 0.59 system.time(nepal_shp <- read_sf( "D:/R/mapdata/china.topojson" ))用戶 系統 流逝 0.02 0.00 0.01

是不是看完這個性能大比拼之後大吃一驚，為sf包的超強IO能力所折服，sf包是一個非常強大的包，實現了基於simple features的所有特性，如果你了解一點兒Postgis的話，你會發現作者把大部分空間運算的函數名稱設計的和Postgis中的函數一模一樣，這就意味著你無論是只了解過sf包函數，或者只了解過Postgis函數，都可以低成本的遷移到兩一個平台，因為同名函數往往功能一致。

如果你要想將sf包導入的數據模型轉換為普通的數據框模式，僅僅只需使用其提供的as(sf,』Spatial』)函數一次轉化即可，當然sf有自己的ggplot2通道函數geom_sf()，這意味著你不必多此一舉。（當然對於sf不甚熟悉，習慣於使用geom_polygon來實現地理信息可視化的小夥伴兒，可以採取這種辦法，但是仍然要推薦大家學習sf包，因為它代表著未來）。

R語言-gis數據批量入庫：

#定義讀寫函數：task <- function(filename,conn){ #此處為寫入本地gis數據（可以是任意格式，可以使用任意一種導入工具） map_data <- readOGR(filename,use_iconv=TRUE,encoding = "utf-8",stringsAsFactors = FALSE) file_name <- sub(.json,,basename(filename)) #此處是寫入資料庫的函數，可以使用sf包、rgdal包以及RPostgreSQL包提供的寫出函數。 writeOGR(obj = map_data ,dsn = conn,driver = "PostgreSQL",layer=file_name,encoding="gbk",overwrite_layer = TRUE) } #此處使用l_ply函數創建批量執行任務Project_io <- function(path){ setwd(path) input_list = list.files(path) conn <- "PG:dbname=mytest host=localhost port=5432 user=postgres password=708965" l_ply(input_list,task,conn)}#啟動任務Project_io("D:/R/mapdata/Province")

Python-gis數據批量入庫：

import geopandas as gpdimport pandas as pdfrom sqlalchemy import create_enginefrom geoalchemy2 import Geometry,WKTElementimport numpy as npimport osimport reimport jsondef dict_to_json(data): for i in data: if isinstance (data[i].tolist()[0],dict): data[i] = data[i].map(lambda x:json.dumps(x)) return data#數據寫入函數：def write_gis(path): map_data = gpd.GeoDataFrame.from_file(path) map_data = dict_to_json(map_data) map_data[geometry] = map_data[geometry].apply(lambda x: WKTElement(x.wkt,4326)) map_data.drop([center,parent], axis = 1, inplace=True) map_data.to_sql( name = re.split(\.,path)[0], con = engine, if_exists= replace, dtype = {geometry:Geometry(geometry_type =POLYGON,srid = 4326)} ) return None#創建批量任務def to_do(file_path,username,password,dbname): os.chdir(file_path) link = "postgresql://{0}:{1}@localhost:5432/{2}".format(username,password,dbname) engine = create_engine(link,encoding = utf-8) file_list = os.listdir() map(lambda x: write_gis(x),file_list) return None#執行任務計劃if __name__ == __main__: file_path = D:/R/mapdata/Province username = postgres password = ***** dbname = mytest to_do(file_path,username,password,dbname) print(DODE)

推薦閱讀：