用 R 處理 PostgreSQL/PostGIS 資料 (1)

前言

前幾天寫的用 ggplot 畫 boxplot 並提及使用 RPostgreSQL 來連接 PostgreSQL查詢,今天我們試著來用更開放的標準來處理 PostgreSQL/PostGIS 資料庫,也就是用 ODBC (open database connectivity) 。另外,也簡單介紹 gdal (geospatial data abstraction library) 在 R 中處理 GIS 圖層的方式,以下都是在泛 unix 下執行,若你是 Windows 的使用者,你可以試試看,但我沒試過不保證會不會有其他的問題。

預習內容

SQL 語法

安裝及設定(版本都只是參考,你可以依據需求安裝)

  • PostgreSQL (版本為 9.3.4) -- 資料庫管理系統
  • PostGIS (版本為 2.1.3) -- 讓 PostgreSQL 支援空間資料庫
  • unixODBC -- unix ODBC
  • psqlodbc -- PostgreSQL 的 ODBC driver
  • gdal -- Geospatial Data Abstraction Library,處理空間資料(含向量及影像)的函式庫

如果你是用 MacOS X 的使用者,我會推薦你另外先安裝 homebrew(請先依照 homebrew 的指示把一些 developer tools 如 Xcode command line tools 等先安裝好。)

使用 homebrew 安裝上述的軟體

$ brew install postgresql postgis unixodbc psqlodbc

因為我們希望 gdal 能夠支援 postgresql,以及完整安裝,所以我們用下面的參數來安裝 gdal:

$ brew install gdal --complete --postgresql

如果你是 FreeBSD 的使用者,而且你使用的是 9.0 以上的版本,我會推薦使用 pkgng (即 pkg ),但 postgresql-odbc 預設 pkg 版本是 9.2,所以你可能要從 ports 安裝

# pkg install postgresql9 unixODBC postgis

# cd /usr/ports/databases/postgresql-odbc ; make install clean

另外,我們也希望裝完整的 gdal 格式支援,所以從 ports 安裝,並選擇所需要的格式

# cd /usr/ports/graphics/gdal; make install clean

如果你是 GNU/Linux Debian 或 Ubuntu 或 kFreeBSD 的使用者,你可以用 aptitude 或 apt-get 來安裝上述套件。

R 的相關套件安裝

以上的 R package 通常都會有 binary 可以安裝,但是若沒有提供則必須要從原始碼自己安裝起。以下舉 MacOS X Mavericks 從原始碼安裝為例:

