logo标签列表

从SpatialPolygonsDataFrame中提取坐标

rgeospatialslot
11

有没有只有我在从SpatialPolygonsDataFrame对象中提取多边形的坐标时遇到问题？我能够提取对象的其他插槽（ID、plotOrder），但不能提取坐标（coords）。我不知道我做错了什么。请查看下面的R会话，其中bdryData是具有两个多边形的SpatialPolygonsDataFrame对象。

> bdryData
An object of class "SpatialPolygonsDataFrame"
Slot "data":
  ID GRIDCODE
0  1        0
1  2        0

Slot "polygons":
[[1]]
An object of class "Polygons"
Slot "Polygons":
[[1]]
An object of class "Polygon"
Slot "labpt":
[1] 415499.1 432781.7

Slot "area":
[1] 0.6846572

Slot "hole":
[1] FALSE

Slot "ringDir":
[1] 1

Slot "coords":
         [,1]     [,2]
[1,] 415499.6 432781.2
[2,] 415498.4 432781.5
[3,] 415499.3 432782.4
[4,] 415499.6 432781.2



Slot "plotOrder":
[1] 1

Slot "labpt":
[1] 415499.1 432781.7

Slot "ID":
[1] "0"

Slot "area":
[1] 0.6846572


[[2]]
An object of class "Polygons"
Slot "Polygons":
[[1]]
An object of class "Polygon"
Slot "labpt":
[1] 415587.3 432779.4

Slot "area":
[1] 20712.98

Slot "hole":
[1] FALSE

Slot "ringDir":
[1] 1

Slot "coords":
           [,1]     [,2]
  [1,] 415499.6 432781.2
  [2,] 415505.0 432781.8
  [3,] 415506.5 432792.6
  [4,] 415508.9 432792.8
  [5,] 415515.0 432791.5
  [6,] 415517.7 432795.6
  [7,] 415528.6 432797.7
  [8,] 415538.8 432804.2
  [9,] 415543.2 432805.8
 [10,] 415545.1 432803.6
 [11,] 415547.1 432804.7
 [12,] 415551.7 432805.8
 [13,] 415557.5 432812.3
 [14,] 415564.2 432817.1
 [15,] 415568.5 432823.9
 [16,] 415571.0 432826.8
 [17,] 415573.2 432828.7
 [18,] 415574.1 432829.7
 [19,] 415576.2 432830.7
 [20,] 415580.2 432833.8
 [21,] 415589.6 432836.0
 [22,] 415593.1 432841.0
 [23,] 415592.2 432843.7
 [24,] 415590.6 432846.6
 [25,] 415589.0 432853.3
 [26,] 415584.8 432855.3
 [27,] 415579.7 432859.8
 [28,] 415577.7 432866.2
 [29,] 415575.6 432868.1
 [30,] 415566.7 432880.7
 [31,] 415562.7 432887.5
 [32,] 415559.2 432889.1
 [33,] 415561.5 432890.7
 [34,] 415586.2 432889.7
 [35,] 415587.1 432888.6
 [36,] 415588.5 432890.2
 [37,] 415598.2 432888.7
 [38,] 415599.1 432887.7
 [39,] 415601.2 432886.7
 [40,] 415603.1 432885.7
 [41,] 415605.2 432884.7
 [42,] 415606.1 432882.7
 [43,] 415607.2 432880.7
 [44,] 415608.3 432878.3
 [45,] 415612.2 432874.8
 [46,] 415614.7 432871.9
 [47,] 415617.1 432870.7
 [48,] 415622.4 432868.2
 [49,] 415622.0 432862.4
 [50,] 415624.2 432855.4
 [51,] 415633.2 432845.3
 [52,] 415639.0 432841.1
 [53,] 415642.8 432832.9
 [54,] 415647.5 432828.7
 [55,] 415654.3 432820.3
 [56,] 415654.1 432816.5
 [57,] 415658.2 432812.8
 [58,] 415661.9 432808.6
 [59,] 415663.5 432808.7
 [60,] 415668.1 432803.5
 [61,] 415676.5 432801.3
 [62,] 415679.1 432802.7
 [63,] 415680.1 432802.7
 [64,] 415681.1 432802.7
 [65,] 415682.2 432802.7
 [66,] 415685.8 432804.7
 [67,] 415691.8 432802.2
 [68,] 415693.6 432798.9
 [69,] 415696.2 432777.0
 [70,] 415689.8 432773.5
 [71,] 415683.7 432771.6
 [72,] 415680.2 432766.7
 [73,] 415679.0 432765.6
 [74,] 415676.8 432753.7
 [75,] 415671.4 432747.7
 [76,] 415662.7 432747.2
 [77,] 415658.7 432750.0
 [78,] 415657.0 432746.3
 [79,] 415654.1 432743.7
 [80,] 415652.3 432739.8
 [81,] 415649.6 432739.6
 [82,] 415648.0 432739.7
 [83,] 415641.9 432736.4
 [84,] 415633.4 432736.9
 [85,] 415630.2 432734.7
 [86,] 415622.3 432733.6
 [87,] 415614.4 432726.5
 [88,] 415617.1 432719.1
 [89,] 415612.5 432718.1
 [90,] 415610.0 432720.9
 [91,] 415606.2 432716.6
 [92,] 415603.2 432713.9
 [93,] 415601.4 432710.0
 [94,] 415580.3 432708.7
 [95,] 415545.1 432709.7
 [96,] 415543.5 432711.5
 [97,] 415534.0 432715.7
 [98,] 415527.1 432713.7
 [99,] 415521.1 432711.6
