使用多边形面积计算R中的人口密度
Use polygon area to calculate the population density in R
我是 R 的新手，我正在尝试制作密度 choroplet 地图。我想使用 rgal:
计算 R 中行政边界的人口密度
path <- "recintos_municipales_inspire_peninbal_etrs89.shp"
lineas_limite <- readOGR(path, use_iconv=TRUE, encoding = "UTF-8")
我看到了 rgal 如何创建一个数据结构，我可以在其中通过 lineas_limite@data 轻松访问 shapefile 属性 table 并通过 [=16= 访问多边形].
我在 lineas_limite@data 中有一个名为 population 的字段，我想再次使用它作为 QGIS 中可用的本机变量，例如 $area。
然后我想创建一个新的 density 变量来在等值线图中表示它：
coast@data <- coast@data %>%
mutate(
dens = population / area
)
这是最好的制作方法吗？
数据来自西班牙国家地理研究所（link来源）。
数据集名称为市界线。该平台将下载几个文件，我正在使用的文件是 recintos_provinciales_inspire_peninbal_etrs89.
另一方面，我已经用 shapefile 创建了一个 GitHub repo，但是那些已经有一个区域字段，这是我感兴趣的并且是我用 QGIS 创建的（所以不要请注意，因为我正在尝试从 R 生成它。
提前致谢！
概览
要计算每个多边形的面积，您需要做两件事：
- 存储每个多边形的边界（此处显示为 77 个矩阵的列表）
- 使用
geosphere::areaPolygon() 函数计算并存储每个多边形的面积。
从那里，您可以继续计算每个多边形的人口密度。
可重现的例子
这是一个计算每个 City of Chicago's 77 community areas（即街区）面积的示例。
# install necessary packages
install.packages( pkgs = c( "geosphere", "sp", "rgdal" ) )
# load necessary packages
library( geosphere )
library( sp )
library( rgdal )
# store Chicago current community area
# GeoJSON URL as a character vector
geojson_comarea_url <- "https://data.cityofchicago.org/api/geospatial/cauq-8yn6?method=export&format=GeoJSON"
# transform URL character vector into spatial dataframe
comarea606 <- readOGR( dsn = geojson_comarea_url
, layer = "OGRGeoJSON"
, stringsAsFactors = FALSE
, verbose = FALSE # to hide progress message after object is created
)
## WRONG WAY is trusting the area slot within each polygon
wrong.list.of.polygon.areas <-
data.frame(
comarea606_Name = comarea606$community
, comarea606_Inherit_Area = sapply( X = slot( object = comarea606, name = "polygons" )
, FUN = function(i)
i@area )
, stringsAsFactors = FALSE
)
# view the wrong results
wrong.list.of.polygon.areas
# comarea606_Name comarea606_Inherit_Area
# 1 DOUGLAS 0.0004632396
# 2 OAKLAND 0.0001702833
# 3 FULLER PARK 0.0002004699
# 4 GRAND BOULEVARD 0.0004881245
# 5 KENWOOD 0.0002926158
# 6 LINCOLN SQUARE 0.0007200422
# 7 WASHINGTON PARK 0.0004263996
# 8 HYDE PARK 0.0004538953
# 9 WOODLAWN 0.0005816583
# 10 ROGERS PARK 0.0005175561
# 11 JEFFERSON PARK 0.0006546563
# 12 FOREST GLEN 0.0008997264
# 13 NORTH PARK 0.0007094071
# 14 ALBANY PARK 0.0005402586
# 15 PORTAGE PARK 0.0011116825
# 16 IRVING PARK 0.0009040059
# 17 DUNNING 0.0010449760
# 18 MONTCLARE 0.0002780923
# 19 BELMONT CRAGIN 0.0011001642
# 20 WEST RIDGE 0.0009936912
# 21 HERMOSA 0.0003287441
# 22 AVONDALE 0.0005576454
# 23 LOGAN SQUARE 0.0010089020
# 24 HUMBOLDT PARK 0.0010128203
# 25 WEST TOWN 0.0012858100
# 26 AUSTIN 0.0020082609
# 27 WEST GARFIELD PARK 0.0003636864
# 28 EAST GARFIELD PARK 0.0005429495
# 29 NEAR WEST SIDE 0.0015968996
# 30 NORTH LAWNDALE 0.0009014540
# 31 UPTOWN 0.0006568040
# 32 SOUTH LAWNDALE 0.0012889745
# 33 LOWER WEST SIDE 0.0008213691
# 34 NEAR SOUTH SIDE 0.0005013217
# 35 ARMOUR SQUARE 0.0002796204
# 36 NORWOOD PARK 0.0012862114
# 37 NEAR NORTH SIDE 0.0007728518
# 38 LOOP 0.0004668873
# 39 SOUTH SHORE 0.0008228658
# 40 CHATHAM 0.0008277123
# 41 AVALON PARK 0.0003504538
# 42 SOUTH CHICAGO 0.0009378263
# 43 BURNSIDE 0.0001708577
# 44 CALUMET HEIGHTS 0.0004908500
# 45 ROSELAND 0.0013497948
# 46 NORTH CENTER 0.0005755101
# 47 PULLMAN 0.0005929323
# 48 SOUTH DEERING 0.0030522425
# 49 EAST SIDE 0.0008365337
# 50 WEST PULLMAN 0.0009980793
# 51 RIVERDALE 0.0009880640
# 52 HEGEWISCH 0.0014658303
# 53 GARFIELD RIDGE 0.0011864234
# 54 ARCHER HEIGHTS 0.0005629105
# 55 BRIGHTON PARK 0.0007640013
# 56 MCKINLEY PARK 0.0003970284
# 57 LAKE VIEW 0.0008796889
# 58 BRIDGEPORT 0.0005869749
# 59 NEW CITY 0.0013551840
# 60 WEST ELSDON 0.0003291759
# 61 GAGE PARK 0.0006167374
# 62 CLEARING 0.0007157968
# 63 WEST LAWN 0.0008280548
# 64 CHICAGO LAWN 0.0009886720
# 65 WEST ENGLEWOOD 0.0008847861
# 66 ENGLEWOOD 0.0008617058
# 67 GREATER GRAND CROSSING 0.0009942937
# 68 LINCOLN PARK 0.0008905024
# 69 ASHBURN 0.0013621619
# 70 AUBURN GRESHAM 0.0010564790
# 71 BEVERLY 0.0008922944
# 72 WASHINGTON HEIGHTS 0.0008004438
# 73 MOUNT GREENWOOD 0.0007594687
# 74 MORGAN PARK 0.0009230990
# 75 OHARE 0.0037524990
# 76 EDGEWATER 0.0004890111
# 77 EDISON PARK 0.0003194218
## RIGHT WAY is calculating the area yourself
# store each polygon's boundaries
list.of.polygon.boundaries <- sapply( X = slot( object = comarea606, name = "polygons" )
, FUN = function(i)
sp::coordinates( obj = slot( object = i, name = "Polygons")[[1]] )
)
# calculate each polygon's area in square meters
comarea606@data$Square_Meter_Area <- sapply( X = list.of.polygon.boundaries, FUN = geosphere::areaPolygon )
# view the right results
comarea606@data[ c( "community", "Square_Meter_Area" ) ]
# community Square_Meter_Area
# 1 DOUGLAS 4273829
# 2 OAKLAND 1571301
# 3 FULLER PARK 1850268
# 4 GRAND BOULEVARD 4504952
# 5 KENWOOD 2700750
# 6 LINCOLN SQUARE 6628739
# 7 WASHINGTON PARK 3936533
# 8 HYDE PARK 4190262
# 9 WOODLAWN 5370998
# 10 ROGERS PARK 4762104
# 11 JEFFERSON PARK 6026453
# 12 FOREST GLEN 8280516
# 13 NORTH PARK 6529977
# 14 ALBANY PARK 4974190
# 15 PORTAGE PARK 10237543
# 16 IRVING PARK 8325095
# 17 DUNNING 9624357
# 18 MONTCLARE 2561945
# 19 BELMONT CRAGIN 10135661
# 20 WEST RIDGE 9144227
# 21 HERMOSA 3028810
# 22 AVONDALE 5136605
# 23 LOGAN SQUARE 9295521
# 24 HUMBOLDT PARK 9334890
# 25 WEST TOWN 11850755
# 26 AUSTIN 18511300
# 27 WEST GARFIELD PARK 3353110
# 28 EAST GARFIELD PARK 5005851
# 29 NEAR WEST SIDE 14724107
# 30 NORTH LAWNDALE 8313567
# 31 UPTOWN 6047440
# 32 SOUTH LAWNDALE 11891299
# 33 LOWER WEST SIDE 7576149
# 34 NEAR SOUTH SIDE 4623602
# 35 ARMOUR SQUARE 2579489
# 36 NORWOOD PARK 11839116
# 37 NEAR NORTH SIDE 7123217
# 38 LOOP 4304581
# 39 SOUTH SHORE 7600312
# 40 CHATHAM 7647583
# 41 AVALON PARK 3237793
# 42 SOUTH CHICAGO 8664835
# 43 BURNSIDE 1578918
# 44 CALUMET HEIGHTS 4535904
# 45 ROSELAND 12477753
# 46 NORTH CENTER 5300414
# 47 PULLMAN 5481215
# 48 SOUTH DEERING 28222390
# 49 EAST SIDE 7732987
# 50 WEST PULLMAN 9231061
# 51 RIVERDALE 9140344
# 52 HEGEWISCH 13559960
# 53 GARFIELD RIDGE 10952400
# 54 ARCHER HEIGHTS 5195326
# 55 BRIGHTON PARK 7050570
# 56 MCKINLEY PARK 3663260
# 57 LAKE VIEW 8102329
# 58 BRIDGEPORT 5415321
# 59 NEW CITY 12507893
# 60 WEST ELSDON 3038932
# 61 GAGE PARK 5693468
# 62 CLEARING 6609603
# 63 WEST LAWN 7647276
# 64 CHICAGO LAWN 9130323
# 65 WEST ENGLEWOOD 8170431
# 66 ENGLEWOOD 7957145
# 67 GREATER GRAND CROSSING 9183452
# 68 LINCOLN PARK 8204655
# 69 ASHBURN 12584519
# 70 AUBURN GRESHAM 9760656
# 71 BEVERLY 8247709
# 72 WASHINGTON HEIGHTS 7398212
# 73 MOUNT GREENWOOD 7021929
# 74 MORGAN PARK 8535503
# 75 OHARE 34379638
# 76 EDGEWATER 4501053
# 77 EDISON PARK 2939134
# end of script #
不信任每个多边形中的继承 'area' 插槽
根据 ?sp::Polygons-class 的文档，每个多边形中的 area 槽不是实际区域。
Object of class "numeric"; the gross total planar area of the Polygon list but not double-counting holes (changed from 0.9-58 - islands are summed, holes are ignored rather than subtracted); these values are used to make sure that polygons of a smaller area are plotted after polygons of a larger area, does not respect projection as objects of this class have no projection defined.
它看起来更像是一个占位符，而不是您要查找的实际值，这就是为什么我使用 geosphere::areaPolygon() 函数来获取该值的原因。
更新
OP 的数据是通过 GitHub repository 提供的。下面是与上面相同的逻辑，只是使用了不同的数据。
# install necessary packages
install.packages( pkgs = c( "geosphere", "sp", "rgdal" ) )
# load necessary packages
library( geosphere )
library( sp )
library( rgdal )
# load necessary data
# download a .zip file of the GitHub repository
download.file( url = "https://github.com/LuisSevillano/lineas-limite/archive/master.zip"
, destfile = "lineas-limite-master.zip"
) #39.2MB ~ 1 minute DL time on my crappy wifi
# unzip the entire GitHub repository
unzip( zipfile = "lineas-limite-master.zip" )
# set working directory to be
# inside the sub-folder where
# the shapefiles of interest are located
setwd( dir = "/Users/cristiannuno/Downloads/lineas-limite-master/recintos_municipales_inspire_peninbal_etrs89/" )
# create spatial polygon data frame
lineas_limite <- rgdal::readOGR( dsn = getwd()
, layer = "recintos_municipales_inspire_peninbal_etrs89"
, stringsAsFactors = FALSE
)
# store each polygon's boundaries
list.of.polygon.boundaries <- sapply( X = slot( object = lineas_limite, name = "polygons" )
, FUN = function(i)
sp::coordinates( obj = slot( object = i, name = "Polygons")[[1]] )
)
# calculate each polygon's area in square meters
lineas_limite@data$Square_Meter_Area <- sapply( X = list.of.polygon.boundaries
, FUN = geosphere::areaPolygon
)
# View the first six results
head(
lineas_limite@data[ c( "NAMEUNIT", "Square_Meter_Area" ) ]
)
# NAMEUNIT Square_Meter_Area
# 0 Abla 45244459
# 1 Abrucena 83677254
# 2 Adra 89694443
# 3 Albánchez 35133840
# 4 Alboloduy 69734542
# 5 Pradejón 31755520
# end of script #
资源
我鼓励您通读 ?sp::Polygons-class 文档以及 SpatialPointsDataFrame 属性和 R 中的运算符，以了解有关空间多边形及其使用方法的更多信息在 R.
Session 信息
R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] rgdal_1.2-16 geosphere_1.5-7 sp_1.2-7
loaded via a namespace (and not attached):
[1] Rcpp_0.12.15 lattice_0.20-35 digest_0.6.14 mime_0.5
[5] grid_3.4.3 R6_2.2.2 xtable_1.8-2 magrittr_1.5
[9] leaflet_1.1.0.9000 tools_3.4.3 htmlwidgets_1.0 crosstalk_1.0.0
[13] shiny_1.0.5 httpuv_1.3.5 yaml_2.1.16 compiler_3.4.3
[17] htmltools_0.3.6
RStudio 版本
$citation
To cite RStudio in publications use:
RStudio Team (2016). RStudio: Integrated Development for R. RStudio,
Inc., Boston, MA URL http://www.rstudio.com/.
A BibTeX entry for LaTeX users is
@Manual{,
title = {RStudio: Integrated Development Environment for R},
author = {{RStudio Team}},
organization = {RStudio, Inc.},
address = {Boston, MA},
year = {2016},
url = {http://www.rstudio.com/},
}
$mode
[1] "desktop"
$version
[1] ‘1.1.419’
我是 R 的新手，我正在尝试制作密度 choroplet 地图。我想使用 rgal:
path <- "recintos_municipales_inspire_peninbal_etrs89.shp"
lineas_limite <- readOGR(path, use_iconv=TRUE, encoding = "UTF-8")
我看到了 rgal 如何创建一个数据结构，我可以在其中通过 lineas_limite@data 轻松访问 shapefile 属性 table 并通过 [=16= 访问多边形].
我在 lineas_limite@data 中有一个名为 population 的字段，我想再次使用它作为 QGIS 中可用的本机变量，例如 $area。
然后我想创建一个新的 density 变量来在等值线图中表示它：
coast@data <- coast@data %>%
mutate(
dens = population / area
)
这是最好的制作方法吗？
数据来自西班牙国家地理研究所（link来源）。
数据集名称为市界线。该平台将下载几个文件，我正在使用的文件是 recintos_provinciales_inspire_peninbal_etrs89.
另一方面，我已经用 shapefile 创建了一个 GitHub repo，但是那些已经有一个区域字段，这是我感兴趣的并且是我用 QGIS 创建的（所以不要请注意，因为我正在尝试从 R 生成它。
提前致谢！
概览
要计算每个多边形的面积，您需要做两件事：
- 存储每个多边形的边界（此处显示为 77 个矩阵的列表）
- 使用
geosphere::areaPolygon()函数计算并存储每个多边形的面积。
从那里，您可以继续计算每个多边形的人口密度。
可重现的例子
这是一个计算每个 City of Chicago's 77 community areas（即街区）面积的示例。
# install necessary packages
install.packages( pkgs = c( "geosphere", "sp", "rgdal" ) )
# load necessary packages
library( geosphere )
library( sp )
library( rgdal )
# store Chicago current community area
# GeoJSON URL as a character vector
geojson_comarea_url <- "https://data.cityofchicago.org/api/geospatial/cauq-8yn6?method=export&format=GeoJSON"
# transform URL character vector into spatial dataframe
comarea606 <- readOGR( dsn = geojson_comarea_url
, layer = "OGRGeoJSON"
, stringsAsFactors = FALSE
, verbose = FALSE # to hide progress message after object is created
)
## WRONG WAY is trusting the area slot within each polygon
wrong.list.of.polygon.areas <-
data.frame(
comarea606_Name = comarea606$community
, comarea606_Inherit_Area = sapply( X = slot( object = comarea606, name = "polygons" )
, FUN = function(i)
i@area )
, stringsAsFactors = FALSE
)
# view the wrong results
wrong.list.of.polygon.areas
# comarea606_Name comarea606_Inherit_Area
# 1 DOUGLAS 0.0004632396
# 2 OAKLAND 0.0001702833
# 3 FULLER PARK 0.0002004699
# 4 GRAND BOULEVARD 0.0004881245
# 5 KENWOOD 0.0002926158
# 6 LINCOLN SQUARE 0.0007200422
# 7 WASHINGTON PARK 0.0004263996
# 8 HYDE PARK 0.0004538953
# 9 WOODLAWN 0.0005816583
# 10 ROGERS PARK 0.0005175561
# 11 JEFFERSON PARK 0.0006546563
# 12 FOREST GLEN 0.0008997264
# 13 NORTH PARK 0.0007094071
# 14 ALBANY PARK 0.0005402586
# 15 PORTAGE PARK 0.0011116825
# 16 IRVING PARK 0.0009040059
# 17 DUNNING 0.0010449760
# 18 MONTCLARE 0.0002780923
# 19 BELMONT CRAGIN 0.0011001642
# 20 WEST RIDGE 0.0009936912
# 21 HERMOSA 0.0003287441
# 22 AVONDALE 0.0005576454
# 23 LOGAN SQUARE 0.0010089020
# 24 HUMBOLDT PARK 0.0010128203
# 25 WEST TOWN 0.0012858100
# 26 AUSTIN 0.0020082609
# 27 WEST GARFIELD PARK 0.0003636864
# 28 EAST GARFIELD PARK 0.0005429495
# 29 NEAR WEST SIDE 0.0015968996
# 30 NORTH LAWNDALE 0.0009014540
# 31 UPTOWN 0.0006568040
# 32 SOUTH LAWNDALE 0.0012889745
# 33 LOWER WEST SIDE 0.0008213691
# 34 NEAR SOUTH SIDE 0.0005013217
# 35 ARMOUR SQUARE 0.0002796204
# 36 NORWOOD PARK 0.0012862114
# 37 NEAR NORTH SIDE 0.0007728518
# 38 LOOP 0.0004668873
# 39 SOUTH SHORE 0.0008228658
# 40 CHATHAM 0.0008277123
# 41 AVALON PARK 0.0003504538
# 42 SOUTH CHICAGO 0.0009378263
# 43 BURNSIDE 0.0001708577
# 44 CALUMET HEIGHTS 0.0004908500
# 45 ROSELAND 0.0013497948
# 46 NORTH CENTER 0.0005755101
# 47 PULLMAN 0.0005929323
# 48 SOUTH DEERING 0.0030522425
# 49 EAST SIDE 0.0008365337
# 50 WEST PULLMAN 0.0009980793
# 51 RIVERDALE 0.0009880640
# 52 HEGEWISCH 0.0014658303
# 53 GARFIELD RIDGE 0.0011864234
# 54 ARCHER HEIGHTS 0.0005629105
# 55 BRIGHTON PARK 0.0007640013
# 56 MCKINLEY PARK 0.0003970284
# 57 LAKE VIEW 0.0008796889
# 58 BRIDGEPORT 0.0005869749
# 59 NEW CITY 0.0013551840
# 60 WEST ELSDON 0.0003291759
# 61 GAGE PARK 0.0006167374
# 62 CLEARING 0.0007157968
# 63 WEST LAWN 0.0008280548
# 64 CHICAGO LAWN 0.0009886720
# 65 WEST ENGLEWOOD 0.0008847861
# 66 ENGLEWOOD 0.0008617058
# 67 GREATER GRAND CROSSING 0.0009942937
# 68 LINCOLN PARK 0.0008905024
# 69 ASHBURN 0.0013621619
# 70 AUBURN GRESHAM 0.0010564790
# 71 BEVERLY 0.0008922944
# 72 WASHINGTON HEIGHTS 0.0008004438
# 73 MOUNT GREENWOOD 0.0007594687
# 74 MORGAN PARK 0.0009230990
# 75 OHARE 0.0037524990
# 76 EDGEWATER 0.0004890111
# 77 EDISON PARK 0.0003194218
## RIGHT WAY is calculating the area yourself
# store each polygon's boundaries
list.of.polygon.boundaries <- sapply( X = slot( object = comarea606, name = "polygons" )
, FUN = function(i)
sp::coordinates( obj = slot( object = i, name = "Polygons")[[1]] )
)
# calculate each polygon's area in square meters
comarea606@data$Square_Meter_Area <- sapply( X = list.of.polygon.boundaries, FUN = geosphere::areaPolygon )
# view the right results
comarea606@data[ c( "community", "Square_Meter_Area" ) ]
# community Square_Meter_Area
# 1 DOUGLAS 4273829
# 2 OAKLAND 1571301
# 3 FULLER PARK 1850268
# 4 GRAND BOULEVARD 4504952
# 5 KENWOOD 2700750
# 6 LINCOLN SQUARE 6628739
# 7 WASHINGTON PARK 3936533
# 8 HYDE PARK 4190262
# 9 WOODLAWN 5370998
# 10 ROGERS PARK 4762104
# 11 JEFFERSON PARK 6026453
# 12 FOREST GLEN 8280516
# 13 NORTH PARK 6529977
# 14 ALBANY PARK 4974190
# 15 PORTAGE PARK 10237543
# 16 IRVING PARK 8325095
# 17 DUNNING 9624357
# 18 MONTCLARE 2561945
# 19 BELMONT CRAGIN 10135661
# 20 WEST RIDGE 9144227
# 21 HERMOSA 3028810
# 22 AVONDALE 5136605
# 23 LOGAN SQUARE 9295521
# 24 HUMBOLDT PARK 9334890
# 25 WEST TOWN 11850755
# 26 AUSTIN 18511300
# 27 WEST GARFIELD PARK 3353110
# 28 EAST GARFIELD PARK 5005851
# 29 NEAR WEST SIDE 14724107
# 30 NORTH LAWNDALE 8313567
# 31 UPTOWN 6047440
# 32 SOUTH LAWNDALE 11891299
# 33 LOWER WEST SIDE 7576149
# 34 NEAR SOUTH SIDE 4623602
# 35 ARMOUR SQUARE 2579489
# 36 NORWOOD PARK 11839116
# 37 NEAR NORTH SIDE 7123217
# 38 LOOP 4304581
# 39 SOUTH SHORE 7600312
# 40 CHATHAM 7647583
# 41 AVALON PARK 3237793
# 42 SOUTH CHICAGO 8664835
# 43 BURNSIDE 1578918
# 44 CALUMET HEIGHTS 4535904
# 45 ROSELAND 12477753
# 46 NORTH CENTER 5300414
# 47 PULLMAN 5481215
# 48 SOUTH DEERING 28222390
# 49 EAST SIDE 7732987
# 50 WEST PULLMAN 9231061
# 51 RIVERDALE 9140344
# 52 HEGEWISCH 13559960
# 53 GARFIELD RIDGE 10952400
# 54 ARCHER HEIGHTS 5195326
# 55 BRIGHTON PARK 7050570
# 56 MCKINLEY PARK 3663260
# 57 LAKE VIEW 8102329
# 58 BRIDGEPORT 5415321
# 59 NEW CITY 12507893
# 60 WEST ELSDON 3038932
# 61 GAGE PARK 5693468
# 62 CLEARING 6609603
# 63 WEST LAWN 7647276
# 64 CHICAGO LAWN 9130323
# 65 WEST ENGLEWOOD 8170431
# 66 ENGLEWOOD 7957145
# 67 GREATER GRAND CROSSING 9183452
# 68 LINCOLN PARK 8204655
# 69 ASHBURN 12584519
# 70 AUBURN GRESHAM 9760656
# 71 BEVERLY 8247709
# 72 WASHINGTON HEIGHTS 7398212
# 73 MOUNT GREENWOOD 7021929
# 74 MORGAN PARK 8535503
# 75 OHARE 34379638
# 76 EDGEWATER 4501053
# 77 EDISON PARK 2939134
# end of script #
不信任每个多边形中的继承 'area' 插槽
根据 ?sp::Polygons-class 的文档，每个多边形中的 area 槽不是实际区域。
Object of class "numeric"; the gross total planar area of the Polygon list but not double-counting holes (changed from 0.9-58 - islands are summed, holes are ignored rather than subtracted); these values are used to make sure that polygons of a smaller area are plotted after polygons of a larger area, does not respect projection as objects of this class have no projection defined.
它看起来更像是一个占位符，而不是您要查找的实际值，这就是为什么我使用 geosphere::areaPolygon() 函数来获取该值的原因。
更新
OP 的数据是通过 GitHub repository 提供的。下面是与上面相同的逻辑，只是使用了不同的数据。
# install necessary packages
install.packages( pkgs = c( "geosphere", "sp", "rgdal" ) )
# load necessary packages
library( geosphere )
library( sp )
library( rgdal )
# load necessary data
# download a .zip file of the GitHub repository
download.file( url = "https://github.com/LuisSevillano/lineas-limite/archive/master.zip"
, destfile = "lineas-limite-master.zip"
) #39.2MB ~ 1 minute DL time on my crappy wifi
# unzip the entire GitHub repository
unzip( zipfile = "lineas-limite-master.zip" )
# set working directory to be
# inside the sub-folder where
# the shapefiles of interest are located
setwd( dir = "/Users/cristiannuno/Downloads/lineas-limite-master/recintos_municipales_inspire_peninbal_etrs89/" )
# create spatial polygon data frame
lineas_limite <- rgdal::readOGR( dsn = getwd()
, layer = "recintos_municipales_inspire_peninbal_etrs89"
, stringsAsFactors = FALSE
)
# store each polygon's boundaries
list.of.polygon.boundaries <- sapply( X = slot( object = lineas_limite, name = "polygons" )
, FUN = function(i)
sp::coordinates( obj = slot( object = i, name = "Polygons")[[1]] )
)
# calculate each polygon's area in square meters
lineas_limite@data$Square_Meter_Area <- sapply( X = list.of.polygon.boundaries
, FUN = geosphere::areaPolygon
)
# View the first six results
head(
lineas_limite@data[ c( "NAMEUNIT", "Square_Meter_Area" ) ]
)
# NAMEUNIT Square_Meter_Area
# 0 Abla 45244459
# 1 Abrucena 83677254
# 2 Adra 89694443
# 3 Albánchez 35133840
# 4 Alboloduy 69734542
# 5 Pradejón 31755520
# end of script #
资源
我鼓励您通读 ?sp::Polygons-class 文档以及 SpatialPointsDataFrame 属性和 R 中的运算符，以了解有关空间多边形及其使用方法的更多信息在 R.
Session 信息
R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] rgdal_1.2-16 geosphere_1.5-7 sp_1.2-7
loaded via a namespace (and not attached):
[1] Rcpp_0.12.15 lattice_0.20-35 digest_0.6.14 mime_0.5
[5] grid_3.4.3 R6_2.2.2 xtable_1.8-2 magrittr_1.5
[9] leaflet_1.1.0.9000 tools_3.4.3 htmlwidgets_1.0 crosstalk_1.0.0
[13] shiny_1.0.5 httpuv_1.3.5 yaml_2.1.16 compiler_3.4.3
[17] htmltools_0.3.6
RStudio 版本
$citation
To cite RStudio in publications use:
RStudio Team (2016). RStudio: Integrated Development for R. RStudio,
Inc., Boston, MA URL http://www.rstudio.com/.
A BibTeX entry for LaTeX users is
@Manual{,
title = {RStudio: Integrated Development Environment for R},
author = {{RStudio Team}},
organization = {RStudio, Inc.},
address = {Boston, MA},
year = {2016},
url = {http://www.rstudio.com/},
}
$mode
[1] "desktop"
$version
[1] ‘1.1.419’