代码之家 › 专栏 › 技术社区 › Rich Pauloo Yan Zhang

多边形的表面积加权空间连接

sf sp spatial r

2

Rich Pauloo Yan Zhang · 技术社区 · 7 年前

问题陈述

我有两组多边形,我想把定量特征从一组多边形连接到另一组。

例如,考虑Yolo County的Multipolygon, Yolo 。我要从字段中的所有功能中聚合tract-level数据 estimate that fit in to the polygon of the city of davis, davis 。


   
    
     
      结果应该是
      
       Davis的多边形,带有一个新字段
       
        Estimate
        
         That is the
         
          Surface Area Weighted
          
           Estimate
           
            of all the features in
            
             Yolo
             
              that fall within
              
               Davis.
               
                
                 。如何在
                 
                  sp或
                  
                   sf中执行此操作?
                  
                 
                
               
              
             
            
           
          
         
        
       
      
     
     
      
       
        
         
          
          
           重复性示例
          
          
           Davis Polygon城市(
           
            Davis
            
             )下载自
             
              this website,file:
              
               citylimits.zip
              
              。
             
            
           
          
          
           
            
             软件包
图书馆(Tidycenses)
图书馆(tidyverse)
图书馆(光栅)

#获取约罗县的土地水平数据
yolo<-获取acs(state=“ca”,county=“yolo”,geography=“tract”,
variables=“B19013 U 001”,geometry=true)


#戴维斯市形状文件
Davis<-光栅::shapefile(“Davis_dbo_citylimits.shp”)。
Davis<-Davis%>%spTransform(,St_CRS(Yolo)$`Proj4String`%>%CRS())
戴维斯<-St_as_SF(戴维斯)
yolo<-yolo%>%st_转换(st_CRS(Davis)$项目4字符串`)

γ图
GGPROTH()+
geom_sf(数据=yolo,aes(填充=估计))。+
geom_sf(data=Davis,alpha=0.3,color=“red”)。+
坐标SF(xlim=c(-121.6,-121.9),ylim=c(38.5,38.6))
< /代码>
            
            
             
             
             
              
             
            
            
           
           
            
             
              注意
             
             :我已经看到了
             
              this so post
             
             
              。死链接使它不可简化,它使用神秘的软件包。yolo
             
            
            
             . 我想从该领域的所有特性中汇总tract-level数据
             
              estimate
             
             在戴维斯市的多边形内,
             
              davis
             
             .
            
           
           
            
             结果应该是一个多边形
             
              戴维斯
             
             用新字段
             
              估计
             
             那就是
             
              表面积加权
              
               估计
              
              在所有功能中
              
               约洛
              
              属于
              
               戴维斯
              
             
             . 我该怎么做呢
             
              sp
             
             或
             
              sf
             
             ?
            
            
            
             重复性示例
            
           
           
            
             戴维斯多边形市(
             
              戴维斯
             
             )下载自
            
            
             this website, file:
             
              CityLimits.zip
             
            
            .
           
           # packages
library(tidycensus)
library(tidyverse)
library(raster)

# get tract level data for yolo county
yolo  <- get_acs(state = "CA", county = "Yolo", geography = "tract", 
                 variables = "B19013_001", geometry = TRUE)


# city of davis shapefile
davis <- raster::shapefile("Davis_DBO_CityLimits.shp")
davis <- davis %>% spTransform(., st_crs(yolo)$`proj4string` %>% crs())
davis <- st_as_sf(davis)
yolo <- yolo %>% st_transform(st_crs(davis)$`proj4string`) 

# plot
ggplot() + 
  geom_sf(data = yolo, aes(fill = estimate)) +
  geom_sf(data = davis, alpha = 0.3, color = "red") +
  coord_sf(xlim=c(-121.6, -121.9), ylim = c(38.5, 38.6))

           
            
             
            
           
           
           
           
           
            
             注释
            
            我见过
            
             this SO post
            
            . 死链接使它不可简化,它使用神秘的软件包。

1 回复 | 直到 7 年前

1

4

Robert Hijmans 7 年前

library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p$value <- 1:length(p)
b <- as(extent(6, 6.4, 49.76, 50), 'SpatialPolygons')
b <- SpatialPolygonsDataFrame(b, data.frame(bid = 1))
crs(b) <- crs(p)
plot(p)
plot(b, add=T, border='red', lwd=2)

i <- intersect(b, p)    
i$AREA <- area(i) / 1000000
aw <- sum(i$AREA * i$value) / sum(i$AREA)
aw
# 5.086891

sp

a <- aggregate(p['value'], b, FUN=sum, areaWeighted=TRUE)
a$value
# 5.085438

sf

library(sf)
pf <- as(p, 'sf')
bf <- as(b, 'sf')

x <- sf::st_interpolate_aw(pf['value'], bf, extensive=F)
x$value
# 5.086891