[R-sig-Geo] Pointers for querying postGIS and bringing spatial geometries into R

[Barry Rowlingson]

On Tue, Oct 9, 2012 at 4:03 AM, Lyndon Estes <lyndon.estes at gmail.com> wrote:

> library(RPostgreSQL)
> drv <- dbDriver("PostgreSQL")
> con <- dbConnect(drv, dbname = "mydb", user = "me", password = "mypw")
> sql <- "SELECT geom from sa1kgrid where id=9"
> kmlGeom <- dbGetQuery(con, sql)
> (But now I am getting this error:
> In postgresqlExecStatement(conn, statement, ...) :
>   RS-DBI driver warning: (unrecognized PostgreSQL field type geometry
> (id:16504) in column 0))
>

Thats a warning, not an error - if you look at the returned value you
are probably getting an encoded version of the geometry.

If you put a little conversion salt on your SQL you can get something
in text that should be parsable:

 sql
[1] "SELECT ST_AsEWKT(the_geom) from ccsm_polygons limit 2"
> dbGetQuery(con,sql)

                           st_asewkt
1 MULTIPOLYGON(((-179.29688 88.23322,-179.29688 90,-177.89063
90,-177.89063 88.23322,-179.29688 88.23322)))
2       MULTIPOLYGON(((-177.89063 90,-176.48438 90,-176.48438
88.23322,-177.89063 88.23322,-177.89063 90)))

There must be an R WKT reader somewhere...

Or as GML, which you could fiddle about with the XML package in R:

> sql <- "SELECT ST_AsGML(the_geom) from ccsm_polygons limit 2"
> dbGetQuery(con,sql)





                    st_asgml
1 <gml:MultiPolygon><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>-179.296879999999987,88.233220000000003
-179.296879999999987,90 -177.890629999999987,90
-177.890629999999987,88.233220000000003
-179.296879999999987,88.233220000000003</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon>
2
<gml:MultiPolygon><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>-177.890629999999987,90
-176.484379999999987,90 -176.484379999999987,88.233220000000003
-177.890629999999987,88.233220000000003
-177.890629999999987,90</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon>

Or maybe GeoJSON is better:

> sql <- "SELECT ST_AsGeoJSON(the_geom) from ccsm_polygons limit 2"
> dbGetQuery(con,sql)



st_asgeojson
1 {"type":"MultiPolygon","coordinates":[[[[-179.296879999999987,88.233220000000003],[-179.296879999999987,90],[-177.890629999999987,90],[-177.890629999999987,88.233220000000003],[-179.296879999999987,88.233220000000003]]]]}
2
{"type":"MultiPolygon","coordinates":[[[[-177.890629999999987,90],[-176.484379999999987,90],[-176.484379999999987,88.233220000000003],[-177.890629999999987,88.233220000000003],[-177.890629999999987,90]]]]}
>

Either way, you can conceivably get the coordinates out that way. When
I tried the rgdal method:

> library(rgdal)
> mydf <- readOGR("PG:dbname='mydf' user='me' password='mypw'",
> layer="sa1kgrid")
> (I have also tried various permutations of the dbname, quotes, etc., but
> all get me this:
>
> Error in ogrInfo(dsn = dsn, layer = layer, input_field_name_encoding =
> input_field_name_encoding)

  - that worked perfectly for me

 > ogrInfo("PG:dbname='gis'", layer='ccsm_polygons')
Source: "PG:dbname='gis'", layer: "ccsm_polygons"
Driver: PostgreSQL number of rows 32768
Feature type: wkbMultiPolygon with 2 dimensions
NA
Number of fields: 4
        name type length typeName
1       area    0      2  Integer
2  perimeter    0      3  Integer
3 ccsm_polys    0      0  Integer
4 ccsm_pol_1    0      0  Integer

 - so I suspect there's something in your database that readOGR
doesn't like - possibly one of the columns is a type that ogr can't
map to something...


> Oh, one other thing that seems to be possibly relevant, when I try with
> straight up ogr2ogr to access the postGIS database, I get this result:
>
> This is simply from doing this:
>
>> ogr2ogr
> ogr2ogr: error while loading shared libraries: libjpeg.so.8: cannot open
> shared object file: No such file or directory

 From the shell command line? Well that's just broken!

Barry