Dear Roger,
Thanks a lot for your very fast answer! Here are the requested outputs:
Le 30/11/2012 10:39, Roger Bivand a ?crit :
On Fri, 30 Nov 2012, Mathieu Basille wrote:
Dear R users,
I have a working PostGIS database that I routinely access through R,
using readOGR or readGDAL. I have this point layer, however, that I can't
read in R:
= "gis.world_ctrdcoast")
OGR data source with driver: PostgreSQL
Source: "PG:dbname=florida host=localhost user=pguser", layer:
"gis.world_ctrdcoast"
with 35398 features and 2 fields
Error in readOGR("PG:dbname=florida host=localhost user=pguser", layer =
"gis.world_ctrdcoast") :
Error getting field 0
The layer is basically a table with just 2 fields: geom (a point feature)
and gid (a serial primary key). I can export this layer without any
problem using ogr2ogr, and the layer looks fine to me (it displays fine
in QGIS):
ogr2ogr -f "ESRI Shapefile" ctrd.shp PG:"dbname=florida host=localhost
user=pguser" "gis.world_ctrdcoast"
Any idea where this is going wrong? Thanks in advance for any suggestion.
Mathieu.
$ ogrinfo "PG:dbname=florida host=localhost user=pguser"
gis.world_ctrdcoast -summary
INFO: Open of `PG:dbname=florida host=localhost user=pguser'
using driver `PostgreSQL' successful.
Layer name: gis.world_ctrdcoast
Geometry: Unknown (any)
Feature Count: 35398
Extent: (-179.916500, -59.083000) - (179.968500, 83.535000)
Layer SRS WKT:
(unknown)
FID Column = gid
Geometry Column = geom
I noticed there was a call to ogrInfo in the readOGR function. After this
line, here is what's reported:
Browse[2]> ogr_info
Source: "PG:dbname=florida host=localhost user=pguser", layer:
"gis.world_ctrdcoast"
Driver: PostgreSQL number of rows 35398
Feature type: wkbPoint with 2 dimensions
NA
Number of fields: 0
[1] name type length typeName
<0 rows> (or 0-length row.names)
Does it make a difference if you add an extra
field?
As a matter of fact, it does! If I add a serial field named 'bla', it works
perfectly:
> ctrdcoast <- readOGR("PG:dbname=florida host=localhost user=pguser",
layer = "gis.world_ctrdcoast")
OGR data source with driver: PostgreSQL
Source: "PG:dbname=florida host=localhost user=pguser", layer:
"gis.world_ctrdcoast"
with 35398 features and 1 fields
Feature type: wkbPoint with 2 dimensions
> summary(ctrdcoast)
Object of class SpatialPointsDataFrame
Coordinates:
min max
coords.x1 -179.9165 179.9685
coords.x2 -59.0830 83.5350
Is projected: NA
proj4string : [NA]
Number of points: 35398
Data attributes:
Min. 1st Qu. Median Mean 3rd Qu. Max.
1 8850 17700 17700 26550 35400
> head(ctrdcoast at data)
bla
1 1
2 2
3 3
4 4
5 5
6 6
I'm wondering whether in this case the GID is not being seen as a
field as such, but the FID of the point. Possibly ogr2ogr steps round this.
It looks like a correct guess, given the output of ogrinfo. In addition,
the exported shapefile from ogr2ogr present one field (named FID), while if
I load the table into QGIS, I can see the field with its correct name
(gid). Admittedly, while I'm pretty sure that QGIS uses GDAL to load
tables, I have no idea how it deals with it.
I don't have a PostGIS system running, so may need help debugging. If you
debug(readOGR), what happens here:
if (drop_unsupported_fields) {
iflds <- as.integer((1:ogr_info$nitems)-1)
iflds <- iflds[keep]
fldnms <- ogr_info$iteminfo$name[keep]
if (any(!keep)) warning(paste("Fields dropped:",
paste(ogr_info$iteminfo$name[!keep], collapse=" ")))
} else {
if (any(!keep)) stop(paste("Unsupported field type:",
paste(ogr_info$iteminfo$typeName[!keep], collapse=" ")))
iflds <- as.integer((1:ogr_info$nitems)-1)
fldnms <- ogr_info$iteminfo$name
}
I think ogr_info$nitems is 1, so we get 0, but if ogr_info$nitems is 0, (0,
-1).
Does that mean that we necessarily need a field in the table (in addition
to geom + gid)? It's not hard to 'fix' by adding a nonsensical field, but I
find it quite counterintuitive.
Thanks again for the help!
Mathieu.
Roger
sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C LC_TIME=fr_FR.UTF-8
[4] LC_COLLATE=fr_FR.UTF-8 LC_MONETARY=fr_FR.UTF-8
LC_MESSAGES=fr_FR.UTF-8
[7] LC_PAPER=C LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=fr_FR.UTF-8
LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] raster_2.0-31 rgdal_0.7-22 sp_1.0-2
loaded via a namespace (and not attached):
[1] compiler_2.15.1 grid_2.15.1 lattice_0.20-6 tools_2.15.1
# SELECT PostGIS_Full_Version();
postgis_full_version
-------------------------------------------------------------------------
POSTGIS="2.1.0SVN r10597" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23
September 2009" GDAL="GDAL 2.0dev, released 2011/12/29" LIBXML="2.8.0"
LIBJSON="UNKNOWN" TOPOLOGY RASTER
~$ whoami
Mathieu Basille, PhD
~$ locate --details
University of Florida \\
Fort Lauderdale Research and Education Center
(+1) 954-577-6314
http://ase-research.org/basille
~$ fortune
? Le tout est de tout dire, et je manque de mots
Et je manque de temps, et je manque d'audace. ?
-- Paul ?luard