reading PostGIS table into sp data frame
Sorry, me again. The PostGIS driver might be missing from the pre-compiled version of rgdal for Windows, but the ODBC driver isn't; so that should be an alternative to using a driver specific for PostGIS. But as already mentioned, the following R statement
p <- readOGR(dsn="ODBC:PostgreSQL VM", "geo.cs10d")
results in the following output: Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding) : Multiple # dimensions: I checked the log files of the postgresql server; there seems to be an sql statement asking for a non-existing field: 192.168.56.101(49436): evberghe: 2012-12-22 12:27:00 EST: LOG: connection authorized: user=evberghe database=obis 192.168.56.101(49436): evberghe: 2012-12-22 12:27:00 EST: ERROR: column "geometry_type" does not exist at character 41 192.168.56.101(49436): evberghe: 2012-12-22 12:27:00 EST: STATEMENT: SELECT f_table_name, f_geometry_column, geometry_type FROM geometry_columns I assume rgdal is checking whether all records have the same number of dimensions and the same geometry type. The geometry type is in a field 'type', not in a field 'geometry_type'. Is it because the PostGIS driver is not available that the ODBC driver doesn't know the correct field in which to look for the geometry type? I could ask the PostgreSQL server to log everything and post the relevant sections of the log if this would help. Edward -----Original Message----- From: Roger Bivand [mailto:Roger.Bivand at nhh.no] Sent: 21 December 2012 13:31 To: Edward Vanden Berghe Cc: r-sig-geo at r-project.org Subject: RE: [R-sig-Geo] reading PostGIS table into sp data frame
On Fri, 21 Dec 2012, Edward Vanden Berghe wrote:
The GDAL was installed from www.gisinternals.com/sdk; I installed the 32 bits version, though my Windows version is 64: the table on gisinternals.com stated that there was a build problem with the 64-bits version. rgdal was installed from CRAN, as a binary, version 0.7-20, build 2.15.1. R version is 2.15.2 (2012-10-26) -- "Trick or Treat", platform x86_64-w64-mingw32/x64 (64-bit). The output from ogrDrivers() is pasted below. If I understand you correctly, and since I am using the standard CRAN Windows binary package, I might be out of luck?
Not out of luck, the users of the CRAN Windows and OSX rgdal binaries are very lucky indeed, thanks to hard work by Uwe Ligges, Brian Ripley and Simon Urbanek. At least many of the most commonly used drivers are present.
There never is a complete set of drivers ionstalled in any GDAL, some depend on proprietary dynamically linked objects.
file.show(system.file("README.windows", package="rgdal"))
displays the accumulated - but *not* updated - wisdom with regard to installing source rgdal on Windows, linked to an external GDAL DLL. If you try to follow this route, please pass back your experience with regard to installing - it is likely that much has changed. You may find that OSGeo4W is another source for a Win32 GDAL binary, but I don't see from the documentation whether it supports PostGIS.
Installing Linux, PostgreSQL, PostGIS, GDAL, Proj.4, etc., then R and rgdal under Linux will probably require less jumping through hoops than building rgdal against an external 64-bit GDAL under windows.
Hope this helps,
Roger
Cheers,
Edward
Output from ogrDrivers():
name write
1 AeronavFAA FALSE
2 ARCGEN FALSE
3 AVCBin FALSE
4 AVCE00 FALSE
5 BNA TRUE
6 CSV TRUE
7 DGN TRUE
8 DXF TRUE
9 EDIGEO FALSE
10 ESRI Shapefile TRUE
11 Geoconcept TRUE
12 GeoJSON TRUE
13 Geomedia FALSE
14 GeoRSS TRUE
15 GML TRUE
16 GMT TRUE
17 GPSBabel TRUE
18 GPSTrackMaker TRUE
19 GPX TRUE
20 HTF FALSE
21 Idrisi FALSE
22 KML TRUE
23 MapInfo File TRUE
24 Memory TRUE
25 MSSQLSpatial TRUE
26 ODBC TRUE
27 OpenAir FALSE
28 PCIDSK TRUE
29 PDS FALSE
30 PGDump TRUE
31 PGeo FALSE
32 REC FALSE
33 S57 TRUE
34 SDTS FALSE
35 SEGUKOOA FALSE
36 SEGY FALSE
37 SUA FALSE
38 SVG FALSE
39 TIGER TRUE
40 UK .NTF FALSE
41 VFK FALSE
42 VRT FALSE
43 XPlane FALSE
-----Original Message-----
From: Roger Bivand [mailto:Roger.Bivand at nhh.no]
Sent: 20 December 2012 20:58
To: Edward Vanden Berghe
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] reading PostGIS table into sp data frame
On Thu, 20 Dec 2012, Edward Vanden Berghe wrote:
Thanks! Unfortunately I have no luck trying to get readOGR to work directly with the PostGIS database. I tried different variations:
Did you check whether your rgdal/GDAL supports the PostGIS driver? Did you look at the output of ogrDrivers()? If you did not install rgdal as a source package built against a GDAL that has the PostGIS driver on your platform (maybe OSGeo4W??), but are using the standard CRAN Windows binary package, then it is simply the case that the driver is not supported. On OSX using the Kyngchaos frameworks, or on Linux, this would be simple, but it isn't simple on Windows. You need to establish how rgdal was installed and report that. Hope this clarifies, Roger
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
-- Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no
-- Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no