Project

General

Profile

« Previous | Next » 

Revision b74bf08c

Added by Adam Wilson almost 11 years ago

adding uncorrected MCD output in wgs84

View differences:

climate/research/cloud/MOD09/MOD09_Visualize.R
6 6

  
7 7
## read in global coasts for nice plotting
8 8
library(maptools)
9
library(doMC)
9 10

  
11
registerDoMC(6)
12
#### Evaluate MOD35 Cloud data
13
#mc=brick("~/acrobates/adamw/projects/cloud/data/mod09.nc",varname="CF")
14
f=list.files("/mnt/data2/projects/cloud/mcd09tif",pattern=paste(".*[O].*[.]tif$",sep=""),full=T)
15
f=f[c(4:12,1:3)]   
16
mc=stack(f,bands=1)
17
names(mc)=month.name
18

  
19
#NAvalue(mc)=-1
10 20
#coast=spTransform(coast,CRS(projection(mod35)))
11 21
land=readShapePoly("/mnt/data/jetzlab/Data/environ/global/gshhg/GSHHS_shp/c/GSHHS_c_L1.shp",force_ring=TRUE)
12 22
projection(land)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
13
CP <- as(extent(-180, 180, -60, 84), "SpatialPolygons")
23
CP <- as(extent(-180, 180, -60, 66), "SpatialPolygons")
14 24
proj4string(CP) <- CRS(proj4string(land))
15
coast=as(land[land$area>50,],"SpatialLines")
16

  
17

  
18
#### Evaluate MOD35 Cloud data
19
mc=brick("~/acrobates/adamw/projects/cloud/data/mod09.nc",varname="CF")
20
NAvalue(mc)=-1
21 25

  
22
cols=colorRampPalette(c("#000000","#00FF00","#FF0000"))#"black","blue","red"))
23
for(i in 1:156){
24
png(paste("output/mod09_fullanimation_",i,".png",sep=""),width=2000,height=1000)
25
  print(i)
26
  r=mm[[i]]
27
  print(levelplot(r,col.regions=cols(100),at=seq(0,100,len=100),margin=F,maxpixels=1e6,ylim=c(-60,70),main=paste(names(mod09)[i])))+
28
    layer(sp.lines(coast))
29
dev.off()
30
}
26
## World outline
27
world=readOGR("/mnt/data/jetzlab/Data/environ/global/MODIS/modis_sinusoidal/","worldborder_sinusoidal")
31 28

  
32
#### Evaluate MOD35 Cloud data
33
mmc=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonmean.nc",varname="CF")
34
names(mmc)=month.name
29
coast=crop(as(land[land$area>50,],"SpatialLines"),CP)
30
coast2=spTransform(coast,CRS(projection(mc)))
31
CP2=spTransform(CP,CRS(projection(mc)))
35 32

  
36 33
cols=colorRampPalette(c("#000000","#00FF00","#FF0000"))#"black","blue","red"))
37
png("output/CF_Animation_%03d.png",width=5000,height=4000,res=600,pointsize=96,bg="white")
38
for(i in 1:12){
34
foreach(i=1:12)%dopar% {
39 35
    print(i)
40
    r=mmc[[i]]
41
    print(levelplot(r,col.regions=cols(100),at=seq(1,100,len=100),margin=F,maxpixels=1e6,ylim=c(-60,73),
36
    r=mc[[i]]
37
    d1=(extent(coast2)@ymax-extent(coast2)@ymin)/(extent(coast2)@xmax-extent(coast2)@xmin)  # get y:x ratio to set output
38
    png(paste("output/CF_mean_",sprintf("%02d", i),".png",sep=""),width=1920,height=round(1920*d1),res=300,pointsize=46,bg="white")
39
    print(levelplot(r,col.regions=cols(100),at=seq(1,100,len=100),margin=F,maxpixels=2.1e6,ylim=c(extent(coast2)@ymin,extent(coast2)@ymax),
42 40
                    main=paste(month.name[i]),cex.main=3,scales=list(draw=F),cuts=99,ylab="",xlab="")+
43
                        layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+
44
                        layer(sp.lines(coast,col="black"),under=F))
41
                        layer(panel.polygon(c(extent(world)@xmin,extent(world)@xmin,extent(world)@xmax,extent(world)@xmax),y=c(extent(world)@ymin,extent(world)@ymax,extent(world)@ymax,extent(world)@ymin),col="black"),under=T)+
42
                        layer(sp.lines(coast2,col="black"),under=F)+
43
                        layer(sp.lines(world,col="black"),under=F))
44

  
45
    dev.off()
45 46
}
46
dev.off()
47 47

  
48 48

  
49

  
50

  
51
####################################################################
52
### Regional Comparisons
53
## Compare with worldclim and NPP
54
wc=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/prec_",1:12,".bil",sep="")))
55
#wc_map=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/bio_12.bil",sep="")))
56

  
57
regs=list(
58
  Cascades=extent(c(-122.8,-118,44.9,47)),
59
  Hawaii=extent(c(-156.5,-154,18.75,20.5)),
60
  Boliva=extent(c(-71,-63,-20,-15)),
61
  Venezuela=extent(c(-69,-59,0,7)),
62
  CFR=extent(c(17.75,22.5,-34.8,-32.6)),
63
  Madagascar=extent(c(46,52,-17,-12))
64
  #reg2=extent(c(-81,-70,-4,10))
65
  )
66
# convert to sinusoidal
67
regs2=lapply(regs,function(r){
68
       r2 = as(r, "SpatialPolygons")
69
       proj4string(r2) <- CRS(proj4string(land))
70
       r3=spTransform(r2,CRS(projection(mc)))
71
       return(r3)
72
   })
73

  
74

  
75
## read in GEWEX 1-degree data
76
gewex=brick("data/gewex/CA_PATMOSX_NOAA.nc",varname="a_CA")
77
names(gewex)=month.name#"1-degree Cloud Frequency (PATMOS-x GEWEX AVHRR)"
78

  
79

  
80

  
81

  
82

  
83

  
84
r="Venezuela"
85

  
86
## print global map with box for region
87
print(levelplot(mc[[1]],col.regions=cols(100),at=seq(1,100,len=100),margin=F,maxpixels=2.1e6,ylim=c(extent(coast2)@ymin,extent(coast2)@ymax),
88
                main=paste(month.name[i]),cex.main=3,scales=list(draw=F),cuts=99,ylab="",xlab="")+
89
      layer(panel.polygon(c(extent(world)@xmin,extent(world)@xmin,extent(world)@xmax,extent(world)@xmax),y=c(extent(world)@ymin,extent(world)@ymax,extent(world)@ymax,extent(world)@ymin),col="black"),under=T)+
90
      layer(sp.lines(coast2,col="black"),under=F)+
91
      layer(sp.lines(world,col="black"),under=F)+
92
      layer(sp.lines(regs[[r]],col="blue",lwd=2),under=F))
93

  
94

  
95
                                        # ylab.right = "Cloud Frequency (%)",par.settings = list(layout.widths = list(axis.key.padding = 0.1,axis.left=0.6,ylab.right = 3,right.padding=2)),
96

  
97
    ## crop the data
98
    tgewex=crop(gewex,regs[[r]])
99
    tmap=crop(wc,regs[[r]])
100
    tmc=projectRaster(crop(mc,regs2[[r]]),tmap)
101

  
102
## set plotting parameters
103

  
104
    range(c(cellStats(tmc,range)),cellStats(tgewex,range)*100)
105
range(cellStats(tmap,range))
106

  
107
pcols=colorRampPalette(c("#CCFFFF","#000066","#FF3300"))#"black","blue","red"))
108

  
109

  
110
foreach(i=1:12)%dopar% {
111
    print(i)
112
    png(paste("output/CF_mean_",r,"_",sprintf("%02d", i),".png",sep=""),width=1920,height=1080,res=300,pointsize=46,bg="white")
113
    pars=list(layout.heights=list(key.bottom=2,key.top=1),layout.widths = list(axis.key.padding = 4,axis.left=1))
114
    my.theme = trellis.par.get()
115
    my.theme$strip.background=list(col="transparent")
116
    trellis.par.set(my.theme)
117
    p1=levelplot(tmc[[i]],col.regions=cols(100),at=seq(20,100,len=99),
118
        colorkey=list(space="bottom",width=1,height=1.2,labels=list(labels=c(25,50,75,100),at=c(25,50,75,100)),just="left"),
119
        cuts=99,margin=F,max.pixels=1e6,par.settings = pars,main=paste(month.name[i]),cex.main=3,xlab="",ylab="",scales=list(x=list(at=c(-68,-64,-60))))
120
    p2=levelplot(tgewex[[i]],col.regions=cols(100),at=seq(.20,1,len=99),cuts=99,margin=F,max.pixels=1e6,
121
        colorkey=F,#list(space="bottom",width=.01,height=0.00000001,labels=list(labels="",at="",cex=0.0001)),
122
        par.settings = pars)
123
    p3=levelplot(tmap[[i]],col.regions=pcols(100),cuts=100,at=seq(0,650,len=100),margin=F,maxpixels=1e6,
124
        colorkey=list(space="bottom",height=.5,width=1,labels=list(labels=c(0,300,600),at=c(0,300,600))),xlab="",ylab="",main=names(regs)[r],useRaster=T,
125
        par.settings = pars)
126
    p_all=c("MODIS Cloud (%)"=p1,"PATMOS-x Cloud (%)"=p2,"WorldClim Precip (mm)"=p3,x.same=T,y.same=T,merge.legends=T,layout=c(3,1))
127
    print(p_all)
128
    dev.off()
129
}
climate/research/cloud/MOD09/ee/ee_compile.R
5 5
library(multicore)
6 6
library(foreach)
7 7
library(mgcv)
8
registerDoMC(2)
8
registerDoMC(12)
9 9

  
10 10

  
11 11
## start raster cluster
12
#beginCluster(10)
12
#beginCluster(5)
13 13

  
14 14
setwd("~/acrobates/adamw/projects/cloud")
15 15

  
......
38 38

  
39 39

  
40 40
jobs=expand.grid(month=1:12,sensor=c("MOD09","MYD09"))
41
## Loop over existing months to build composite netcdf files
41
i=1
42

  
43
#jobs=jobs[jobs$sensor=="MYD09",]
44

  
45

  
46
## Loop over data to mosaic tifs and remove problematic data at poles 
42 47
    foreach(i=1:nrow(jobs)) %dopar% {
43 48
        ## get month
44 49
        m=jobs$month[i]
......
47 52
        ## get sensor
48 53
        s=jobs$sensor[i]
49 54
        ## Define output and check if it already exists
50
        tvrt=paste(tmpfs,"/",s,"_",m,".vrt",sep="")
51
        ttif1=paste(tmpfs,"/",s,"_",m,".tif",sep="")
52
        ttif2=paste(tmpfs,"/",s,"_",m,"_wgs84.tif",sep="")
53
        ncfile=paste("data/mcd09nc/",s,"_",m,".nc",sep="")
55
        tvrt=paste(tmpfs,"/",s,"_",sprintf("%02d", m),".vrt",sep="")
56
        ttif1=paste("data/mcd09tif/",s,"_",sprintf("%02d", m),".tif",sep="")
54 57

  
55 58
        ## check if output already exists
56
        if(!rerun&file.exists(ncfile)) return(NA)
59
        if(!rerun&file.exists(ttif1)) return(NA)
57 60
        ## build VRT to merge tiles
58 61
        system(paste("gdalbuildvrt ",tvrt," ",paste(df$path[df$month==m&df$sensor==s],collapse=" ")))
59
        
60
        ## 
62

  
63
       ## 
61 64
        mod=stack(tvrt)
62 65
        names(mod)=c("cf","cfsd","nobs","pobs")
63
        #################################
64
        ## correct for orbital artifacts
65 66

  
66
        ## sample points to fit correction model
67
        reg=extent(-20016036,20016036, -4295789.46149, 4134293.61707)  #banding region
68
        cmod=crop(mod,reg)
69
        names(cmod)=c("cf","cfsd","nobs","pobs")
67
        ## mask no data regions (with less than 1 observation per day within that month)
68
        biasf=function(cf,cfsd,nobs,pobs) {
69
            ## drop data in areass with nobs<1 or pobs<=50 or where sd=0 and cf=0 or 50 (some polar regions)
70
            drop=nobs<=0|pobs<=50|(cf==0&cfsd==0)|(cf==50&cfsd==0)
71
            cf[drop]=NA
72
            cfsd[drop]=NA
73
            return(c(cf=cf,cfsd=cfsd,nobs=nobs,pobs=pobs))}
74

  
75
                
76
        mod2=overlay(mod,fun=biasf,unstack=TRUE,filename=ttif1,format="GTiff",
77
            dataType="INT1U",overwrite=T,NAflag=255, options=c("COMPRESS=LZW","BIGTIFF=YES"))
78

  
79
        tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS ",s,"GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).",
80
            "The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ",
81
            "Band Descriptions: 1) Mean Monthly Cloud Frequency 2) Four Times the Standard Deviation of Mean Monthly Cloud Frequency",
82
            " 3) Mean number of daily observations for each pixel 4) Proportion of days with at least one observation '"),
83
              "TIFFTAG_DOCUMENTNAME='Collection 5 ",s," Cloud Frequency'",
84
              paste("TIFFTAG_DATETIME='2013",sprintf("%02d", m),"15'",sep=""),
85
              "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")
