在R中使用由多个多边形组成的shapefile从NetCDF提取数据？

Question

我有一个NetCDF文件，其中包含了38年的降水数据。我使用我的区域的shapefile来提取我感兴趣的区域的数据。当我使用下面的代码时，我得到了与多边形内落在网格点对应的数据。然而，我希望根据shapefile中的每个多边形（本例中为12个），提取数据并使用多边形的名称保存数据框。以下是我的示例代码。

library(ncdf4)
library(rgdal)
library(raster)
library(maptools)
library(GISTools)

NC = brick("Daily_Pcp.nc")
Subbasin_SHP=readOGR("my_Shapefile.shp")
crs(NC)
crs(Subbasin_SHP)
SHP=spTransform(Subbasin_SHP, crs(NC))

Polygon_Names=my_shapefile$Subbasin # there are 12 polygon (subbasins) in my shapefile with specific name

PCP_Data=mask(NC, Subbasin_SHP)
DF=as.data.frame(PCP_Data, xy=TRUE)
DF_insidepoints=DF[complete.cases(DF),]
write.csv(DF_insidepoints, "DataForEntireShapefile.csv")

这是一个包含所有多边形的shapefile地图。我希望将落在特定多边形中的所有点数据保存在该多边形的名称下。总共应使用mask函数获得12个文件，其中每个文件都有对应于落在该多边形内的网格点的数据。

plot(Subbasin_SHP, col="gray", border="blue", axes=TRUE, pbg="white")
pointLabel(coordinates(Subbasin_SHP), labels = Subbasin_SHP$Subbasin) 

- CForClimate

---

有人知道这个问题的答案吗？一种解决方法是将shapefile拆分为不同的多边形，然后提取每个多边形的数据并保存，但这是完成此任务非常费时的方式（我用过这种方法）。一定有更有效率的解决方法。 - CForClimate

Answer 1

你可以使用raster包中的extract函数结合stack函数。

NC=raster::stack("Daily_Pcp.nc") # Stack raster just like brick you did
Subbasin_SHP=readOGR("my_Shapefile.shp") # Read shape file
df=raster::extract(NC,Subbasin_SHP,fun=mean, na.rm=T) #Extract dataframe

数据框中的其他更改由您决定。例如，行和列名称。

- Suwas Ghimire