2014-10-20 51 views
0

我曾经在R中使用过rts。我有一组栅格,我需要在一段时间内找到平均值。这是我迄今所做的:R中一般均值和子集的问题

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) 

我有子集化手段十二月至二月很难再产生一个平均的光栅。这是whatI已经能够做的子集每个月份:

###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)) 

这看起来像周围ANS心不是让我给我所需要的最终结果很长的路要走。请问最简单的选择是什么?

回答

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) 

不幸的是,RTS包中的R并不强劲做某些事情呢。