postgis connection

[fxi]

I'm looking for a way to quickly extract valid spatial data from postgis, and transform geometry if necessary, to display formated result in R leaflet package. For now, i'm using geojson with system("ogr2ogr ...") and -sql argument, and read the output using jsonlite. Minimalist and reliable, but not integrated in R. Other solutions (rgdal/ gdalUtils) did not work for me (always returning a shp..) and readOGR is quite slow and produce a sp object. This seems to be a bottleneck. Could be a direct connection a good improvement/enhancement ? Something like:

geojson_read.postgis(dbCon,query,t_srs)

I'm also experimenting with PGRestAPI, mvt and leaflet plugin Leaflet.MapboxVectorTile. This is incredibly quick and reliable, but for simple query, it's quite convoluted.

Thanks !

[sckott]

hi @fxi Thanks for the message. I haven't used postgis - I'm playing with it now to become more familiar.

What data output do you get when doing a query with postgis? Can you not use the postgres R client https://github.com/rstats-db/RPostgres to query and get data back?

Good timing because we are interested in making it easier to interact with the SQLite spatial client https://github.com/ropensci/spsqlite

[fxi]

Hi, Thanks for your quick reply !

(Sorry for my frenchglish )

First of all, I'm not at all an expert in PostGIS, but I've learned some basic tricks I could share.

You can extract a variety of geometry type from PostGIS: http://postgis.net/docs/reference.html#Geometry_Outputs

If you export geometry as geojson with something like:

SELECT ST_AsGeoJSON(the_geom) from my table;

You get, for each row, a individual geojson feature. So, you have to create the structure yourself, in R or with PostGIS: example with postgis

With Rpostgres, you will have in all case, a data.frame containing the output of the query. There is maybe a way to get raw data, but I didn't find anything on this topic: when you query a DB, you expect a table, that's quite normal.

As I said, ogr2ogr could do the job, but it's seems to be an unnecessary step for me. And you need to write something on your disk. readOGR from rgdal could help, but you loose control of the query. If you are using postgis, there is certainly a reason: big spatial data. So extracting 120MB of geojson doesn't seems to be optimal: the query is not an option.

ogr2ogr attempt:

#' Geojson from postGIS base
#' @param dbInfo Named list with dbName,host,port,user and password
#' @param query PostGIS spatial sql querry.
#' @return geojson list
#' @export
dbGetGeoJSON<-function(dbInfo,query){
  # NOTE: check package geojsonio for topojson and geojson handling.
  # https://github.com/ropensci/geojsonio/issues/61
  dsn <- sprintf("PG:dbname='%s' host='%s' port='%s' user='%s' password='%s'",
    d$dbname,d$host,d$port,d$user,d$password
    )
  tmp <- paste0(tempfile(),".geojson")
  print(tmp)
  system(sprintf("ogr2ogr -f GeoJSON '%s' '%s' -sql '%s' -t_srs '%s'",tmp,dsn,query,"EPSG:4326"))
  return(jsonlite::fromJSON(tmp))
}

usage:

d <- list(host='localhost', dbname='spatial_db', port='5432', user='myusername', password='mypassword')
spatialObj<-dbGetGeoJSON(dbInfo=d,"SELECT * FROM spatial_table")

I've rewritten (see original code) a function that prepare an "unlogged" table (optimised for temp query) and then, use readOGR to get everything back, in SpatialData* format:

#' Transfert postgis feature by sql query to sp object
#' @param dbInfo Named list with dbName,host,port,user and password.
#' @param query PostGIS spatial sql querry.
#' @return spatial object.
#' @export
dbGetSp <- function(dbInfo,query) {
  if(!require('rgdal')|!require(RPostgreSQL))stop('missing rgdal or RPostgreSQL')
  d <- dbInfo
  tmpTbl <- sprintf('tmp_table_%s',round(runif(1)*1e5))
  dsn <- sprintf("PG:dbname='%s' host='%s' port='%s' user='%s' password='%s'",
    d$dbname,d$host,d$port,d$user,d$password
    )
  drv <- dbDriver("PostgreSQL")
  con <- dbConnect(drv, dbname=d$dbname, host=d$host, port=d$port,user=d$user, password=d$password)
  tryCatch({
    sql <- sprintf("CREATE UNLOGGED TABLE %s AS %s",tmpTbl,query)
    res <- dbSendQuery(con,sql)
    nr <- dbGetInfo(res)$rowsAffected
    if(nr<1){
      warning('There is no feature returned.');
      return()
    }
    sql <- sprintf("SELECT f_geometry_column from geometry_columns WHERE f_table_name='%s'",tmpTbl)
    geo <- dbGetQuery(con,sql)
    if(length(geo)>1){
      tname <- sprintf("%s(%s)",tmpTbl,geo$f_geometry_column[1])
    }else{
      tname <- tmpTbl;
    }
    out <- readOGR(dsn,tname)
    return(out)
  },finally={
    sql <- sprintf("DROP TABLE %s",tmpTbl)
    dbSendQuery(con,sql)
    dbClearResult(dbListResults(con)[[1]])
    dbDisconnect(con)
  })
}