[100,] 415505.6 432710.6
[101,] 415501.3 432710.9
[102,] 415499.3 432708.7
[103,] 415495.6 432711.6
[104,] 415482.6 432726.2
[105,] 415477.2 432734.0
[106,] 415478.1 432737.7
[107,] 415479.2 432739.7
[108,] 415480.9 432743.4
[109,] 415486.5 432751.2
[110,] 415493.2 432760.7
[111,] 415494.1 432762.7
[112,] 415498.1 432767.9
[113,] 415497.2 432770.7
[114,] 415490.6 432773.2
[115,] 415493.2 432775.6
[116,] 415496.0 432778.7
[117,] 415499.2 432779.7
[118,] 415499.6 432781.2



Slot "plotOrder":
[1] 1

Slot "labpt":
[1] 415587.3 432779.4

Slot "ID":
[1] "1"

Slot "area":
[1] 20712.98



Slot "plotOrder":
[1] 2 1

Slot "bbox":
       min      max
x 415477.2 415696.2
y 432708.7 432890.7

Slot "proj4string":
CRS arguments:
 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000
+datum=OSGB36 +units=m +no_defs +ellps=airy
+towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894

从bdryData中选择第二个多边形

> bdryData@polygons[[2]]
An object of class "Polygons"
Slot "Polygons":
[[1]]
An object of class "Polygon"
Slot "labpt":
[1] 415587.3 432779.4

Slot "area":
[1] 20712.98

Slot "hole":
[1] FALSE

Slot "ringDir":
[1] 1

