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
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.