86
        system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",ttif1," ",paste("-mo ",tags,sep="",collapse=" "),sep=""))
87
        
88
        rm(mod,mod2)
89
        writeLines(paste("Month:",m," Sensor:",s," Finished"))
90

  
91
    }
92

  
93

  
94
## Create combined (MOD+MYD) uncorrected mean CF
95
#    foreach(i=1:12) %dopar% {
96
#        ## get files
97
#        f=list.files("/mnt/data2/projects/cloud/mcd09tif",pattern=paste(".*[O|Y].*_",i,"[.]tif$",sep=""),full=T)
98
#        ## Define output and check if it already exists
99
#        tmcd=paste("/mnt/data2/projects/cloud/mcd09tif/MCD09_",sprintf("%02d", i),".tif",sep="")
100
#        ## check if output already exists
101
#        ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata 255 -dstnodata 255 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",
102
#            "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5")
103
#        system(paste("gdalwarp -overwrite -r average -co COMPRESS=LZW -co ZLEVEL=9  ",ops," ",paste(f,collapse=" ")," ",tmcd))
104
#        ## update metadata
105
#        tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS MOD09GA and MYD09GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).",
106
#            "The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ",
107
#            "Band Descriptions: 1) Mean Monthly Cloud Frequency 2) Four Times the Standard Deviation of Mean Monthly Cloud Frequency",
108
#            " 3) Mean number of daily observations for each pixel 4) Proportion of days with at least one observation '"),
109
#            "TIFFTAG_DOCUMENTNAME='Collection 5 MCD09 Cloud Frequency'",
110
#            paste("TIFFTAG_DATETIME='2013",sprintf("%02d", i),"15'",sep=""),
111
#              "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")
112
#        system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",tmcd," ",paste("-mo ",tags,sep="",collapse=" "),sep=""))
113
#        writeLines(paste("Finished month",i))
114
#    }
115

  
116

  
117
### Perform bias correction
118
foreach(i=1:nrow(jobs)) %dopar% {
119
        ## get month
120
        m=jobs$month[i]
121
        date=df$date[df$month==m][1]
122
        print(date)
123

  
124
        ttif2=paste(tmpfs,"/",s,"_",m,"_wgs84.tif",sep="")
125
        ncfile=paste("data/mcd09nc/",s,"_",m,".nc",sep="")
126

  
70 127
        modpts=sampleRandom(cmod, size=10000, na.rm=TRUE, xy=T, sp=T)
128
        rm(cmod)  #remove temporary raster to save space
71 129
        modpts=modpts[modpts$nobs>0,]  #drop data from missing tiles
72 130
        ### fit popbs correctionmodel (accounting for spatial variation in cf)
73 131
        modlm1=bam(cf~s(x,y)+pobs,data=modpts@data)
74 132
        summary(modlm1)
75 133
        modbeta1=coef(modlm1)["pobs"]
76 134

  
77

  
78
        #modlm2=bam(cf~s(x,y),data=modpts@data)
79
        #resid=residuals(modlm2)#,newdata=data.frame(x=modpts$x,y=modpts$y,pobs=100))
80
        #pred=predict(modlm1,newdata=data.frame(x=modpts$x,y=modpts$y,pobs=100))
81
        #plot(resid~pobs,data=modpts@data);abline(lm(resid~modpts$pobs))
82
        #plot(modlm1);points(modpts)
83

  
84
        #writeLines(paste(date,"       slope:",round(modbeta1,4)))
135
        writeLines(paste(date,"       slope:",round(modbeta1,4)))
136
        ## Smooth data to remove large-grain variability
137
#        system(paste("pkfilter -f mean -dx 99 -dy 99 -ot Byte -i ",df$path[df$month==m&df$sensor==s][1]," -o ",ttif1))
85 138

  
86 139
        ## mask no data regions (with less than 1 observation per day within that month)
87 140
        ## use model above to correct for orbital artifacts
......
101 154

  
102 155
        mod2=overlay(cmod,fun=biasf,unstack=TRUE,filename=ttif1,format="GTiff",
103 156
            dataType="INT1U",overwrite=T,NAflag=255, options=c("COMPRESS=LZW", "BIGTIFF=YES"))
104
                
105
#        system(paste("pksetmask -i ",modvrt," -m ",nObs,
106
#                     " --operator='<' --msknodata 1 --nodata 255 --operator='>' --msknodata 10 --nodata 10 -o ",modtif1,sep=""))
107
        
157

  
108 158
        ## warp to wgs84
109 159
        ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata 255 -dstnodata 255 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",
110 160
            "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5")