Slot "coords":
           [,1]     [,2]
  [1,] 415499.6 432781.2
  [2,] 415505.0 432781.8
  [3,] 415506.5 432792.6
  [4,] 415508.9 432792.8
  [5,] 415515.0 432791.5
  [6,] 415517.7 432795.6
  [7,] 415528.6 432797.7
  [8,] 415538.8 432804.2
  [9,] 415543.2 432805.8
 [10,] 415545.1 432803.6
 [11,] 415547.1 432804.7
 [12,] 415551.7 432805.8
 [13,] 415557.5 432812.3
 [14,] 415564.2 432817.1
 [15,] 415568.5 432823.9
 [16,] 415571.0 432826.8
 [17,] 415573.2 432828.7
 [18,] 415574.1 432829.7
 [19,] 415576.2 432830.7
 [20,] 415580.2 432833.8
 [21,] 415589.6 432836.0
 [22,] 415593.1 432841.0
 [23,] 415592.2 432843.7
 [24,] 415590.6 432846.6
 [25,] 415589.0 432853.3
 [26,] 415584.8 432855.3
 [27,] 415579.7 432859.8
 [28,] 415577.7 432866.2
 [29,] 415575.6 432868.1
 [30,] 415566.7 432880.7
 [31,] 415562.7 432887.5
 [32,] 415559.2 432889.1
 [33,] 415561.5 432890.7
 [34,] 415586.2 432889.7
 [35,] 415587.1 432888.6
 [36,] 415588.5 432890.2
 [37,] 415598.2 432888.7
 [38,] 415599.1 432887.7
 [39,] 415601.2 432886.7
 [40,] 415603.1 432885.7
 [41,] 415605.2 432884.7
 [42,] 415606.1 432882.7
 [43,] 415607.2 432880.7
 [44,] 415608.3 432878.3
 [45,] 415612.2 432874.8
 [46,] 415614.7 432871.9
 [47,] 415617.1 432870.7
 [48,] 415622.4 432868.2
 [49,] 415622.0 432862.4
 [50,] 415624.2 432855.4
 [51,] 415633.2 432845.3
 [52,] 415639.0 432841.1
 [53,] 415642.8 432832.9
 [54,] 415647.5 432828.7
 [55,] 415654.3 432820.3
 [56,] 415654.1 432816.5
 [57,] 415658.2 432812.8
 [58,] 415661.9 432808.6
 [59,] 415663.5 432808.7
 [60,] 415668.1 432803.5
 [61,] 415676.5 432801.3
 [62,] 415679.1 432802.7
 [63,] 415680.1 432802.7
 [64,] 415681.1 432802.7
 [65,] 415682.2 432802.7
 [66,] 415685.8 432804.7
 [67,] 415691.8 432802.2
 [68,] 415693.6 432798.9
 [69,] 415696.2 432777.0
 [70,] 415689.8 432773.5
 [71,] 415683.7 432771.6
 [72,] 415680.2 432766.7
 [73,] 415679.0 432765.6
 [74,] 415676.8 432753.7
 [75,] 415671.4 432747.7
 [76,] 415662.7 432747.2
 [77,] 415658.7 432750.0
 [78,] 415657.0 432746.3
 [79,] 415654.1 432743.7
 [80,] 415652.3 432739.8
 [81,] 415649.6 432739.6
 [82,] 415648.0 432739.7
 [83,] 415641.9 432736.4
 [84,] 415633.4 432736.9
 [85,] 415630.2 432734.7
 [86,] 415622.3 432733.6
 [87,] 415614.4 432726.5
 [88,] 415617.1 432719.1
 [89,] 415612.5 432718.1
 [90,] 415610.0 432720.9
 [91,] 415606.2 432716.6
 [92,] 415603.2 432713.9
 [93,] 415601.4 432710.0
 [94,] 415580.3 432708.7
 [95,] 415545.1 432709.7
 [96,] 415543.5 432711.5
 [97,] 415534.0 432715.7
 [98,] 415527.1 432713.7
 [99,] 415521.1 432711.6
[100,] 415505.6 432710.6
[101,] 415501.3 432710.9
[102,] 415499.3 432708.7
[103,] 415495.6 432711.6
[104,] 415482.6 432726.2
[105,] 415477.2 432734.0
[106,] 415478.1 432737.7
[107,] 415479.2 432739.7
[108,] 415480.9 432743.4
[109,] 415486.5 432751.2
[110,] 415493.2 432760.7
[111,] 415494.1 432762.7
[112,] 415498.1 432767.9
[113,] 415497.2 432770.7
[114,] 415490.6 432773.2
[115,] 415493.2 432775.6
[116,] 415496.0 432778.7
[117,] 415499.2 432779.7
[118,] 415499.6 432781.2



Slot "plotOrder":
[1] 1

Slot "labpt":
[1] 415587.3 432779.4

Slot "ID":
[1] "1"

Slot "area":
[1] 20712.98
提取插槽
> bdryData@polygons[[2]]@ID 
[1] "1"

> bdryData@polygons[[2]]@plotOrder
[1] 1

但是坐标存在问题。

> bdryData@polygons[[2]]@coords
Error: no slot of name "coords" for this object of class "Polygons"

非常感谢您的帮助。

谢谢。

6个回答
15

最后，我弄清楚了我没有正确解析输出的问题。 正确的方法是 bdryData@polygons[[2]]@Polygons[[1]]@coords。 注意命令中的差异polygons(Polygons和polygons)，而且花费了我很长时间才找到。