# 從 source 安裝 RODBC,因為預設會用 /usr/lib 去找 MacOS X 內建的 iODBC,
# 但實際執行會有問題,所以我們希望在編譯的時候去使用 homebrew 所安裝的 unixodbc,
# 其預設目錄在 /usr/local/lib
install.packages('RODBC', type="source", configure.args = '--with-odbc-lib=/usr/local/lib')
,同時也指定 homebrew 所安裝的 gdal-config 目錄
# 從 source 安裝 rgdal ,編譯時所引用的函式庫也是 homebrew 所安裝的 gdal 和 proj
install.packages('rgdal', type='source', 
    configure.args='--with-gdal-config=/usr/local/bin/gdal-config 
                                    --with-proj-include=/usr/local/include')
# 從 source 安裝 RPostgreSQL
install.packages('RPostgreSQL', type='source')

另一個從 shell 安裝的方法是

R CMD INSTALL _pkgname_ --configure-args='_args_'

開始!

用 RPostgreSQL 查詢資料

library(RPostgreSQL)

drv = dbDriver("PostgreSQL")
# 建立連線,並使用 dbGetQuery() 取得查詢後的資料表
conn = dbConnect(drv, user= 'user', pwd='secret' dbname='database_name')
res = dbGetQuery(conn, "SELECT * FROM schema.table_name;")
# 或是用 dbSendQuery(),再用 fetch 取回資料
res2 = dbSendQuery(conn, "SELECT * FROM schema.table_name;")
res_df = fetch(res2, n=-1)

使用 dbGetQuery() 或是 dbSendQuery()/fetch() 時,取回來的值都是 DataFrame,所以就可以使用下列的方式來篩選資料

# 把 res 中欄位名稱為 ColumnName 中數值等於一的行(row)取出來
res[res$ColumnName==1,] 
# 取出 res 第一行
res[1,]
# 取出第一至第五列
res[,1:5]

用 RODBC 查詢資料

在使用 RODBC 時,必須要先設定好 odbcinst.ini 及 odbc.ini 這兩個 ODBC 設定檔,先在家目錄(/Users/yourusername 或 /home/yourusername,'~' 代表家目錄) 底下建立這兩個檔案

touch ~/.odbcinst.ini ~/.odbc.ini

接下來用編輯器編輯(vim/emacs/ee/gedit ...挑個喜歡的編輯器,內建的 TextEdit.app 也可以)這兩個設定檔。假設 postgresql odbc driver 是安裝在 /usr/local/lib 底下,在 odbcint.ini 檔案要加入 Driver (在本例中為 /usr/local/bin/psqlodbcw.so)路徑。

~/.odbcinst.ini
[PostgreSQL]
Description     = PostgreSQL driver 
Driver          = /usr/local/lib/psqlodbcw.so
FileUsage       = 1

一樣的,如果你的 PostgreSQL 安裝後沒有特別去更改什麼,就照下面設定稍微修改即可( your_database_name 填你實際的資料庫名稱,Username 預設用 postgres,如果你有另外新增 user ,則用你新增的 username,Protocol 則填寫你的 postgresql 版本主分支(ex: 9.3)

~/.odbc.ini
[your_database_name]
Driver = PostgreSQL
Database = your_database_name
Servername = localhost
Port = 5432
Username = postgres
Protocol = 9.3
ReadOnly = 0

如果上面都設定完成後,可以在終端機(Terminal)中輸入 isql -v (詳細輸出) dbname 來看能否正常用 ODBC 連接 PostgreSQL 資料庫,若出現下面畫面,代表已經成功了!

mutolisp@brachypodium ~ (master*) $ /usr/local/bin/isql -v dbname
+---------------------------------------+
| Connected!                            |
|                                       |
| sql-statement                         |
| help [tablename]                      |
| quit                                  |
|                                       |
+---------------------------------------+
SQL> 

你也可以直接寫一些 SQL 測試看看是否正確回傳你的查詢,例如

SELECT * FROM table WHERE condition LIMIT 10;

接下來的步驟則是在 RStudio(或 Rconsole 中)來使用 RODBC 查詢 PostgreSQL 資料庫

  1. 先確認 RODBC 的 odbc 資料來源是否存在,如果用 odbcDataSources() 查詢有出現上面你設定的資料庫名稱就代表順利連結 odbc 來源。
    > library(RODBC)
    > odbcDataSources()
      your_database_name 
    "PostgreSQL" 
    
  2. 接下來要建立和資料庫的連結
    > conn = odbcConnect("your_database_name")
    
  3. 查詢資料表,在這邊是使用 sqlQuery() 來查詢,如果之後都不會再用查詢了,記得使用 odbcClose() 關閉連結
    query = "SELECT * FROM table WHERE condition LIMIT 10;"
    res = sqlQuery(conn, query)
    odbcClose(conn)
    
    註:odbcQuery() 是 low-level ODBC function,如果你要使用 odbcQuery(),必須要先建立連接,查詢後再手動用 odbcFetchRows 抓回查詢回傳資料,但是回傳的資料型態是 list,和用 sqlQuery() 回傳的 DataFrame 是不同的,例如:
    conn = odbcConnect("your_database_name")
    query = "SELECT * FROM table WHERE condition;"
    qres = odbcQuery(conn, query)
    res = odbcFetchRows(conn, max=1000)
    

實際比較了一下 RPostgreSQL 和 RODBC 這兩種,RODBC 設定和使用上是比較麻煩的。RPostgreSQL 呼叫 DBI 去和不同的資料庫溝通,語法也相對上簡單,假設要和 SQLite 資料庫溝通,就可以直接安裝 RSQLite,作法也都是類似的

資料庫中繼介面 支援的資料庫 查詢概要
RPostgreSQL PostgreSQL dbConnect(drv), dbGetQuery()
RODBC 只要是支援 unixODBC,都可以使用 odbcConnect(), sqlQuery(), odbcClose()

用 rgdal 載入 PostGIS

  1. 先在 PostgreSQL 資料庫中建立 PostGIS 支援,用 "CREATE EXTENSION postgis;",相關匯入的作法也可以參考 http://mutolisp.logdown.com/posts/192197

    > library(RPostgreSQL)
    > library(rgdal)
    Loading required package: sp
    rgdal: version: 0.8-16, (SVN revision 498)
    Geospatial Data Abstraction Library extensions to R successfully loaded
    Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
    Path to GDAL shared files: /usr/local/Cellar/gdal/1.11.0/share/gdal
    Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
    Path to PROJ.4 shared files: (autodetected)
    > drv = dbDriver(drvName="PostgreSQL")
    > conn = dbConnect(drv, user='your_user_name', dbname='your_database_name')
    > dbSendQuery(conn, "create extension postgis;")
    
  2. 使用 ogrListLayers() 列出有空間幾何欄位的表格

    # 看一下支援的 driver 列表
    > ogrDrivers()
    # 先設定連結的 driver
    > dsn = "PG:dbname=nvdimp user=your_user_name password=secret"
    # 列出有空間幾何欄位的資料表
    > ogrListLayers(dsn) 
    [1] "nvdimp_plots"          "nvdimp_brbqplots_info" "nvdimp_tradplots_info"
    [4] "hc.syntaxa_groups"     "nvdimp_uplots"         "hc.syntaxa_view"   
    
  3. 繪圖

    syntaxa = rgdal::readOGR(dsn, layer= 'hc.syntaxa_view', verbose = TRUE)
    plot(syntaxa, pch=21, col='black', fill='orange')
    box()
    axis(side=1)
    axis(side=2)
    

References

Comments

comments powered by Disqus