3

It is easy to get Vectorlayers into R from a PostgreSQL+PostGIS solution.

library(rgdal)
dsn="PG:dbname=plots host=localhost user=test password=test"
plots = readOGR(dsn,"Plots")

But how to do the same for raster layers stored in-db and/or out-db? I am a bit confused if this is even possible as I haven't found any help or similar questions on the net.

At least readGDAL finds the driver, but returns error messages that I am unsure how to handle.

>dsn="PG:dbname=plots host=localhost user=test password=test port=5432"
>readGDAL(dsn,"map")

PG:dbname=test host=localhost user=test password=test port=5432 has GDAL driver PostGISRaster 
and has 0 rows and 0 columns
Error in offset[2]/(cellsize[1]/abs(gt[2])) : 
  non-numerical arguement in binary operator
additional: Warnings:
1: In dim(x) : no bands in dataset
2: In dim(x) : no bands in dataset

The imported raster should be fine as I can load it via QGIS. If this is possible can someone provide a working example?


OS: Debian Jessie + PostgreSQL 9.4dev + PostGIS 2.. R.version > 3.0.0

GDAL version:

rgdal: version: 0.9-1, (SVN revision 518)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
Path to GDAL shared files: /usr/share/gdal/1.10
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: (autodetected)

GDALInfo (even if pointed directly to table):

>dsn="PG:dbname=plots host=localhost user=test password=test port=5432 table=map"
>rgdal::GDALinfo(dsn)
Error in .local(.Object, ...) : Error retrieving raster metadata
  • what GDAL version? What does gdalinfo say about the datasource? – Mike T Oct 13 '14 at 20:55
  • @MikeT Edited Question. – Curlew Oct 13 '14 at 21:23
  • Outside R, what does gdalinfo from the command line show? What is the schema of map? – Mike T Oct 13 '14 at 22:09
  • gdalinfo "PG:host=localhost port=5432 dbname='plots' user='test' password='test' schema='public' table=map" also refuses to return the raster metadata – Curlew Oct 13 '14 at 22:17
  • Okay, I just tested a plain raster2psql and it works. It has something to do with my custom map table as explained here gis.stackexchange.com/questions/118006/… – Curlew Oct 13 '14 at 22:20
3

So assuming that the raster is correctly loaded in PostGIS you can get your raster via readGDAL in R in the following way:

library(raster)
library(rgdal)

dsn="PG:dbname='plots' host=localhost user='test' password='test' port=5432 schema='gisdata' table='map' mode=2"

ras <- readGDAL(dsn) # Get your file as SpatialGridDataFrame
ras2 <- raster(ras,1) # Convert the first Band to Raster
plot(ras2)

Note the mode-tag at the end which is important depending on how you store your data!

2

If this is still relevant, at the University of Florida, David Bucklin and I have released a rpostgis package that provides bi-directional transfer between PostGIS and R for vector and raster data. The package does not rely on GDAL (and rgdal), and should be platform independent.

rpostgis relies on RPostgreSQL to establish a functional connection between R and the database. You can then use it to import PostGIS raster data type into R using the function pgGetRast, for instance:

library(rpostgis)
conn <- dbConnect(drv = "PostgreSQL", host = "localhost", dbname = "plots", user = "test", password = "test")
my_raster <- pgGetRast(conn, "map")

If the table is not in the public schema (which would be a good idea anyway), you will need to declare it using c("schema", "table") instead. The function assumes that the raster tiles are stored in the column "rast" by default, but you can change that with the argument rast. Now, depending on the size and other considerations, this may be significantly slower (but a lot more flexible) than using readGDAL. We are still working on it, but this is the cost of providing a "pure R" solution. You can also use the boundary argument if you are only interested in a subset of the entire raster (which will significantly increase the loading time).

Note also that there is pgGetGeom for points/lines/polygons, instead of using readOGR.

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.