1
这段代码是为哪个多边形设计的？我无法理解[[2]]和[[1]]引用不匹配的原因。 - LoveMeow
2
这不是可推广的，就像@Ms.Meow一样，我不确定2从哪里来。 - forlooper
2
这对我来说完美地工作了。我的列表中只有一个多边形，因此我用[[1]]替换了[[2]]，用我的SpatialPolygons对象替换了bdryData，然后就能提取坐标了。 - colonelforbin97
6
使用sp包中的coordinates()函数。它应该以列表格式给出值。
你也可以从shapefile中获取多边形属性。
mfile = readOGR(dsn=dsn,layer=layername)
polys = attr(mfile,'polygons')
npolys = length(polys)
for (i in 1:npolys){
  poly = polys[[i]]
  polys2 = attr(poly,'Polygons')
  npolys2 = length(polys2)
  for (j in 1:npolys2){
     #do stuff with these values
     coords = coordinates(polys2[[j]])
  }
}
谢谢。但是真正的问题在于我没有正确解析它，花了我很长时间才解决。请检查我的答案。 - rm167
1
coordinates() 函数检索的是 labpt 插槽，而不是多边形坐标。 - cmbarbu
@cmbarbu所说的是正确的，coordinates()不起作用。 - forlooper
5
这对我来说也花了一些时间才想明白。下面是我写的函数，适用于我的情况。sp.df 应该是一个 SpatialPolygonsDataFrame。
extractCoords <- function(sp.df)
{
    results <- list()
    for(i in 1:length(sp.df@polygons[[1]]@Polygons))
    {
        results[[i]] <- sp.df@polygons[[1]]@Polygons[[i]]@coords
    }
    results <- Reduce(rbind, results)
    results
}
2
这个问题也在gis.stackexchange上得到了解答，在这里。我做了一个例子，测试了@mdsumner提到的所有选项在这里。还可以看看在这里。
library(sp)
library(sf)
#> Warning: package 'sf' was built under R version 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(raster)
library(spbabel)
#> Warning: package 'spbabel' was built under R version 3.5.3
library(tmap)

library(microbenchmark)
library(ggplot2)

# Prepare data
data(World)
# Convert from sf to sp objects
atf_sf <- World[World$iso_a3 == "ATF", ]
atf_sp <- as(atf_sf, "Spatial")
atf_sp
#> class       : SpatialPolygonsDataFrame 
#> features    : 1 
#> extent      : 5490427, 5660887, -6048972, -5932855  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
#> variables   : 15
#> # A tibble: 1 x 15
#>   iso_a3 name  sovereignt continent area  pop_est pop_est_dens economy
#>   <fct>  <fct> <fct>      <fct>     <S3:>   <dbl>        <dbl> <fct>  
#> 1 ATF    Fr. ~ France     Seven se~ 7257~     140       0.0193 6. Dev~
#> # ... with 7 more variables: income_grp <fct>, gdp_cap_est <dbl>,
#> #   life_exp <dbl>, well_being <dbl>, footprint <dbl>, inequality <dbl>,
#> #   HPI <dbl>

# Try various functions:

raster::geom(atf_sp)
#>       object part cump hole       x        y
#>  [1,]      1    1    1    0 5550200 -5932855
#>  [2,]      1    1    1    0 5589907 -5964836
#>  [3,]      1    1    1    0 5660887 -5977490
#>  [4,]      1    1    1    0 5656160 -5996685
#>  [5,]      1    1    1    0 5615621 -6042456
#>  [6,]      1    1    1    0 5490427 -6048972
#>  [7,]      1    1    1    0 5509148 -5995424
#>  [8,]      1    1    1    0 5536900 -5953683
#>  [9,]      1    1    1    0 5550200 -5932855

ggplot2::fortify(atf_sp)
#> Regions defined for each Polygons
#>      long      lat order  hole piece id group
#> 1 5550200 -5932855     1 FALSE     1  8   8.1
#> 2 5589907 -5964836     2 FALSE     1  8   8.1
#> 3 5660887 -5977490     3 FALSE     1  8   8.1
#> 4 5656160 -5996685     4 FALSE     1  8   8.1
#> 5 5615621 -6042456     5 FALSE     1  8   8.1
#> 6 5490427 -6048972     6 FALSE     1  8   8.1
#> 7 5509148 -5995424     7 FALSE     1  8   8.1
#> 8 5536900 -5953683     8 FALSE     1  8   8.1
#> 9 5550200 -5932855     9 FALSE     1  8   8.1

spbabel::sptable(atf_sp)
#> # A tibble: 9 x 6
#>   object_ branch_ island_ order_       x_        y_
#>     <int>   <int> <lgl>    <int>    <dbl>     <dbl>
#> 1       1       1 TRUE         1 5550200. -5932855.
#> 2       1       1 TRUE         2 5589907. -5964836.
#> 3       1       1 TRUE         3 5660887. -5977490.
#> 4       1       1 TRUE         4 5656160. -5996685.
#> 5       1       1 TRUE         5 5615621. -6042456.
#> 6       1       1 TRUE         6 5490427. -6048972.
#> 7       1       1 TRUE         7 5509148. -5995424.
#> 8       1       1 TRUE         8 5536900. -5953683.
#> 9       1       1 TRUE         9 5550200. -5932855.

