计算 shapefile 中的密度
Calculating density within a shapefile
我正在尝试计算 shapefile 中的密度,但我很确定我做错了。这个想法是按密度找出哪些地理区域的销售额最高。
Here is a link to the file that I use (testdata.shp)
library(sf)
sample <- st_read("testdata.shp")
sample$area <- st_area(sample$geometry)
density_calc <-sample %>% st_buffer(0) %>% group_by(areas) %>% summarise(`Sales (density)` = sum(sales)/sum(area))
这里是 shapefile 的详细信息:
Geometry set for 2106 features
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -120.0065 ymin: 35.00184 xmax: -114.0396 ymax: 42.00221
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
我想我的问题是,我真的不知道什么是对什么是错,所以我不知道我是否做对了。
抱歉,如果这不是最广泛的问题,我只是不太记得我的高中几何!
raster
包有助于使此计算变得非常简单,就像在 R
中使用 data.frame 一样:
library(raster)
list.files(workDir)
test_shp <- shapefile(file.path(workDir, 'testdata.shp'))
names(test_shp)
#[1] "distrct" "sbdstrc" "terrtry"
#[4] "region" "turf" "sales"
#[7] "leads" "cnvrsns" "areas"
sum(is.na(test_shp$sales)) #note that 346 polygons have no sales data
#get the area as square kilometers
test_shp$km2 <- area(test_shp) / 10000
#calc the sales density
test_shp$sales_density <- test_shp$sales / test_shp$km2
#calculate the 25th, 50th, and 75th percentile of all polygons
quartiles <- quantile(test_shp$sales_density, probs=c(0.25, 0.5, 0.75), na.rm=TRUE)
#plot the result, coloring by which percentile the sales density is for a given polygon
plot(test_shp, col=ifelse(is.na(test_shp$sales_density), 'gray', ifelse(test_shp$sales_density >= quartiles[3], 'dark green', ifelse(test_shp$sales_density >= quartiles[2], 'light green', ifelse(test_shp$sales_density >= quartiles[1], 'yellow', 'red')))), border='transparent') (eg. >75th, 50-75th, etc.)
#add the legend
legend('bottomleft', legend=c('Q4', 'Q3', 'Q2', 'Q1', 'No data'), pch=15, col=c('dark green', 'light green', 'yellow', 'red', 'gray'))
我正在尝试计算 shapefile 中的密度,但我很确定我做错了。这个想法是按密度找出哪些地理区域的销售额最高。
Here is a link to the file that I use (testdata.shp)
library(sf)
sample <- st_read("testdata.shp")
sample$area <- st_area(sample$geometry)
density_calc <-sample %>% st_buffer(0) %>% group_by(areas) %>% summarise(`Sales (density)` = sum(sales)/sum(area))
这里是 shapefile 的详细信息:
Geometry set for 2106 features
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -120.0065 ymin: 35.00184 xmax: -114.0396 ymax: 42.00221
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
我想我的问题是,我真的不知道什么是对什么是错,所以我不知道我是否做对了。
抱歉,如果这不是最广泛的问题,我只是不太记得我的高中几何!
raster
包有助于使此计算变得非常简单,就像在 R
中使用 data.frame 一样:
library(raster)
list.files(workDir)
test_shp <- shapefile(file.path(workDir, 'testdata.shp'))
names(test_shp)
#[1] "distrct" "sbdstrc" "terrtry"
#[4] "region" "turf" "sales"
#[7] "leads" "cnvrsns" "areas"
sum(is.na(test_shp$sales)) #note that 346 polygons have no sales data
#get the area as square kilometers
test_shp$km2 <- area(test_shp) / 10000
#calc the sales density
test_shp$sales_density <- test_shp$sales / test_shp$km2
#calculate the 25th, 50th, and 75th percentile of all polygons
quartiles <- quantile(test_shp$sales_density, probs=c(0.25, 0.5, 0.75), na.rm=TRUE)
#plot the result, coloring by which percentile the sales density is for a given polygon
plot(test_shp, col=ifelse(is.na(test_shp$sales_density), 'gray', ifelse(test_shp$sales_density >= quartiles[3], 'dark green', ifelse(test_shp$sales_density >= quartiles[2], 'light green', ifelse(test_shp$sales_density >= quartiles[1], 'yellow', 'red')))), border='transparent') (eg. >75th, 50-75th, etc.)
#add the legend
legend('bottomleft', legend=c('Q4', 'Q3', 'Q2', 'Q1', 'No data'), pch=15, col=c('dark green', 'light green', 'yellow', 'red', 'gray'))