3

I would like to be able to take my SQL queries for PostGIS and make each row for the geom attribute part of an sp object so I can plot it with R's leaflet bindings. I can't really share what my database but I will share how I learned how to connect and what I am trying to do to bind each row with a geom into a single sp object. When I use readWKT for a single row the object class is SpatialPolygons.

#make my query to return multiple rows with just the geom attribute
query <- "SELECT ST_AsText(geom) FROM assembly"

#connect to db, send query, recieve table, and clear the result
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, host = "localhost", port = "5432", 
                 dbname = "anc", user = "postgres")
res <- dbSendQuery(con, query)
df <- dbFetch(res) ; dbClearResult(res)

Here is where I give it a shot to bring the rows together into a list object similar to what I've found with this answer. 

#preallocate vector for combining spatial objects into a list
n <- dim(df)[1]
map_obj <- vector(mode = "list", length = n)

for(i in seq(n)) {
  map_obj[i] <- readWKT(df[i,1])
}

#Put the spatial objects together? 
Polygons(map_obj,1:n)

An output of str(readWKT(df[1,1])):

Formal class 'SpatialPolygons' [package "sp"] with 4 slots
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 1684530 2627594
  .. .. .. .. .. .. ..@ area   : num 3.54e+08
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:820, 1:2] 1690044 1690045 1690046 1690046 1690047 ...
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 1684530 2627594
  .. .. .. ..@ ID       : chr "1"
  .. .. .. ..@ area     : num 3.54e+08
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] 1670993 2612705 1698389 2642795
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA
Community
  • 1
cylondude
  • 336
  • 2
  • 11

1 Answers1

1

You can do this with a set of *apply functions:

df<-dbGetQuery(con,'select st_astext(geom) from schema.table where geom is not null;')

map_obj<-mapply(function(x) readWKT(x), x=df[,1])

Spol <- SpatialPolygons(lapply(1:length(map_obj), function(i) {
    pol <- slot(map_obj[[i]], "polygons")[[1]]
    slot(pol, "ID") <- as.character(i)
    pol
}))

If you wanted to preserve a unique ID from the database table, you could modify it as follows:

df<-dbGetQuery(con,'select gid, st_astext(geom) from schema.table where geom is not null;')

map_obj<-mapply(function(x,y) readWKT(x,y), x=df[,2], y=df[,1])

Spol <- SpatialPolygons(lapply(1:length(map_obj), function(i) {
    pol <- slot(map_obj[[i]], "polygons")[[1]]
    slot(pol, "ID") <- slot(slot(map_obj[[i]], "polygons")[[1]],"ID") ##assign original ID to polygon
    pol
}))

I found myself needing to do this a lot, so I wrote a package in R to wrap these functions, along with handling the projection and data from the table (to create Spatial*DataFrames). You can install it from Github with a function from the package devtools:

library(devtools)
install_github('dnbucklin/pgis2r')

The pgis2spol() function could be used after you establish your connection to the database, for example:

library(pgis2r)
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, host = "localhost", port = "5432", dbname = "anc", user = "postgres")
Spol_df<-pgis2spol(con,'schema.tablename',geom='geom',gid='unique_id_column',proj=TRUE)
David
  • 21
  • 2