栅格数据到多边形以计算 R 中的平均值

Raster data to polygons to compute averages in R

我对计算全球平均降水量很感兴趣。我的数据包含在 R 中的 RasterStack 对象中,是全球分布的降水数据。但是,我只想隔离陆地区域以分别计算它们的平均值。我想做的是以某种方式隔离那些中心与陆地重叠的网格单元,然后计算年均值。我已经首先创建了一个名为 "RCP1pctCO2Mean" 的栅格堆栈对象,它包含感兴趣的平均值但是是全局分布的数据。共有 138 层,每一层代表一年。

因此,我想做的是:

  1. 将栅格转换为多边形,但在转换时保留原始的 138 个栅格图层
  2. 指定哪些网格单元是陆地(即落在陆地上的那些网格单元 centers/coordinates)
  3. 计算这些陆地网格单元的年平均降水量。

我的栅格堆栈对象 "RCP1pctCO2Mean" 具有以下属性:

class       : RasterStack 
dimensions  : 64, 128, 8192, 138  (nrow, ncol, ncell, nlayers)
resolution  : 2.8125, 2.789327  (x, y)
extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin,  
ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,       
layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   
layer.12,   layer.13,   layer.14,   layer.15, ... 
min values  : 0.42964514, 0.43375653, 0.51749371, 0.50838983, 0.45366730,    
0.53099146, 0.49757186, 0.45697752, 0.41382199, 0.46082401, 0.45516687, 
0.51857087, 0.41005131, 0.45956529, 0.47497867, ... 
max values  :   96.30350,  104.08584,   88.92751,   97.49373,   89.57201,   
90.58570,   86.67651,   88.33519,   96.94720,  101.58247,   96.07792,   
93.21948,   99.59785,   94.26218,   90.62138, ...  

之前,我尝试通过指定经度和纬度范围来隔离特定区域以获得该区域的均值和中值,就像这样:

>expansion1<-expand.grid(103:120, 3:15) #This is a range of longitudes and then latitudes
>lonlataaa<-extract(RCP1pctCO2Mean, expansion1)
>Columnaaa<-colMeans(lonlataaa)

#Packages loaded        
library(raster)
library(maps)
library(maptools)
library(rasterVis)

但是,使用上述方法,太多的水会与陆地面积混合,如果我缩小陆地上的 latitude/longitude 范围,我可能会错过太多的土地,无法有意义地计算平均值。

因此,有了这个 RasterStack,是否可以创建一个条件来告诉 R,如果每个网格单元的 "center point" 或质心(每个网格单元中心代表一个特定的 latitude/longitude坐标)恰好落在陆地上,那么它会被视为陆地(即这将是 TRUE - 如果不是,则为 FALSE,或者可能以某种方式使用 0s 或 1s)?即使一个格子碰巧有水和陆地混合在一起,但格子的中心point/centroid在陆地上,那也会被认为是陆地。我也想对特定国家/地区执行此操作,以同样的方式计算全国平均水平。

我希望保留个人 138 layers/years,以便可以对所有相关网格单元格的所有第 1 年进行平均,然后对所有第 2 年、所有第 3 年、所有第 4 年等进行平均.(稍后创建时间序列)。我不确定这是否是正确的方法,但我首先做的是采用 "RCP1pctCO2Mean" RasterStack 并使用以下方法创建了一个 SpatialPolygonsDataframe:

trans <- raster::rasterToPolygon(RCP1pctCO2Mean)    
trans

class      : SpatialPolygonsDataFrame 
features    : 8192 
extent      : -181.4062, 178.5938, -89.25846, 89.25846  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
variables  : 138
names      :          layer.1,          layer.2,          layer.3,          layer.4,          layer.5,          layer.6,          layer.7,          layer.8,          layer.9,          layer.10,        layer.11,          layer.12,          layer.13,          layer.14,          layer.15, ... 
min values  : 0.429645141288708, 0.433756527757047, 0.517493712584313, 0.508389826053265, 0.453667300300907, 0.530991463885754,  0.4975718597839, 0.456977516231847, 0.413821991744321, 0.460824014230889, 0.45516687008843, 0.518570869929649, 0.410051312472821, 0.459565291388527, 0.474978673081429, ... 
max values  :  96.3034965348338,  104.085840813594,  88.9275127089197,  97.4937326695693,  89.5720053000712,  90.5857030396531, 86.6765123781949,  88.3351859796546,  96.947199473011,  101.582468961459, 96.0779212204739,  93.2194802269814,  99.5978503841538,  94.2621847475029,  90.6213755054263, ... 

第一步有意义吗?这是否保留了原始栅格数据,但现在是多边形形式?如果是这样,是否有可能以某种方式仅隔离一个全陆地多边形,然后以某种方式指定哪些单元格被视为陆地,然后计算网格单元格每年的平均值?

我以前曾寻找过任何程序示例来执行此操作,但我还没有遇到任何太具体的事情。

感谢您抽出宝贵时间,任何帮助都将非常宝贵!

您不需要将栅格转换为多边形。最好使用 shapefile 来掩盖您的栅格并仅保留陆地上的点。如果您的降水数据已经采用名为 "RCP1pctCO2Mean" 的 rasterStack 格式,那么您可以这样做:

require(maptools)
require(raster)
data(wrld_simpl) ##this loads a coarse-resolution global shapefile from the maptools package


##GLOBAL EXAMPLE
##isolate only points on land
##wrld_simpl and RCP1pctCO2Mean have the same CRS so can be used directly
all_land_vals = mask(RCP1pctCO2Mean,wrld_simpl)

##get mean value for all points but for each layer separately 
cellStats(all_land_vals, stat='mean')


##COUNTRY EXAMPLE
##isolate only points on Brazil
country=wrld_simpl[which(wrld_simpl$NAME=="Brazil"),]
country_vals = mask(RCP1pctCO2Mean,country)

##get mean value for all points in given country but for each layer separately 
cellStats(country_vals, stat='mean')

N.B。您的降水数据非常粗糙,因此您可能需要先 disaggregate 它。此外,您直接计算平均值的方法存在缺陷,因为 latlong 是一个地理坐标系,并且像素离赤道越远,它实际代表的土地面积就越小。为了克服这个技术上的降水值应该加权以反映距赤道的距离。