【问题标题】:How to overlay raster and point layers and calculate summary statistics如何叠加栅格和点图层并计算汇总统计数据
【发布时间】:2018-07-19 20:18:46
【问题描述】:

我正在寻找一些关于基本地理空间统计的建议。 我正在使用来自 Worldpop 的栅格文件,表示巴西每 100 平方米的人口。我有另一个带有巴西医院坐标的 latlong 数据集。

我想做以下事情:

  1. 识别距离医院 1 公里以内、2-10 公里和 10 公里以上的区域(并可能创建多边形)
  2. 计算以上每个区域的人数

我想提供一些可重现的示例,但光栅文件非常大。有一些instructions on how to do this 带有两个单独的纬度/经度点列表,但我不知道如何使用光栅文件来执行此操作。

有什么想法吗?

【问题讨论】:

    标签: r geospatial raster


    【解决方案1】:

    示例数据

    library(raster)
    bra <- getData('GADM', country="BRA", level=1)
    r <- raster(bra, res=1)
    values(r) <- 1:ncell(r)
    r <- mask(r, bra)
    pts <- coordinates(bra)
    # plot(r)
    # points(pts)
    

    解决方案

    b1 <- extract(r, pts, buffer=100000)  # 100 km
    b2 <- extract(r, pts, buffer=200000)  # 200 km
    pop1 <- sapply(b1, sum)
    pop2 <- sapply(b2, function(i)sum(i, na.rm=TRUE)) - pop1
    

    查看区域

    spts <- SpatialPoints(pts, proj4string=crs(bra))
    buf1 <- buffer(spts, width=100000, dissolve=FALSE)
    buf2 <- buffer(spts, width=200000, dissolve=FALSE)
    
    # adding IDs so that they can also be used in "extract"
    buf1 <- SpatialPolygonsDataFrame(buf1, data.frame(id1=1:length(buf1)))
    buf2 <- SpatialPolygonsDataFrame(buf2, data.frame(id2=1:length(buf2)))
    
    # To combine buf1 and buf2 you could do
    # buf <- (buf2-buf1) + buf1
    # but in this example there are overlapping buffers, so I do
    bb <- list()
    for (i in 1:length(buf1)) {
        bb[[i]] <- (buf2[i,]-buf1[i,]) + buf1[i,]
    }
    buf <- do.call(bind, bb)
    
    plot(r)
    plot(buf,  col=c("red", "blue"), add=TRUE)
    

    现在你可以做

    z <- extract(r, buf, fun=sum, na.rm=TRUE)
    z <- cbind(data.frame(buf), z)
    head(z)
    

    pop1 和 pop2 得到与上面相同的结果

    head(pop1)
    head(pop2)
    

    【讨论】:

    • 谢谢罗伯特H。几个问题:我有兴趣根据与最近医院的距离来识别地图中的不同区域。我想映射它并可能将它放在多边形中 - 例如1区10km),并计算这些区域中有多少人。这可能吗?
    • 我添加了另一个区域
    猜你喜欢
    • 1970-01-01
    • 2023-01-17
    • 1970-01-01
    • 2019-08-18
    • 1970-01-01
    • 2023-01-05
    • 1970-01-01
    • 1970-01-01
    • 2023-03-24
    相关资源
    最近更新 更多