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

Fri Oct 19 17:02:28 CEST 2012

Hi Edzer,

Thanks for the tip, and sorry for the slow response.

I tried this, based on the approach from the stpg vignette.  It looks
like it is getting closer, but still not quite successful.

Here's my approach (sorry I can't give a toy postgis example).  This
is trying to load in a set of 609 polygons.

# Libraries
library(RPostgreSQL)
library(rgdal)

dbname <- "mydb"
user <- "uname"
password <- "mypw"
OGRstring <- paste("PG:dbname=", dbname, " user=", user, " password=",
password, sep = "")
sitepolys <- readOGR(OGRstring, layer="sites", driver = "PostgreSQL")

> sitepolys <- readOGR(OGRstring, layer="newqaqc_sites")
OGR data source with driver: PostgreSQL
Source: "PG:dbname=mydb user=uname password=mypw", layer: "sites"
with 1 features and 2 fields
Feature type: NA with 2 dimensions
Error in readOGR(OGRstring, layer = "sites") :
  Incompatible geometry: 1853189987
In addition: Warning message:
In readOGR(OGRstring, layer = "sites") :
  Deleted feature IDs: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47,
48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112,
113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126,
127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140,
141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154,
155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168,
169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182,
183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196,
197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210,
211, 212, 213, 214, 215, 216, 217, [... truncated]

So, it looks like it could be working, but there is something
unappealing about geometries.

Still, it works fine when I load the same data from a shapefile:

qaqc.polys <- readOGR(dsn = "sites.shp", layer = "sites")

So maybe it got broken when being brought into PostGIS.

If I get readOGR working with password and user name from postGIS, I
will post back.  Otherwise, thanks again for the suggestions.

Best, Lyndon

On Wed, Oct 10, 2012 at 9:30 AM, Edzer Pebesma
<edzer.pebesma at uni-muenster.de> wrote:
> Lyndon,
>
> as a side note, readOGR with the postGIS driver does work nicely, as
> shown e.g. here:
>
> http://wiki.intamap.org/index.php/PostGIS
>
> (where no username / passwd was set); more recent (but deeper hidden) I
> reported attempts with username / passwd here:
> http://cran.r-project.org/web/packages/spacetime/vignettes/stpg.pdf (in
> this example, readOGR was not used because only the attribute tables
> resulting from spatial and/or temporal selections in postGIS were retrieved)
>
> On 10/10/2012 02:01 PM, Lyndon Estes wrote:
>> Hi to everyone,
>>
>> Thanks to all of you for the advice on this, which has solved my problem
>> quite nicely now.  Based on your suggestions, I figured out the following:
>>
>> ogr2ogr was broken for me because I had not completely removed a prior
>> custom install of gdal that I had set up for hdf4 support, which required
>> libjpeg. Once removed, I can use ogr2ogr fine from the linux shell.
>>
>> Since I want to select a particular record from the postgis database, and
>> want to keep the object in R without writing it to disk for the time being
>> (which would rule out a system call to ogr2ogr as an option), I am using
>> the RpostgreSQL option together with rgeos::readWKT.  I am basing this on
>> the idea that it is not possible to use sql statements in readOGR, and thus
>> can't select an individual geometry from a table (although I see there was
>> some talk about this in 2010:
>> http://r-sig-geo.2731867.n2.nabble.com/passing-SQL-through-readOGR-td3816889.html).
>>
>>
>> library(RPostgreSQL)
>> # CRS
>> drv <- dbDriver("PostgreSQL")
>> con <- dbConnect(drv, dbname = "mydb", user = "me", password = "mypw")
>> sql <- "SELECT ST_AsEWKT(geom) from sa1kgrid where id=9"
>> kmlGeom <- dbGetQuery(con, sql)
>> kmlGeom
>> 1 SRID=97490;MULTIPOLYGON(((570779.01553729
>> -2401338.07340772,571779.01553729 -2401338.07340772,571779.01553729
>> -2402338.07340772,570779.01553729 -2402338.07340772,570779.01553729
>> -2401338.07340772)))
>> # Comes back as data.frame with SRID, which needs to be stripped out for
>> readWKT to work
>> geomStr <- unlist(strsplit(kmlGeom[1, ], ";"))
>> geomPoly <- readWKT(geomStr[-grep("SRID", geomStr)])
>>
>> That seems to work nicely, so thanks for the solution!
>>
>> Best, Lyndon
>>
>>
>>
>>
>>
>> On Tue, Oct 9, 2012 at 5:21 PM, Edzer Pebesma <edzer.pebesma at uni-muenster.de
>>> wrote:
>>
>>>
>>>
>>> On 10/09/2012 12:04 PM, Barry Rowlingson wrote:
>>>> There must be an R WKT reader somewhere...
>>>
>>> rgeos::readWKT
>>> --
>>> Edzer Pebesma
>>> Institute for Geoinformatics (ifgi), University of Münster
>>> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
>>> 8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
>>> http://www.52north.org/geostatistics      e.pebesma at wwu.de
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of Münster
> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
> 8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
> http://www.52north.org/geostatistics      e.pebesma at wwu.de
>