在多列上聚合 sf 对象
Aggregating sf object over multiple columns
我正在处理一些地理参考工业数据。坐标系指某产业项目所在城市的质心。
structure(list(date = c("01dec2013", "01jul2003", "01nov2008",
"01dec2017", "01dec2017", "01dec2003"), company = c("Shwe Taung",
"PetroChina Exploration and Development", "Repsol SA", "Repsol SA",
"Ipsen Pharmaceutical", "Ceva Laval"), parent_company = c("Shwe Taung",
"China National Petroleum (CNPC)", "Repsol SA", "Repsol SA",
"Ipsen Pharmaceutical", "Ceva Sante Animale"), Website = c("www.shwetaunggroup.com",
"www.cnpc.com.cn", "www.repsol.com", "www.repsol.com", "www.ipsen.com",
"www.ceva.com"), revenues_usd_ml = c(NA, 394554.53, 53215.45,
53215.45, 1760.671, 967.152), Headcount = c(NA, 1396144L, 24634L,
24634L, NA, 3500L), r_d_exp = c(NA, NA, 77.67, 77.67, NA, NA),
est_year = c(NA, 1988L, 1927L, 1927L, 1929L, 1989L), o_country = c("Myanmar",
"China", "Spain", "Spain", "France", "France"), o_state = c("Rangoon (Yangon)",
"Beijing Municipality", "Comunidad de Madrid", "Comunidad de Madrid",
"Ile-de-France", "Sud-Ouest (FR)"), o_admin = c("Not Specified",
"Not Specified", "Madrid", "Madrid", "Ile-de-France", "Not Specified"
), o_city = c("Rangoon (Yangon)", "Beijing", "Madrid", "Madrid",
"Paris", "Not Specified"), country = c("Algeria", "Algeria",
"Algeria", "Algeria", "Algeria", "Algeria"), state = c("Adrar",
"Adrar", "Adrar", "Adrar", "Adrar", "Adrar"), region = c("Not Specified",
"Not Specified", "Not Specified", "Not Specified", "Not Specified",
"Not Specified"), city = c("Adrar", "Adrar", "Reggane", "Reggane",
"Sidi Abdallah", "Sidi Abdallah"), free_zone = c("", "",
"", "", "", ""), relocation = c("", "", "", "", "", ""),
sector = c("Building materials", "Coal, oil & gas", "Coal, oil & gas",
"Coal, oil & gas", "Pharmaceuticals", "Healthcare"), sub_sector = c("Cement & concrete products",
"Oil & gas extraction", "Oil & gas extraction", "Oil & gas extraction",
"Pharmaceutical preparations", "Other (Healthcare)"), cluster = c("Construction",
"Energy", "Energy", "Energy", "Life sciences", "Life sciences"
), activity = c("Manufacturing", "Extraction", "Extraction",
"Extraction", "Manufacturing", "Manufacturing"), fdi_jobs = c(351L,
145L, 235L, 227L, 150L, 45L), est_fdi_jobs = c("Yes", "Yes",
"Yes", "Yes", "No", "No"), capital = c(139.9, 350, 565, 299.7,
29.55, 2.5), est_capital = c("Yes", "No", "No", "Yes", "No",
"No"), fdi_type = c("New", "New", "New", "Expansion", "New",
"New"), fdi_status = c("Announced", "Announced", "Announced",
"Opened", "Announced", "Opened"), source = c("geocode", "geocode",
"opencagegeo", "opencagegeo", "opencagegeo", "opencagegeo"
), year = c(2013L, 2003L, 2008L, 2017L, 2017L, 2003L), code_d = c(12L,
12L, 12L, 12L, 12L, 12L), income_d = c("MIDLW", "MIDLW",
"MIDLW", "MIDLW", "MIDLW", "MIDLW"), continent_d = c("Africa",
"Africa", "Africa", "Africa", "Africa", "Africa"), lang_d = c("Arabic",
"Arabic", "Arabic", "Arabic", "Arabic", "Arabic"), landlocked = c(0L,
0L, 0L, 0L, 0L, 0L), iso_d = c("DZA", "DZA", "DZA", "DZA",
"DZA", "DZA"), isic = c(26L, 11L, 11L, 11L, 24L, 85L), isic4 = c(2695L,
1110L, 1110L, 1110L, 2411L, 8519L), sector_eora = c("Petroleum, Chemical and Non-Metallic Mineral Products",
"Mining and Quarrying", "Mining and Quarrying", "Mining and Quarrying",
"Petroleum, Chemical and Non-Metallic Mineral Products",
"Mining and Quarrying"), n_proj = c(1, 1, 1, 1, 1, 1), geometry = structure(list(
structure(c(-0.287488, 27.87338), class = c("XY", "POINT",
"sfg")), structure(c(-0.287488, 27.87338), class = c("XY",
"POINT", "sfg")), structure(c(0.1751507, 26.7207267), class = c("XY",
"POINT", "sfg")), structure(c(0.1751507, 26.7207267), class = c("XY",
"POINT", "sfg")), structure(c(5.0152099, 35.726173), class = c("XY",
"POINT", "sfg")), structure(c(5.0152099, 35.726173), class = c("XY",
"POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(xmin = -0.287488,
ymin = 26.7207267, xmax = 5.0152099, ymax = 35.726173), class = "bbox"), crs = structure(list(
input = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",
wkt = "BOUNDCRS[\n SOURCECRS[\n GEOGCRS[\"unknown\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",6326]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8901]],\n CS[ellipsoidal,2],\n AXIS[\"longitude\",east,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"latitude\",north,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]],\n TARGETCRS[\n GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]],\n ABRIDGEDTRANSFORMATION[\"Transformation from unknown to WGS84\",\n METHOD[\"Geocentric translations (geog2D domain)\",\n ID[\"EPSG\",9603]],\n PARAMETER[\"X-axis translation\",0,\n ID[\"EPSG\",8605]],\n PARAMETER[\"Y-axis translation\",0,\n ID[\"EPSG\",8606]],\n PARAMETER[\"Z-axis translation\",0,\n ID[\"EPSG\",8607]]]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(date = NA_integer_,
company = NA_integer_, parent_company = NA_integer_, Website = NA_integer_,
revenues_usd_ml = NA_integer_, Headcount = NA_integer_, r_d_exp = NA_integer_,
est_year = NA_integer_, o_country = NA_integer_, o_state = NA_integer_,
o_admin = NA_integer_, o_city = NA_integer_, country = NA_integer_,
state = NA_integer_, region = NA_integer_, city = NA_integer_,
free_zone = NA_integer_, relocation = NA_integer_, sector = NA_integer_,
sub_sector = NA_integer_, cluster = NA_integer_, activity = NA_integer_,
fdi_jobs = NA_integer_, est_fdi_jobs = NA_integer_, capital = NA_integer_,
est_capital = NA_integer_, fdi_type = NA_integer_, fdi_status = NA_integer_,
source = NA_integer_, year = NA_integer_, code_d = NA_integer_,
income_d = NA_integer_, continent_d = NA_integer_, lang_d = NA_integer_,
landlocked = NA_integer_, iso_d = NA_integer_, isic = NA_integer_,
isic4 = NA_integer_, sector_eora = NA_integer_, n_proj = NA_integer_
), .Label = c("constant", "aggregate", "identity"), class = "factor"), row.names = c(NA,
6L), class = c("sf", "data.frame"))
我想统计每个部门（变量：sector）每个点城市（变量：geometry）的项目数（变量： n_proj）。
我的想法是aggregate两列：
FDI_sf_ag1 <- aggregate(FDI_sf["n_proj"] ~ FDI_sf["geometry"] + FDI_sf["sector"], FDI_sf, sum)
但是它返回了这个错误：
invalid type (list) for variable 'FDI_sf["n_proj"]'
然后，我重新尝试取消列出列：
FDI_sf_ag1 <- aggregate(unlist(FDI_sf["n_proj"]) ~ cbind(unlist(FDI_sf["geometry"]),unlist(FDI_sf["sector"])), FDI_sf, sum)
但它返回另一个错误：
Error in aggregate.data.frame(lhs, mf[-1L], FUN = FUN, ...) : arguments must have same length In addition: Warning message: In cbind(unlist(FDI_sf["geometry"]), unlist(FDI_sf["sector"])) : number of rows of result is not a multiple of vector length (arg 1)
问题：如何聚合 n_project 几何（sfc 对象）和扇区？
aggregate(fdi_sf[, "n_proj"],
with(fdi_sf, list(sector=sector, city=city)),
sum)
Simple feature collection with 5 features and 3 fields
Attribute-geometry relationship: 0 constant, 1 aggregate, 2 identity
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -0.287488 ymin: 26.72073 xmax: 5.01521 ymax: 35.72617
CRS: +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
sector city n_proj
1 Building materials Adrar 1
2 Coal, oil & gas Adrar 1
3 Coal, oil & gas Reggane 2
4 Healthcare Sidi Abdallah 1
5 Pharmaceuticals Sidi Abdallah 1
geometry
1 POINT (-0.287488 27.87338)
2 POINT (-0.287488 27.87338)
3 POINT (0.1751507 26.72073)
4 POINT (5.01521 35.72617)
5 POINT (5.01521 35.72617)
我正在处理一些地理参考工业数据。坐标系指某产业项目所在城市的质心。
structure(list(date = c("01dec2013", "01jul2003", "01nov2008",
"01dec2017", "01dec2017", "01dec2003"), company = c("Shwe Taung",
"PetroChina Exploration and Development", "Repsol SA", "Repsol SA",
"Ipsen Pharmaceutical", "Ceva Laval"), parent_company = c("Shwe Taung",
"China National Petroleum (CNPC)", "Repsol SA", "Repsol SA",
"Ipsen Pharmaceutical", "Ceva Sante Animale"), Website = c("www.shwetaunggroup.com",
"www.cnpc.com.cn", "www.repsol.com", "www.repsol.com", "www.ipsen.com",
"www.ceva.com"), revenues_usd_ml = c(NA, 394554.53, 53215.45,
53215.45, 1760.671, 967.152), Headcount = c(NA, 1396144L, 24634L,
24634L, NA, 3500L), r_d_exp = c(NA, NA, 77.67, 77.67, NA, NA),
est_year = c(NA, 1988L, 1927L, 1927L, 1929L, 1989L), o_country = c("Myanmar",
"China", "Spain", "Spain", "France", "France"), o_state = c("Rangoon (Yangon)",
"Beijing Municipality", "Comunidad de Madrid", "Comunidad de Madrid",
"Ile-de-France", "Sud-Ouest (FR)"), o_admin = c("Not Specified",
"Not Specified", "Madrid", "Madrid", "Ile-de-France", "Not Specified"
), o_city = c("Rangoon (Yangon)", "Beijing", "Madrid", "Madrid",
"Paris", "Not Specified"), country = c("Algeria", "Algeria",
"Algeria", "Algeria", "Algeria", "Algeria"), state = c("Adrar",
"Adrar", "Adrar", "Adrar", "Adrar", "Adrar"), region = c("Not Specified",
"Not Specified", "Not Specified", "Not Specified", "Not Specified",
"Not Specified"), city = c("Adrar", "Adrar", "Reggane", "Reggane",
"Sidi Abdallah", "Sidi Abdallah"), free_zone = c("", "",
"", "", "", ""), relocation = c("", "", "", "", "", ""),
sector = c("Building materials", "Coal, oil & gas", "Coal, oil & gas",
"Coal, oil & gas", "Pharmaceuticals", "Healthcare"), sub_sector = c("Cement & concrete products",
"Oil & gas extraction", "Oil & gas extraction", "Oil & gas extraction",
"Pharmaceutical preparations", "Other (Healthcare)"), cluster = c("Construction",
"Energy", "Energy", "Energy", "Life sciences", "Life sciences"
), activity = c("Manufacturing", "Extraction", "Extraction",
"Extraction", "Manufacturing", "Manufacturing"), fdi_jobs = c(351L,
145L, 235L, 227L, 150L, 45L), est_fdi_jobs = c("Yes", "Yes",
"Yes", "Yes", "No", "No"), capital = c(139.9, 350, 565, 299.7,
29.55, 2.5), est_capital = c("Yes", "No", "No", "Yes", "No",
"No"), fdi_type = c("New", "New", "New", "Expansion", "New",
"New"), fdi_status = c("Announced", "Announced", "Announced",
"Opened", "Announced", "Opened"), source = c("geocode", "geocode",
"opencagegeo", "opencagegeo", "opencagegeo", "opencagegeo"
), year = c(2013L, 2003L, 2008L, 2017L, 2017L, 2003L), code_d = c(12L,
12L, 12L, 12L, 12L, 12L), income_d = c("MIDLW", "MIDLW",
"MIDLW", "MIDLW", "MIDLW", "MIDLW"), continent_d = c("Africa",
"Africa", "Africa", "Africa", "Africa", "Africa"), lang_d = c("Arabic",
"Arabic", "Arabic", "Arabic", "Arabic", "Arabic"), landlocked = c(0L,
0L, 0L, 0L, 0L, 0L), iso_d = c("DZA", "DZA", "DZA", "DZA",
"DZA", "DZA"), isic = c(26L, 11L, 11L, 11L, 24L, 85L), isic4 = c(2695L,
1110L, 1110L, 1110L, 2411L, 8519L), sector_eora = c("Petroleum, Chemical and Non-Metallic Mineral Products",
"Mining and Quarrying", "Mining and Quarrying", "Mining and Quarrying",
"Petroleum, Chemical and Non-Metallic Mineral Products",
"Mining and Quarrying"), n_proj = c(1, 1, 1, 1, 1, 1), geometry = structure(list(
structure(c(-0.287488, 27.87338), class = c("XY", "POINT",
"sfg")), structure(c(-0.287488, 27.87338), class = c("XY",
"POINT", "sfg")), structure(c(0.1751507, 26.7207267), class = c("XY",
"POINT", "sfg")), structure(c(0.1751507, 26.7207267), class = c("XY",
"POINT", "sfg")), structure(c(5.0152099, 35.726173), class = c("XY",
"POINT", "sfg")), structure(c(5.0152099, 35.726173), class = c("XY",
"POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(xmin = -0.287488,
ymin = 26.7207267, xmax = 5.0152099, ymax = 35.726173), class = "bbox"), crs = structure(list(
input = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0",
wkt = "BOUNDCRS[\n SOURCECRS[\n GEOGCRS[\"unknown\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",6326]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8901]],\n CS[ellipsoidal,2],\n AXIS[\"longitude\",east,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"latitude\",north,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]],\n TARGETCRS[\n GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]],\n ABRIDGEDTRANSFORMATION[\"Transformation from unknown to WGS84\",\n METHOD[\"Geocentric translations (geog2D domain)\",\n ID[\"EPSG\",9603]],\n PARAMETER[\"X-axis translation\",0,\n ID[\"EPSG\",8605]],\n PARAMETER[\"Y-axis translation\",0,\n ID[\"EPSG\",8606]],\n PARAMETER[\"Z-axis translation\",0,\n ID[\"EPSG\",8607]]]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(date = NA_integer_,
company = NA_integer_, parent_company = NA_integer_, Website = NA_integer_,
revenues_usd_ml = NA_integer_, Headcount = NA_integer_, r_d_exp = NA_integer_,
est_year = NA_integer_, o_country = NA_integer_, o_state = NA_integer_,
o_admin = NA_integer_, o_city = NA_integer_, country = NA_integer_,
state = NA_integer_, region = NA_integer_, city = NA_integer_,
free_zone = NA_integer_, relocation = NA_integer_, sector = NA_integer_,
sub_sector = NA_integer_, cluster = NA_integer_, activity = NA_integer_,
fdi_jobs = NA_integer_, est_fdi_jobs = NA_integer_, capital = NA_integer_,
est_capital = NA_integer_, fdi_type = NA_integer_, fdi_status = NA_integer_,
source = NA_integer_, year = NA_integer_, code_d = NA_integer_,
income_d = NA_integer_, continent_d = NA_integer_, lang_d = NA_integer_,
landlocked = NA_integer_, iso_d = NA_integer_, isic = NA_integer_,
isic4 = NA_integer_, sector_eora = NA_integer_, n_proj = NA_integer_
), .Label = c("constant", "aggregate", "identity"), class = "factor"), row.names = c(NA,
6L), class = c("sf", "data.frame"))
我想统计每个部门（变量：sector）每个点城市（变量：geometry）的项目数（变量： n_proj）。
我的想法是aggregate两列：
FDI_sf_ag1 <- aggregate(FDI_sf["n_proj"] ~ FDI_sf["geometry"] + FDI_sf["sector"], FDI_sf, sum)
但是它返回了这个错误：
invalid type (list) for variable 'FDI_sf["n_proj"]'
然后，我重新尝试取消列出列：
FDI_sf_ag1 <- aggregate(unlist(FDI_sf["n_proj"]) ~ cbind(unlist(FDI_sf["geometry"]),unlist(FDI_sf["sector"])), FDI_sf, sum)
但它返回另一个错误：
Error in aggregate.data.frame(lhs, mf[-1L], FUN = FUN, ...) : arguments must have same length In addition: Warning message: In cbind(unlist(FDI_sf["geometry"]), unlist(FDI_sf["sector"])) : number of rows of result is not a multiple of vector length (arg 1)
问题：如何聚合 n_project 几何（sfc 对象）和扇区？
aggregate(fdi_sf[, "n_proj"],
with(fdi_sf, list(sector=sector, city=city)),
sum)
Simple feature collection with 5 features and 3 fields
Attribute-geometry relationship: 0 constant, 1 aggregate, 2 identity
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -0.287488 ymin: 26.72073 xmax: 5.01521 ymax: 35.72617
CRS: +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
sector city n_proj
1 Building materials Adrar 1
2 Coal, oil & gas Adrar 1
3 Coal, oil & gas Reggane 2
4 Healthcare Sidi Abdallah 1
5 Pharmaceuticals Sidi Abdallah 1
geometry
1 POINT (-0.287488 27.87338)
2 POINT (-0.287488 27.87338)
3 POINT (0.1751507 26.72073)
4 POINT (5.01521 35.72617)
5 POINT (5.01521 35.72617)