climate/research/cloud/MOD09/ee_biastest.R
16 16
## set raster options
17 17
rasterOptions(maxmemory=1e8) 
18 18

  
19
## get CF data
20
##  Get list of available files
21
#df=data.frame(path=list.files("data/mod09cloud",pattern="*.tif$",full=T,recur=T),stringsAsFactors=F)
22
#df[,c("month")]=do.call(rbind,strsplit(basename(df$path),"_|[.]|-"))[,c(6)]
23
#df$date=as.Date(paste(2013,"_",df$month,"_15",sep=""),"%Y_%m_%d")
24
#df=df[order(df$date),]
19
#reg=extent(c(-1434564.00523,1784892.95369, 564861.173869, 1880991.3772))
20
#reg=extent(c(-2187230.72881, 4017838.07688,  -339907.592509, 3589340.63805))  #all sahara
21
#reg=extent(c(-12020769.9608, 10201058.231,  -3682105.25271, 3649806.08382))  #large equitorial
25 22

  
26
#d=stack(df$path,bands=1)
27
#names(d)=df$month
28 23

  
29
reg=extent(c(-1434564.00523,1784892.95369, 564861.173869, 1880991.3772))
30
reg=extent(c(-2187230.72881, 4017838.07688,  -339907.592509, 3589340.63805))  #all sahara
31
reg=extent(c(-12020769.9608, 10201058.231,  -3682105.25271, 3649806.08382))  #large equitorial
32

  
33
#reg=extent(c(-151270.082307, 557321.958761, 1246351.8356, 1733123.75947))
34

  
35
d=stack("data/mcd09/2014113_combined_20002013_MOD09_1-img-0000000000-0000000000.tif")
36

  
37

  
38
names(d)=c("cf","cfsd","nobs","pobs")
39
#cd=coordinates(d)
40
#cd$x=as.integer(cd[,"x"]); cd$y=as.integer(cd[,"y"])
41
#cdr=rasterFromXYZ(cbind(cd,cd))
42
#names(cdr)
43
#d=stack(d,cdr)
44
#obs=crop(raster("data/mcd09/2014113_combined_20002013_MOD09_1-img-0000000000-0000000000.tif",band=4),reg)
45
#levelplot(stack(obs,d))
46

  
47
#plot(obs,d)
48
 
