0

我有一段时间在 R 中使用 rts。我有一组栅格,我需要找到一段时间内的平均值。这是我到目前为止所做的:

fun1G<-function(x){
 attr1<-"HDFEOS/GRIDS/ColumnAmountNO2/Data Fields/ColumnAmountNO2Trop"
 y<-h5read(x,attr1)
 y[y<=0]=NA
 z<-t(y)
 R1<-raster (z, xmn=-180,xmx=180,ymn=-90,ymx=90)
 R2<-flip (R1,direction="y")
 proj4string(R2)<-CRS("+proj=utm +ellps=WGS84")
 myproj = "+proj=utm +north +zone=28, 29, 30,31, 32, 33 +ellps=WGS84" 
 WA= readShapeSpatial("West_Africa.shp", proj4string = CRS(myproj))
  WAN<- as(WA, 'SpatialPolygons')
 frm <- extent(c(-19, 19,2,29))
 pfrm <- as(frm, 'SpatialPolygons')

  #Create Buffer around west africa (WAN) to crop to good extents
  BUFF2<-gBuffer(WAN,width=1)
  WAE<-crop(R2,pfrm)}

 regexp2<-"([[:digit:]]{4})([[:alpha:]]{1})([[:digit:]]{4})"
##To extract the dates
 DatesOMI<-lapply(OMI, function(x)stri_extract_first_regex(x,regexp2))

 ##To remove "m" from the dates
 DateO2<-gsub("m", "", DatesOMI)
 DDO2<-data.frame(DateO2)
 DDO3<-DDO2[["DateO2"]]

 #COnvert straight to dates
 Timex<-as.Date(DDO3, "%Y%m%d")

 #### TO create the rts(raster time series) and subset
  RNOM2<-stack(CM2)
  NOrts<-rts(RNOM2,Timex)
  RN0411<-NOrts[["2004-10-30/2011-07-10"]]

  ##TO get the monthly data
  RNMM<-apply.monthly(RN0411,mean,na.rm=T)

我很难将 12 月到 2 月的均值子集化,然后生成平均栅格。这就是我在每个月中能够对子集进行的操作:

 ###Subset for december, January and february
 Decsub<-subset(RNMM, seq(from=3, to=82, by=12))
 Jansub<-subset(RNMM, seq(from=4, to=82, by=12))
 Febsub<-subset(RNMM, seq(from=5, to=82, by=12))

这看起来还有很长的路要走,并没有让我得到我需要的最终结果。请问有哪些最简单的方法可以做到这一点?

4

1 回答 1

0

我最终使用的一种方法是:

MMD<-period.apply(Decsub,FUN=mean,na.rm=T,INDEX=...)
MMJ<-period.apply(Jansub,FUN=mean,na.rm=T,INDEX=...)

其中 INDEX 是每个子集中的栅格总数。写入文件

write.rts(MMD,"Decsub",overwrite=T)
write.rts(CMMW,"Jansub",overwrite=T)
write.rts(CMMW,"Febsub",overwrite=T)

R 将这些文件写入 3 个文件的单独文件夹中,因此您需要将这些文件复制到您的主工作文件夹中。作为单独的栅格带回

D<-raster("Deccub.gri")
J<-raster("Jansub.gri")
F<-raster("Febsub.gri")

堆叠栅格

DJF<-stack(D,J,F)
DJFmean <- calc(DJF,mean)

不幸的是,R 中的 rts 包在做某些事情时还不够健壮。

于 2015-01-16T13:55:02.507 回答