[R-sig-Geo] reading PostGIS table into sp data frame

[Edward Vanden Berghe]

Thanks!

Unfortunately I have no luck trying to get readOGR to work directly with the PostGIS database. I tried different variations:

> sp.object <- readOGR(dsn="PG:dbname=...
> sp.object <- readOGR(dsn="my_dsn" ... [where my_dsn is a Windows system dsn]
> sp.object <- readOGR(dsn="PG:my_dsn" ...

all result in the same error message:

	Error in ogrInfo(dsn = dsn, layer = layer, input_field_name_encoding = input_field_name_encoding) : 
	  Cannot open file

Another variation, 

> sp.object <- readOGR(dsn="ODBC:my_dsn" ...

Results in 

	Error in ogrInfo(dsn = dsn, layer = layer, input_field_name_encoding = input_field_name_encoding) : 
	  Multiple # dimensions:

though all he geometries are 2d, and there is a constraint enforcing this in the database. 

It seems rgdal is properly installed: I can use readOGR to read shape files (and this, actually, resolves my problem, as I can use this as a work-around); also, ogrDrivers() returns a list of 43 drivers (but PG, PostgreSQL or PostGIS are not among these; ESRI Shapefile and ODBC are). 

I have no problems connecting to the PostgreSQL database using either DBI/RPostgreSQL, or RODBC (the latter with the same Windows system DSN that I used for the tests above). I am not familiar with GDAL, don't know how to try and connect to the PostgreSQL database from GDAL directly, rather than through rgdal.

Again, any insight in what might go wrong would be most appreciated.

Cheers,

Edward

Message: 4
Date: Wed, 19 Dec 2012 14:58:52 +0100
From: Raffaele Morelli <raffaele.morelli at gmail.com>
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] reading PostGIS table into sp data frame
Message-ID:
	<CAD4guxOt5fSUPbK4VbbPuZf-oTBrcY-UvNPcguvUYt4qS0+_kQ at mail.gmail.com>
Content-Type: text/plain

2012/12/19 Edward Vanden Berghe <evberghe at gmail.com>

> I wanted to create a global map with squares in lat-lon. I have PostGIS
> tables to define these squares  but I havent been able to figure out an
> efficient way of reading those tables into R. The code I am using now is:
>
>        crs <- CRS("+proj=longlat +ellps=WGS84")
>        s <- paste("select id, st_astext(geom) as geom from geo.cs10d";",
> sep="")
>        r <- dbGetQuery(con, s)
>        p <- readWKT(r$geom[1],id=r$id[1],p4s=crs)
>        for(i in 2:length(r$id)){
>               p <- rbind(p, readWKT(r$geom[i], id=r$id[i], p4s=crs))
>        }
>
> where geo.cs10d is the table with squares, id the primary key of the
> table, and geom the binary geometry field.
>
> The code above works fine for the larger squares, such as 10 degrees, of
> which I only need 648 to cover the globe. For finer resolutions, the above
> takes just too long  I assume because the rbind function rewrites the
> whole sp object each time it executes. Ive seen other R scripts that
> initiate an empty data frame of the correct length to go round similar
> problems with the rbind function; I havent been able to find an equivalent
> for spatial polygons. How can I initiate an empty data frame with the right
> structure, and the right length?
>
> A preferable solution would be if there would be a single function to load
> a complete PostGIS table, rather than having to load the polygons one by
> one in a loop. Is there such a function?
>
> Im using PostgreSQL 8.4, PostGIS 1.5, R 2.15.2, platform
> x86_64-w64-mingw32; IDE is StatET 3.0.1 plugin for Eclipse 3.7.2.
>
> Any help would be much appreciated.
>
> Edward
>

Use postgis st_transfom to convert your table in epsg:4326, I would suggest
to use a view for that.
Then use readOGR (in package rgdal) to read the table/view (geom and
attributes) in a sp object with :

sp.object <- readOGR("PG:dbname=your_db host=your_host user=username
password=xxx", "geo.cs10d ")

regards
-r


-- 
*L'unica speranza di catarsi, ammesso che ne esista una, resta affidata
all'istinto di ribellione, alla rivolta non isterilita in progetti, alla
protesta violenta e viscerale. (V. Evangelisti)
*

	[[alternative HTML version deleted]]