R语言 如何提取静止多个多边形

xvw2m8pv  于 2023-09-27  发布在  其他
关注(0)|答案(1)|浏览(73)

使用下面的代码来给予每个地区的平均值,这个链接CHIRPS tar.gz data用于降水数据,library(geodata)用于行政区。

library(sf)
library(raster)
# library(geodata) for gadm shape files
rwanda_shapefile <- st_read("C:/Users/HP/Desktop/IITA/saved_OUTPUT/CHIRPS/PREP/gadm41_RWA_shp/gadm41_RWA_0.shp")

precipitation_data <- raster("C:/Users/HP/Desktop/IITA/saved_OUTPUT/CHIRPS/PREP/chirps-v2.0.2021.08.tif")

district_precipitation <- extract(precipitation_data, rwanda_shapefile)
  
d_mean_prec <- sapply(district_precipitation, mean, na.rm = TRUE)

d_prec$Mean_Prec <- d_mean_prec

print(d_mean_prec)

得到这些地区的平均降水量

mbjcgjjk

mbjcgjjk1#

如果您查看?raster::extract,它的用法是

extract(x, y, fun=NULL, na.rm=FALSE, 
  exact=FALSE, weights=FALSE,
  normalizeWeights=TRUE, cellnumbers=FALSE,
  small=TRUE, df=FALSE, layer, nl, 
  factors=FALSE, sp=FALSE, ...)

所以你可以在extract()调用中得到分区的平均值。

library(sf)
library(raster)
library(terra)
library(geodata)

Rwanda_precincts1 = geodata::gadm('RWA', level = 1, path = 'maps')# maps a dir in working directory
precip_wld = terra::rast(R.utils::gunzip('~/Downloads/chirps-v2.0.2021.08.tif.gz', remove = FALSE)) 
raster::extract(precip_wld, Rwanda_precincts1, fun = mean, na.rm = TRUE)
  ID chirps-v2.0.2021.08
1  1            47.29495
2  2            33.24807
3  3            25.63329
4  4            47.92463
5  5            28.60866

所以上面的raster::extract返回一个data.frame,其中包含按分区号划分的平均降水量。如果您想更新值,使其仅按分区表示精确度,则可以使用以下方法- hearafter using terra

# clip out the precip from world to just Rwanda
precip_RWA = (crop(precip_wld, Rwanda_precincts1, mask = TRUE)) 
# mask each precinct
RWA_pre1_msk = mask(precip_RWA, Rwanda_precincts1[1])
RWA_pre2_msk = mask(precip_RWA, Rwanda_precincts1[2])
RWA_pre3_msk = mask(precip_RWA, Rwanda_precincts1[3])
RWA_pre4_msk = mask(precip_RWA, Rwanda_precincts1[4])
RWA_pre5_msk = mask(precip_RWA, Rwanda_precincts1[5])
# get extract( values using terra
precip_RWA_df = extract(precip_RWA, Rwanda_precincts1, fun = 'mean', na.rm = TRUE)
# extract from df to vector
precip_vect1 = precip_RWA_df[, 2]
# change RWA_pre1:pre5 cell values to precip_vect1[1:5] values
RWA_pre1_msk[not.na(RWA_pre1_msk)] = precip_vect1[1]
RWA_pre2_msk[not.na(RWA_pre2_msk)] = precip_vect1[2]
RWA_pre3_msk[not.na(RWA_pre3_msk)] = precip_vect1[3]
RWA_pre4_msk[not.na(RWA_pre4_msk)] = precip_vect1[4]
RWA_pre5_msk[not.na(RWA_pre5_msk)] = precip_vect1[5]
# stitch these together with `terra::mosaic(
RWA_precip_mosaic = mosaic(RWA_pre1_msk, RWA_pre2_msk, RWA_pre3_msk, RWA_pre4_msk, RWA_pre5_msk, fun = 'min')
plot(RWA_precip_mosaic)

可能还有其他更好的方法来做到这一点。

相关问题