R：计算由地理点周围的非线性函数加权的多边形面积总和或平均栅格像素值
R: Calculate sum of polygon area or mean raster pixel value weighted by a non-linear function around geographic points
我有 100 个采样点分布在我的研究区域以及描述该研究区域环境背景的空间多边形和栅格。我有兴趣根据非线性函数计算每个采样点周围的面积（对于多边形）或加权平均值（对于栅格）的加权总和，该函数描述环境上下文变量的相对重要性如何随着您的进一步降低远离采样点
这是非线性方程，其中r是距离采样点的距离（km）：w(r)=0.0382131 × exp(-0.49r)
我没能找到以前在空间数据上尝试过的例子。在这一点上，我的策略是在每个采样点周围创建许多大小不断增加的缓冲区，量化每个多边形的面积或栅格的平均值，然后将上述函数应用于生成的数值向量。
这是一个玩具示例，我用上面提到的简单方法计算加权值：
library(raster)
library(mapview)
library(sp)
library(rgeos)
library(maptools)
#loading a raster from the raster package
filename <- system.file("external/test.grd", package="raster")
r <- raster(filename)
r<-projectRaster(r, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
#generating some coordinates
coords_df<-data.frame( long=c(5.75718, 5.74224, 5.73521),lat=c(50.98469, 50.97551, 50.96372))
xy<-SpatialPoints(coords_df, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
xy_df<-SpatialPointsDataFrame(coords_df,data=data.frame(ID=c(1,2,3), row.names=c(1,2,3)), match.ID=T, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
#creating some polygons
poly1_xcoord<-c(5.757752, 5.758310, 5.756508, 5.755950)
poly1_ycoord<-c(50.986693,50.985828,50.985072,50.986017)
poly2_xcoord<-c(5.739311,5.740770,5.741907, 5.740480)
poly2_ycoord<-c(50.976658,50.975942,50.976253,50.976929)
poly3_xcoord<-c(5.734416,5.734759,5.735425,5.735510)
poly3_ycoord<-c(50.964193,50.963706,50.963652,50.964315)
poly4_xcoord<-c(5.737270,5.738643,5.760530, 5.759328)
poly4_ycoord<-c(50.961017,50.960368,50.983555, 50.983231)
poly1_coords <- cbind(poly1_xcoord, poly1_ycoord)
poly1 = Polygon(poly1_coords, hole = F)
poly2_coords <- cbind(poly2_xcoord, poly2_ycoord)
poly2 = Polygon(poly2_coords, hole=F)
poly3_coords <- cbind(poly3_xcoord, poly3_ycoord)
poly3 = Polygon(poly3_coords, hole=F)
poly4_coords <- cbind(poly4_xcoord, poly4_ycoord)
poly4 = Polygon(poly4_coords, hole=F)
polys_spatial = SpatialPolygons(list(Polygons(list(poly1,poly2,poly3,poly4), 1)))
proj4string(polys_spatial) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
#mapview(r) +mapview(xy)+ mapview(polys_spatial)
#bad way of calculating weighted mean for raster
#extracting mean raster value within buffers around each point
raster_vals_100<-sapply(extract(r, xy, buffer=c(100)), mean)
raster_vals_200<-sapply(extract(r, xy, buffer=c(200)), mean)
raster_vals_300<-sapply(extract(r, xy, buffer=c(300)), mean)
#and so on..
raster_vals<-list(raster_vals_100, raster_vals_200, raster_vals_300)
#nonlinear function
weight<-function(x){0.0382131 * exp(-.49*x)}
#distances of buffers
dist<-c(.1,.2,.3)
#calculating weighted mean of raster values
site1_raster_mean<-mean(w(c(dist)) * sapply( raster_vals, function(x) x[1]))
site2_raster_mean<-mean(w(c(dist)) * sapply( raster_vals, function(x) x[2]))
site3_raster_mean<-mean(w(c(dist)) * sapply( raster_vals, function(x) x[3]))
#bad way of calculating weighted sum for polygons
#calculating weighted sum of area for polygons (I know the width argument isn't in proper units, but doesn't affect main question)
site_buffers_100<-gBuffer(xy_df, width=.001, byid=T)
site_buffers_200<-gBuffer(xy_df, width=.002, byid=T)
site_buffers_300<-gBuffer(xy_df, width=.003, byid=T)
#preventing weird orphaned hole issue
slot(polys_spatial, "polygons") <- lapply(slot(polys_spatial, "polygons"), checkPolygonsHoles)
#extracting intersection between polygons and buffers around sites
poly_intersect_100<-raster::intersect(site_buffers_100,polys_spatial)
poly_intersect_200<-raster::intersect(site_buffers_200, polys_spatial)
poly_intersect_300<-raster::intersect(site_buffers_300, polys_spatial)
#summing the area of the intersecting polygons
poly_intersect_100_area<-gArea(poly_intersect_100, byid = TRUE)
poly_intersect_200_area<-gArea(poly_intersect_200, byid = TRUE)
poly_intersect_300_area<-gArea(poly_intersect_300, byid = TRUE)
area_list<-list(poly_intersect_100_area,poly_intersect_200_area,poly_intersect_300_area)
#calculating the weighted sum by site
dist<-c(.1,.2,.3)
site1_polygon_sum<-sum(w(c(dist)) * sapply( area_list, function(x) x[1]))
site2_polygon_sum<-sum(w(c(dist)) * sapply( area_list, function(x) x[2]))
site3_polygon_sum<-sum(w(c(dist)) * sapply( area_list, function(x) x[3]))
我尝试使用连续的距离测量而不是缓冲区，尽管我不确定这对多边形是否有意义，因为本身没有 overlap。可能可以栅格化多边形来做到这一点。
总体思路是根据到多边形栅格像元的距离对栅格值进行加权。我使用 spex 包进行多边形化只是因为它是 fastest:
library(sf)
library(spex)
library(fasterize)
library(lwgeom)
library(mapview)
xy_sf <- st_as_sf(xy_df)
r_poly <- st_as_sf(spex::polygonize(r))
r_poly <- r_poly[,-1]
r_poly$vals <- r[]
r_dists <- st_distance(xy_sf, r_poly)
units(r_dists) <- "km"
r_dists <- data.frame(matrix(t(r_dists), ncol = 3))
r_dists <- data.frame(apply(r_dists, 2, function(x) w(x) * r_poly$vals))
r_dists <- dplyr::bind_cols(r_poly, r_dists)
par(mfrow = c(1, 3))
plot(fasterize(r_dists, r, field = "X1"))
plot(fasterize(r_dists, r, field = "X2"))
plot(fasterize(r_dists, r, field = "X3"))
par(mfrow = c(1,1))
polys_sf <- st_cast(st_as_sf(polys_spatial), "POLYGON")
polys_sf$vals <- st_area(polys_sf)
polys_dists <- matrix(st_distance(xy_sf, polys_sf), nrow = 3)
polys_dists <- w(polys_dists)
xy_sf$polys_weights <- rev(colSums(polys_sf$vals * t(polys_dists)))
xy_sf$raster_weights <- colSums(st_set_geometry(r_dists[,c("X1", "X2", "X3")], NULL), na.rm = TRUE)
mapview(polys_sf) + mapview(xy_sf)
我有 100 个采样点分布在我的研究区域以及描述该研究区域环境背景的空间多边形和栅格。我有兴趣根据非线性函数计算每个采样点周围的面积（对于多边形）或加权平均值（对于栅格）的加权总和，该函数描述环境上下文变量的相对重要性如何随着您的进一步降低远离采样点
这是非线性方程，其中r是距离采样点的距离（km）：w(r)=0.0382131 × exp(-0.49r)
我没能找到以前在空间数据上尝试过的例子。在这一点上，我的策略是在每个采样点周围创建许多大小不断增加的缓冲区，量化每个多边形的面积或栅格的平均值，然后将上述函数应用于生成的数值向量。
这是一个玩具示例，我用上面提到的简单方法计算加权值：
library(raster)
library(mapview)
library(sp)
library(rgeos)
library(maptools)
#loading a raster from the raster package
filename <- system.file("external/test.grd", package="raster")
r <- raster(filename)
r<-projectRaster(r, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
#generating some coordinates
coords_df<-data.frame( long=c(5.75718, 5.74224, 5.73521),lat=c(50.98469, 50.97551, 50.96372))
xy<-SpatialPoints(coords_df, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
xy_df<-SpatialPointsDataFrame(coords_df,data=data.frame(ID=c(1,2,3), row.names=c(1,2,3)), match.ID=T, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
#creating some polygons
poly1_xcoord<-c(5.757752, 5.758310, 5.756508, 5.755950)
poly1_ycoord<-c(50.986693,50.985828,50.985072,50.986017)
poly2_xcoord<-c(5.739311,5.740770,5.741907, 5.740480)
poly2_ycoord<-c(50.976658,50.975942,50.976253,50.976929)
poly3_xcoord<-c(5.734416,5.734759,5.735425,5.735510)
poly3_ycoord<-c(50.964193,50.963706,50.963652,50.964315)
poly4_xcoord<-c(5.737270,5.738643,5.760530, 5.759328)
poly4_ycoord<-c(50.961017,50.960368,50.983555, 50.983231)
poly1_coords <- cbind(poly1_xcoord, poly1_ycoord)
poly1 = Polygon(poly1_coords, hole = F)
poly2_coords <- cbind(poly2_xcoord, poly2_ycoord)
poly2 = Polygon(poly2_coords, hole=F)
poly3_coords <- cbind(poly3_xcoord, poly3_ycoord)
poly3 = Polygon(poly3_coords, hole=F)
poly4_coords <- cbind(poly4_xcoord, poly4_ycoord)
poly4 = Polygon(poly4_coords, hole=F)
polys_spatial = SpatialPolygons(list(Polygons(list(poly1,poly2,poly3,poly4), 1)))
proj4string(polys_spatial) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
#mapview(r) +mapview(xy)+ mapview(polys_spatial)
#bad way of calculating weighted mean for raster
#extracting mean raster value within buffers around each point
raster_vals_100<-sapply(extract(r, xy, buffer=c(100)), mean)
raster_vals_200<-sapply(extract(r, xy, buffer=c(200)), mean)
raster_vals_300<-sapply(extract(r, xy, buffer=c(300)), mean)
#and so on..
raster_vals<-list(raster_vals_100, raster_vals_200, raster_vals_300)
#nonlinear function
weight<-function(x){0.0382131 * exp(-.49*x)}
#distances of buffers
dist<-c(.1,.2,.3)
#calculating weighted mean of raster values
site1_raster_mean<-mean(w(c(dist)) * sapply( raster_vals, function(x) x[1]))
site2_raster_mean<-mean(w(c(dist)) * sapply( raster_vals, function(x) x[2]))
site3_raster_mean<-mean(w(c(dist)) * sapply( raster_vals, function(x) x[3]))
#bad way of calculating weighted sum for polygons
#calculating weighted sum of area for polygons (I know the width argument isn't in proper units, but doesn't affect main question)
site_buffers_100<-gBuffer(xy_df, width=.001, byid=T)
site_buffers_200<-gBuffer(xy_df, width=.002, byid=T)
site_buffers_300<-gBuffer(xy_df, width=.003, byid=T)
#preventing weird orphaned hole issue
slot(polys_spatial, "polygons") <- lapply(slot(polys_spatial, "polygons"), checkPolygonsHoles)
#extracting intersection between polygons and buffers around sites
poly_intersect_100<-raster::intersect(site_buffers_100,polys_spatial)
poly_intersect_200<-raster::intersect(site_buffers_200, polys_spatial)
poly_intersect_300<-raster::intersect(site_buffers_300, polys_spatial)
#summing the area of the intersecting polygons
poly_intersect_100_area<-gArea(poly_intersect_100, byid = TRUE)
poly_intersect_200_area<-gArea(poly_intersect_200, byid = TRUE)
poly_intersect_300_area<-gArea(poly_intersect_300, byid = TRUE)
area_list<-list(poly_intersect_100_area,poly_intersect_200_area,poly_intersect_300_area)
#calculating the weighted sum by site
dist<-c(.1,.2,.3)
site1_polygon_sum<-sum(w(c(dist)) * sapply( area_list, function(x) x[1]))
site2_polygon_sum<-sum(w(c(dist)) * sapply( area_list, function(x) x[2]))
site3_polygon_sum<-sum(w(c(dist)) * sapply( area_list, function(x) x[3]))
我尝试使用连续的距离测量而不是缓冲区，尽管我不确定这对多边形是否有意义，因为本身没有 overlap。可能可以栅格化多边形来做到这一点。
总体思路是根据到多边形栅格像元的距离对栅格值进行加权。我使用 spex 包进行多边形化只是因为它是 fastest:
library(sf)
library(spex)
library(fasterize)
library(lwgeom)
library(mapview)
xy_sf <- st_as_sf(xy_df)
r_poly <- st_as_sf(spex::polygonize(r))
r_poly <- r_poly[,-1]
r_poly$vals <- r[]
r_dists <- st_distance(xy_sf, r_poly)
units(r_dists) <- "km"
r_dists <- data.frame(matrix(t(r_dists), ncol = 3))
r_dists <- data.frame(apply(r_dists, 2, function(x) w(x) * r_poly$vals))
r_dists <- dplyr::bind_cols(r_poly, r_dists)
par(mfrow = c(1, 3))
plot(fasterize(r_dists, r, field = "X1"))
plot(fasterize(r_dists, r, field = "X2"))
plot(fasterize(r_dists, r, field = "X3"))
par(mfrow = c(1,1))
polys_sf <- st_cast(st_as_sf(polys_spatial), "POLYGON")
polys_sf$vals <- st_area(polys_sf)
polys_dists <- matrix(st_distance(xy_sf, polys_sf), nrow = 3)
polys_dists <- w(polys_dists)
xy_sf$polys_weights <- rev(colSums(polys_sf$vals * t(polys_dists)))
xy_sf$raster_weights <- colSums(st_set_geometry(r_dists[,c("X1", "X2", "X3")], NULL), na.rm = TRUE)
mapview(polys_sf) + mapview(xy_sf)