我对计算全球平均降水量感兴趣。我的数据包含在R中的RasterStack对象中,并且是全球分布的降水量数据。但是,我只想隔离陆地区域以分别计算这些区域的平均值。我想做的是以某种方式隔离中心与土地重叠的那些网格单元,然后计算年平均值。我已经首先创建了一个栅格堆栈对象,称为“ RCP1pctCO2Mean”,其中包含感兴趣的平均值,但它是全局分布的数据。共有138层,每层代表一年。
因此,我想做的是以下事情:
我的栅格堆栈对象“ 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, ...
这第一步有意义吗?这是否保留原始栅格数据,但现在为多边形形式?如果是这样,那么是否可以以某种方式隔离一个全地形多边形,然后以某种方式指定哪些单元格被视为陆地,然后计算网格单元中每年的平均值?
我之前已经寻找过任何可以执行此操作的过程示例,但是我还没有遇到任何太具体的事情。
感谢您的时间,任何帮助将非常宝贵!
您无需将栅格转换为多边形。最好使用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是地理坐标系,像素距赤道越远,其实际代表的土地面积就越小。为了从技术上克服这一问题,应该对降水值进行加权以反映与赤道的距离。