Usage:

d <- list(host='localhost', dbname='spatial_db', port='5432', user='myusername', password='mypassword')
spatialObj<-dbGetSp(dbInfo=d,"SELECT * FROM spatial_table")

And I will keep an eye on spsqlite !

Thanks

Fred

[sckott]

Thanks for this, getting to this, just got done with busy hackathon...

[sckott]

hmmm, working on this, but constraint I'm seeing is that we want geojson_read() to always return geojson data (actually, can optionally return spatial class now, but default is geojson as a list) - so if a postgres query a user gives would return points, or individual geojson features (e.g,. {"type":"Point","coordinates":[42.40047,-71.2577]}), we'd want to make sure those are combined into a valid geojson featurecollection

[sckott]

Another problem is that I'd prefer to use the newer Postgres R client https://github.com/rstats-db/RPostgres, but not sure when that will be on CRAN

[sckott]

@fxi see https://github.com/ropensci/geojsonio/blob/master/R/geojson_read.R#L44-L56 and https://github.com/ropensci/geojsonio/blob/master/R/postgis.R

[sckott]

@fxi for your 2nd example above in https://github.com/ropensci/geojsonio/issues/61#issuecomment-109292184 can you share what those tables are: spatial_table and geometry_columns - I'm not sure yet what that example is doing

[fxi]

Thanks !

• I've seen that you have reproduced the exemple from postgisonline. Did you get something interesting ?
• What the advantages of the newer version of RPostgres ?
• The PostGIS geometry_columns table contains geometry info, such as f_geometry_column (name of the column containing the geometry) or coord_dimension (number of spatial dimensions: 2,3 or 4), see see PostGIS user manual, chapter 4.
• spatial_table is your features table name. E.g. dbGetSp(dbInfo=d,"SELECT * FROM locations") return the data in a sp* object. This is not what we want.

[sckott]

I've seen that you have reproduced the exemple from postgisonline. Did you get something interesting ?

I just used that so I could get something tested quickly...

What the advantages of the newer version of RPostgres ?

The newer one will have support moving forward https://github.com/rstats-db/RPostgres

[sckott]

Sorry it's been a while on this. I want to wait until https://github.com/rstats-db/RPostgres is on CRAN, should be hopefully early next year some time.

[fxi]

Thanks a lot for keeping an eye on this.

[fxi]

For future reference : A small tool to convert postgis query to geojson (/topojson) in python. https://gist.github.com/jczaplew/7680118 I did not tried it.

[sckott]

thanks, looks useful

[sckott]

@fxi it has been a long long time, sorry about that. Just pushed up a branch https://github.com/ropensci/geojsonio/tree/geojson_read-postgis with some changes to address your use case. It has been a long time, and now other tools exist in the R world, so maybe this is no longer necessary here.

install remotes::install_github("ropensci/geojsonio@geojson_read-postgis"), see ?geojson_read, an example:

library(geojsonio)
conn <- DBI::dbConnect(RPostgres::Postgres(), dbname = 'postgistest')
state <- "SELECT row_to_json(fc)
 FROM (SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features
 FROM (SELECT 'Feature' As type
    , ST_AsGeoJSON(lg.geog)::json As geometry
    , row_to_json((SELECT l FROM (SELECT loc_id, loc_name) As l
      )) As properties
   FROM locations As lg   ) As f )  As fc;"
json <- geojson_read(conn, query = state, what = "json")
map_leaf(json)

geojson_read on master only returns either a list or an sp package object - I've changed it to return list, sp, or json as character string.

With what = "json" in geojson_read you get back geojson as character, so we can then use map_leaf.character method to make a quick leaflet map

Any thoughts appreciated -- also cc @ateucher

[fxi]

Thanks @sckott !

I don't use it anymore, but I'm curious. It could still be useful.

I tried to install it, but I've failed. I've got many dependencies issues : v8, rgdal, geos, rcpp, mime, curl, jqr .. As soon as I've more time, I will try again.

Do you have a docker image pre-packaged with RPostgres ?

Cheers, F

[sckott]

thanks for the quick reply @fxi - I'll see if I can whip up a docker image

[fxi]

I've managed to build one that should do it

FROM rocker/geospatial

RUN apt-get update \
  && apt-get install -y --no-install-recommends \
    libv8-dev \
    libpq-dev \
    postgresql \
  && install2.r --error \
   RPostgres \
   remotes \
    && R -e "remotes::install_github('ropensci/geojsonio@geojson_read-postgis')"

I will try to test it soon...

[sckott]

i made an attempt here https://github.com/ropensci/geojsonio/commit/f8ad43c7150eaed89605d66b85bf0f930eb9d295 - doesn't have postgres itself

[github-actions[bot]]

This issue has been automatically locked. If you believe you have found a related problem, please file a new issue (with a reprex: https://reprex.tidyverse.org) and link to this issue.