使用 R 在 ESRI Shapefile 中读取错误答案

【问题标题】：Error reading in ESRI Shapefile using R使用 R 在 ESRI Shapefile 中读取错误
【发布时间】：2015-03-10 05:54:15
【问题描述】：

我正在尝试使用 rgdal 库读取 shapefile，但运气不佳。

当我尝试使用以下语法导入时：

geo <- readOGR("/path/to/layer","layer")

我遇到了错误

Error in stopifnot(is.list(srl)) : infinite label point

如何诊断和解决问题？非常感谢。

【问题讨论】：

这个问题应该在这里问：gis.stackexchange.com.
你能在命令行上运行ogrinfo -so /path/to/layer layer吗（在Linux/Mac上最简单）？ shapefile 可能以某种方式具有一些无限坐标...检查该命令的范围输出。
你在这里和那里交叉发布。请删除其中一篇帖子。
接缝合理的范围：(-179.147340, 17.881242) - (179.778465, 71.390482)

标签： r gdal rgdal

【解决方案1】：

TL;DR：可能是自相交或重叠的几何图形，多边形 Shapefile 中的环/孔存在问题。

首先，通过在提示符中输入不带参数的函数名来查看readOGR 的来源。对源代码的“查找”表明错误中的代码和消息不在该函数中。

使用this approach，我找到了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"

那里列出的函数实际上是 library(sp) 的一部分（你可以知道，因为当你打印它们的源代码时它们会说 <environment: namespace:sp>），错误消息看起来像 comes from Polygon 或 @987654330 @类。

在Polygons 中，我们从错误消息中找到代码位：stopifnot(is.list(srl))。查看help(Polygons)，您会看到srl 是“Polygon 类对象的列表”，因此错误意味着传递的不是列表。

所以现在，我们应该回到readOGR 的来源，寻找对Polygons 的调用（共有三个）。两个第一个调用make_Polygonlist（一个C rgdal 函数），第三个组装它自己的列表pllist。最后一个很有趣，因为错误被静音：try(pllist[[j]] <- Polygon(cmat), silent = TRUE)，这意味着这不太可能引发您的错误。剩下两个，然后将您带到this source code，我们看到makePolygonlist 在同一个脚本中调用make_Polygon。

在那里阅读（这不是我的主要语言），我看到以下评论：

//默认基于cmets的孔设置（OGC SFS顺序） // 但如果注释为 NULL 121019，则按响铃顺序

稍后：

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"));

这表明它需要某种 x,y 坐标才能处理多边形中的孔/环。回到help("Polygons-class")我们看到一个参数labpt是：

“数字”类的对象；一对 x, y 坐标给出一个标签 point，最大环组件的标签点

现在我们可以对错误的原因给出一个假设：它与多边形中环/孔的处理有关，并且在创建Polygon 对象时，它需要标签点来正确处理环/孔.但是它找到的标签点不是有限的，会触发错误。

最后，有一个注释：

属于 Polygons 对象的 Polygon 对象不应该相互重叠，或应完全包括在内（如湖泊或岛屿湖泊）。它们不应该是自相交的。

这表明您可能遇到了重叠或自相交的不良几何形状的问题。您可以使用工具（例如在 QGis 中）来检查和修复几何图形。

【讨论】：