49
#pts=sampleStratified(d[["pobs"]], size=10000, exp=10, na.rm=TRUE, xy=FALSE, ext=NULL, sp=FALSE)
50
pts=sampleRandom(d, size=100000, na.rm=TRUE, xy=T, ext=NULL, sp=T)
51
pts=pts[pts$nobs>0,]  #drop data from missing tiles
52

  
53
#dt=data.frame(cf=as.numeric(values(d[["cf"]])),obs=as.numeric(values(d[["pobs"]])))
54
#dt[,c("x","y")]=coordinates(d)
55
#dt=dt[dt$obs>0,]  #drop data from missing tiles
56

  
57
#s=1:nrow(dt)
58
#s=sample(1:nrow(dt),10000)
59
#s=strata(dt, stratanames=dt$obs, size=1000)
60

  
61
lm1=bam(cf~s(x,y)+pobs+nobs,data=pts@data)
62
summary(lm1)
63
beta1=coef(lm1)["pobs"]; beta1
64
beta2=coef(lm1)["nobs"]; beta2
65

  
66
#d2=d
67
#d2[["pobs"]]=100
68

  
69
#pred1=raster::predict(d,lm1,file=paste("data/bias/pred1.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#lm1=bam(as.vector(d[["cf"]]) ~ as.vector(d[["pobs"]]))
70
#pred2=raster::predict(d2,lm1,file=paste("data/bias/pred2.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#lm1=bam(as.vector(d[["cf"]]) ~ as.vector(d[["pobs"]]))
71
#pred3=overlay(pred1,pred2,function(x,y) x-y,file="data/bias/pred3.tif",overwrite=T)
72
#biasc=overlay(d[["cf"]], pred3,fun=sum,file="data/bias/biasc.tif",overwrite=T)
73

  
74

  
75
getbias=function(cf,cfsd,nobs,pobs) return((100-pobs)*beta1+(4-nobs)*beta2)
76
getbias(cf=c(50,40,30),c(1,1,1),c(3,4,5),c(82,87,87))
77
bias=overlay(d,fun=getbias, unstack=TRUE,file=paste("data/bias/bias.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)
78
dc=overlay(d[["cf"]],bias,fun=function(x,y) x+y,
79
    file=paste("data/bias/biasc.tif",sep=""),
80
    format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9"))
81
#levelplot(stack(d,dc))
82

  
83

  
84
library(multicore)
85

  
86
Mode <- function(x) {
87
      ux <- unique(x)
88
        ux[which.max(tabulate(match(x, ux)))]
89
  }
90
#rasterOptions(maxmemory)#=50GB) 
91 24
## set up processing chunks
92
nrw=nrow(d)
25
nrw=nrow(mod)
93 26
nby=20
94 27
nrwg=seq(1,nrw,by=nby)
95 28
writeLines(paste("Processing ",length(nrwg)," groups and",nrw,"lines"))
96 29

  
97 30
## Parallel loop to conduct moving window analysis and quantify pixels with significant shifts across pp or lulc boundaries
98 31
output=mclapply(nrwg,function(ti){
99
      ## Extract focal areas
100
      ngb=c(51,101)
101
      vals_cf=getValuesFocal(d,ngb=ngb,row=ti,nrows=nby)
32
    ## Extract focal areas
33
    nr=51
34
    nc=101
35
    ngb=c(nr,nc)
36
      vals_cf=getValuesFocal(mod[[c("cf","pobs")]],ngb=ngb,row=ti,nrows=nby)
102 37
      vals_obs=getValuesFocal(obs,ngb=ngb,row=ti,nrows=nby)
103 38
      ## extract bias
104 39
      bias=raster(matrix(do.call(rbind,lapply(1:nrow(vals_cf),function(i) {

Also available in: Unified diff