as.data.frame(as(as(atf_sp, "SpatialLinesDataFrame"),"SpatialPointsDataFrame"))
#>     iso_a3                   name sovereignt               continent
#> 8      ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#> 8.1    ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#> 8.2    ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#> 8.3    ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#> 8.4    ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#> 8.5    ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#> 8.6    ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#> 8.7    ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#> 8.8    ATF Fr. S. Antarctic Lands     France Seven seas (open ocean)
#>                area pop_est pop_est_dens              economy
#> 8   7257.455 [km^2]     140   0.01929051 6. Developing region
#> 8.1 7257.455 [km^2]     140   0.01929051 6. Developing region
#> 8.2 7257.455 [km^2]     140   0.01929051 6. Developing region
#> 8.3 7257.455 [km^2]     140   0.01929051 6. Developing region
#> 8.4 7257.455 [km^2]     140   0.01929051 6. Developing region
#> 8.5 7257.455 [km^2]     140   0.01929051 6. Developing region
#> 8.6 7257.455 [km^2]     140   0.01929051 6. Developing region
#> 8.7 7257.455 [km^2]     140   0.01929051 6. Developing region
#> 8.8 7257.455 [km^2]     140   0.01929051 6. Developing region
#>                  income_grp gdp_cap_est life_exp well_being footprint
#> 8   2. High income: nonOECD    114285.7       NA         NA        NA
#> 8.1 2. High income: nonOECD    114285.7       NA         NA        NA
#> 8.2 2. High income: nonOECD    114285.7       NA         NA        NA
#> 8.3 2. High income: nonOECD    114285.7       NA         NA        NA
#> 8.4 2. High income: nonOECD    114285.7       NA         NA        NA
#> 8.5 2. High income: nonOECD    114285.7       NA         NA        NA
#> 8.6 2. High income: nonOECD    114285.7       NA         NA        NA
#> 8.7 2. High income: nonOECD    114285.7       NA         NA        NA
#> 8.8 2. High income: nonOECD    114285.7       NA         NA        NA
#>     inequality HPI Lines.NR Lines.ID Line.NR coords.x1 coords.x2
#> 8           NA  NA        1        8       1   5550200  -5932855
#> 8.1         NA  NA        1        8       1   5589907  -5964836
#> 8.2         NA  NA        1        8       1   5660887  -5977490
#> 8.3         NA  NA        1        8       1   5656160  -5996685
#> 8.4         NA  NA        1        8       1   5615621  -6042456
#> 8.5         NA  NA        1        8       1   5490427  -6048972
#> 8.6         NA  NA        1        8       1   5509148  -5995424
#> 8.7         NA  NA        1        8       1   5536900  -5953683
#> 8.8         NA  NA        1        8       1   5550200  -5932855

# What about speed? raster::geom is the fastest
res <- microbenchmark(raster::geom(atf_sp),
                      ggplot2::fortify(atf_sp),
                      spbabel::sptable(atf_sp),
                      as.data.frame(as(as(atf_sp, "SpatialLinesDataFrame"),
                                       "SpatialPointsDataFrame")))
ggplot2::autoplot(res)
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.

由reprex包（v0.2.1）于2019年3月23日创建

0
此帖子中唯一有效的答案是由作者“repres_package”提供的。如果您想得到正确的答案，请查看该作者推荐的解决方案。如果您想获取多边形数据集的几何信息，则需要寻找多边形要素类中每个顶点的经度和纬度。例如，使用raster::geom()或ggplot2::fortify()等方法，将为您提供包含在spatialpolygonsdataframe中的所有顶点的总数。这就是您想要的。其他作者未能做到这一点。
例如，在我从美国人口普查获得的北卡罗来纳县的spatialpolygonsdataframe中，我有1259547个顶点。通过使用raster::geom(NC_counties)，我会得到一个数据框，其中包含这1259547个顶点的经度和纬度。我也可以使用gglot2::fortify(NC_counties)来获取这1259547个顶点的坐标。所有有效的选项都在“repres_package”的答案中给出。
当我运行了这篇帖子中其他答案中推荐的代码时，我只得到了672个顶点、1041个顶点或1721个顶点的经纬度坐标，与实际需要的1259547个顶点相差超过一百万个。我怀疑这些代码是在插值多边形的质心，而不是多边形的几何形状。
这并没有回答问题。一旦您拥有足够的声望，您将能够评论任何帖子；相反，提供不需要询问者澄清的答案。- 来自审核 - IRTFM
0

ggplot2 的 fortify() 函数可能在某个时候被弃用，因此现在建议采用 broom 包。

library(broom)
broom::tidy(atf_sp)
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接