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

问题描述 投票:0回答:1

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

因此,我想做的是以下事情:

  1. [将栅格转换为多边形,但在转换时保留原始的138个栅格图层
  2. 指定哪些网格单元是陆地(即那些网格单元中心/落在陆地上的坐标)
  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)

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

因此,使用此RasterStack,可以创建一个条件来告诉R如果每个网格单元(每个网格单元中心代表一个特定的纬度/经度坐标)的“中心点”或质心碰巧落在土地,那么将其视为土地(即为TRUE-如果不是,则为FALSE,或者以某种方式使用0或1s)?即使网格单元恰好有水与陆地混合,但是网格的中心点/质心在陆地上,也可以视为陆地。我也想对特定国家/地区进行此操作,以相同的方式计算国家平均值。

我希望保留单个138层/年,以便可以在所有相关网格单元中平均所有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, ... 

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

我之前已经寻找过任何可以执行此操作的过程示例,但是我还没有遇到任何太具体的事情。

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

r r-raster
1个回答
0
投票

您无需将栅格转换为多边形。最好使用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是地理坐标系,像素距赤道越远,其实际代表的土地面积就越小。为了从技术上克服这一问题,应该对降水值进行加权以反映与赤道的距离。

© www.soinside.com 2019 - 2024. All rights reserved.