I'm trying to read in a shapefile using the rgdal library, and am having no luck.
When I try to import using the following syntax:
geo <- readOGR("/path/to/layer","layer")
I'm met with the error
Error in stopifnot(is.list(srl)) : infinite label point
How can I go about diagnosing and fixing the problem? Many thanks.
TL;DR: Likely self-intersecting or overlapping geometries, an issue with rings/holes in the polygon Shapefile.
First off, look at the source for readOGR, by entering the function name without arguments in the prompt. A "find" on the source shows that the code and message in your error are not in that function.
Using this approach, I found the list of functions called by readOGR:
[1] "-" ":" "!" "!=" ".Call"
[6] "(" "[" "[[" "{" "&&"
[11] "%in%" "+" "<" "<-" "<="
[16] "==" ">" "||" "$" "all.equal"
[21] "any" "as.character" "as.integer" "as.logical" "attr"
[26] "attributes" "c" "cat" "cbind" "class"
[31] "comment" "CRS" "data.frame" "do.call" "for"
[36] "function" "gc" "geometry" "getCPLConfigOption" "getGDALVersionInfo"
[41] "iconv" "identical" "if" "ifelse" "integer"
[46] "is.character" "is.na" "is.null" "isTRUE" "lapply"
[51] "length" "Line" "Lines" "list" "make.names"
[56] "match" "match.arg" "max" "message" "missing"
[61] "names" "nchar" "nrow" "ogrFIDs" "ogrInfo"
[66] "paste" "Polygon" "Polygons" "print" "rbind"
[71] "return" "rm" "sapply" "seq" "seq_along"
[76] "setCPLConfigOption" "slot" "sort" "SpatialLines" "SpatialLinesDataFrame"
[81] "SpatialPointsDataFrame" "SpatialPolygons" "SpatialPolygonsDataFrame" "stop" "stopifnot"
[86] "strsplit" "sum" "suppressMessages" "switch" "table"
[91] "try" "unique" "vector" "warning" "which"
There are functions listed there that are actually part of library(sp) (you can tell because they say <environment: namespace:sp> when you print their source), and the error message looks like it comes from the Polygon or Polygons classes.
In Polygons, we find the code bit from the error message: stopifnot(is.list(srl)). Looking at help(Polygons), you see that srl is a "list of Polygon class objects", so the error means that something other than a list is being passed.
So now, we should go back to the source of readOGR and look for calls to Polygons (there are three). The two first call make_Polygonlist (a C rgdal function), the third assembles its own list pllist. That last one is interesting, because the errors are silenced: try(pllist[[j]] <- Polygon(cmat), silent = TRUE), which means this is not likely to be raising your error. That leaves the two others, and brings you to this source code, where we see that makePolygonlist is calling make_Polygon in the same script.
Reading there (and this is not my main language), I see a comment about:
// hole setting based on comments by default (OGC SFS order) // but by ring order if comment NULL 121019
And a bit later:
SET_STRING_ELT(VECTOR_ELT(dimnames, 1), 0, COPY_TO_USER_STRING("x"));
SET_STRING_ELT(VECTOR_ELT(dimnames, 1), 1, COPY_TO_USER_STRING("y"));
which suggests it needs x,y coordinates of some sort in order to handle holes / rings in polygons. Back to help("Polygons-class") we see a parameter labpt which is:
Object of class "numeric"; pair of x, y coordinates giving a label point, the label point of the largest ring component
And now we can give a hypothesis about he cause of the error: it is related to the handling of rings/holes in polygons, and while creating the Polygon objects, it needs label points to correctly handle the rings/holes. However, the label points it finds are not finite, and the error is triggered.
And finally, there is a note:
Polygon objects belonging to an Polygons object should either not overlap one-other, or should be fully included (as lakes or islands in lakes). They should not be self-intersecting.
which suggests that you are likely having issues with a bad geometry, which is overlapping or self-intersecting. You can use tools (for example in QGis) to check and fix the geometry.
