R/mpaprocess.R

#' This function compute a set of indicators based on a series of map of a EMIS/GMIS parameter
#' and draw some maps and plots related to the parameter. 
#'
#' The parameters are monthly data available on EMIS or GMIS.
#' "name" gives the short name of the parameter processed.
#' "timerange" is the time range of the parameter.
#' "pixelnb" gives the number of 4 by 4 km pixels covering the study area: this area is the geographical extent of the MPA polygon (ie a rectangle encompassing the MPA).
#' "fracnbna" is the percentage of missing data in the MPA polygon over the time range. 
#' "fractsna" is the percentage of missing value for the time series of the parameters. This time serie is calculated using the monthly mean of the pixel values INSIDE the MPA polygon.
#' "mean" is the average value of the parameter INSIDE the MPA polygon over the time range.
#' "sd" is the standard deviation of the parameter values INSIDE the MPA polygon over the time range.
#' "min" is the minimum of the parameter values INSIDE the MPA polygon over the time range.
#' "max" is the maximum of the parameter values INSIDE the MPA polygon over the time range.
#' The monthly time series X(t) of the parameter is decomposed as  X(t)=S(t)+T(t)+I(t) (see Vantrepotte and Melin, 2009), where S, T, and I represent, respectively, a seasonal signal, a trend, and an irregular (or residual) component (or remainder). These components are presented in the figure seasonal composition. In this framework:
#' "varseason" is the percentage of variance associated the seasonal component and "vartrend" the percentage of variance associated with the linear component (the trend).
#' A statistical test (Kendall's tau test) is performed on the linear component in order to assess the significance of a monotic trend:
#' "trendtest" gives the p.value associated to this test.
#' "senslope" gives the value of the Sen slope (alternately, Theil slope). This value is the median
#' slope joining all pairs of observations and is expressed by
#' quantity per year.
#
#' References:
#' Vantrepotte, V. & Melin, F. Temporal variability of 10-year global SeaWiFS time-series of phytoplankton chlorophyll a concentration ICES J. Mar. Sci., 2009, 66, 1547-1556
#      
#' @param name A character vector of the shortname of the variable
#' @param resolution A character vector giving the spatial resolution of the data: "4km" or "2km"
#' @param startdate A character vector of the month who begins the series("YYYY-MM")
#' @param enddate A character vector of the month who ends the series("YYYY-MM")
#' @param wdpaid An integer of the ID of the marine protected area corresponding to the wdpaid of the mpa 
#' 
#' @export
#' @return A list of object containing the plot and the statistics
#' @keywords EMIS
#' @keywords wcs-t
#' @examples
#'   \dontrun{
#'      #analysis of the MODIS sea surface temperature at 2 km 2005 and 2006 on the Pantelleria mpa (Italy)
#'	pltstat<-mpaprocess("EMIS_T_SST","2km","2005-01","2007-12",555540558)
#'	#map of the whole series
#'	pltstat[[1]]$pltall
#'	#map of the average SST
#'	pltstat[[1]]$pltmean
#'	#map of the climatology
#'	pltstat[[1]]$pltclim
#'	#some statistics
#'	pltstat[[2]]$stat
#'	#time series decomposition of the parameters averaged on the MPA area
#'	pltstat[[2]]$pltts
#'	#extraction of the MODIS sea surface temperature at 4 km from GMIS for 2005-2007 in the the Bird Island area marine protected area (South Africa)
#'	pltstat<-mpaprocess("GMIS_T_SST","4km","2005-01","2007-12",306180)
#'	#map of the whole series
#'	pltstat[[1]]$pltall
#'	#map of the average SST
#'	pltstat[[1]]$pltmean
#'	#map of the climatology
#'	pltstat[[1]]$pltclim
#'	#some statistics
#'	pltstat[[2]]$stat
#'	#time series decomposition of the parameters averaged on the MPA area
#'	pltstat[[2]]$pltts
#'   } 
#'
mpaprocess<-function(name="EMIS_A_CHLA",resolution="4km",startdate="2005-09",enddate="2005-10",wdpaid=555540558){ 
 #check the mpa
 checkmpa1<- length(which(mpa1@data$wdpaid==wdpaid))
 checkmpa2<- length(which(mpa2@data$wdpaid==wdpaid))
 if((checkmpa1+checkmpa2)!=1){
  	print(paste("The MPA", wdpaid, "do not exist in the mpa database. Please check your wdpaid"))
 }else{
 	#extract the polygon of the mpa corresponding to the wdpaid
 	if(wdpaid<=4000){
  		mpa<-mpa1[mpa1@data$wdpaid==wdpaid,]
  		proj4string(mpa)<-CRS("+proj=longlat +datum=WGS84")
 	}else{
  		mpa<-mpa2[mpa2@data$wdpaid==wdpaid,]
  		proj4string(mpa)<-CRS("+proj=longlat +datum=WGS84")
 	}
 	#download the parameter from EMIS/GMIS
 	imgs<-mpaextract(name=name,resolution=resolution,startdate=startdate,enddate=enddate,wdpaid=wdpaid) 
 	#grap the unit
	if(substr(name,1,4)=="EMIS"){
 		if(resolution=="4km"){data_emis<-data_emis_4km}
        	if(resolution=="2km"){data_emis<-data_emis_2km}
 		idvar<-grep(name,data_emis$shortname)
		if(length(grep("log",data_emis$unit[idvar],ignore.case=TRUE))>0){
                	unite<-substr(data_emis$unit[idvar],7,nchar(data_emis$unit[idvar]))
        	}else{
			unite<-data_emis$unit[idvar]
		}	
 	}
	if(substr(name,1,4)=="GMIS"){
 		if(resolution=="4km"){data_gmis<-data_gmis_4km}
        	if(resolution=="9km"){data_gmis<-data_gmis_9km}
 		idvar<-grep(name,data_gmis$shortname)
		if(length(grep("log",data_gmis$unit[idvar],ignore.case=TRUE))>0){
                	unite<-substr(data_gmis$unit[idvar],7,nchar(data_gmis$unit[idvar]))
        	}else{
			unite<-data_gmis$unit[idvar]
		}	
 	}
  	#plot the data
	plt<-mpaprocessplot(imgs=imgs,mpa=mpa,name=name,unite=unite,logscale=FALSE)
 	stat<-mpaprocessstat(imgs=imgs,mpa=mpa,name=name,unite=unite)
	return(list(plt,stat))
 }
}

#' This function draw figures based on a series of map of a EMIS/GMIS parameter for a marine protected area
#'
#' @param imgs A rasterstack of the series of map extracted using mpaextract 
#' @param mpa A SpatialPolygonsDataFrame of the marine protected area of interest
#' @param name A character vector of the variable name
#' @param unite A character vector of the unit of variable name
#' @param logscale A boolean defining the use of a logscale on a graph (TRUE) or not (FALSE)
#' 
#' @export
#' @return A list of plots
#' @keywords EMIS
#' @keywords GMIS
#' @keywords wcs-t
#' @examples
#'      #analysis of the MODIS sea surface temperature at 2 km 2009 and 2012 on the Pantelleira mpa (Italy)
#'	plt<-mpaprocessplot(imgs=pantelleria_sst,mpa=pantelleria_mpa,name="EMIS_T_SST",unite="oC",logscale=FALSE)
#'	#map of the whole series
#'	plt[[1]]
#'	#map of the average SST
#'	plt[[2]]
#'	#map of the climatology
#'	plt[[3]]
#'	#boxplot of the climatology
#'	plt[[4]]
#'
mpaprocessplot<-function(imgs,mpa,name,unite,logscale){
 if(dim(imgs)[3]<12){
  	print("Series is too short (ie less than 12 months)")
 }else{
 	#precomputation
	 #average
	 imgsmean<-raster::mean(imgs,na.rm=T)
	 # monthly climatology
	 imgsclim<-stackApply(imgs,indices=as.numeric(substr(names(imgs),7,8)),fun=mean,na.rm=TRUE)
	 names(imgsclim)<-paste("Mean",unique(as.numeric(substr(names(imgs),7,8))))
 	#plot

	 titre<-paste(name," (",unite,")",sep="")
	 pltall<-rasterVis::levelplot(imgs,zscaleLog=logscale,contour=T,col.regions=topo.colors(100),names=substr(names(imgs),2,8),main=titre)
	 titre<-paste(name," (",unite,") ",substr(names(imgs)[1],2,8),"-",substr(names(imgs)[dim(imgs)[3]],2,8)," average",sep="")
	 pltmean<-rasterVis::levelplot(imgsmean,margin=F,zscaleLog=logscale,contour=T,col.regions=topo.colors(100),main=titre)
	 titre<-paste(name," (",unite,") monthly climatology",sep="")
	 pltclim<-rasterVis::levelplot(imgsclim,zscaleLog=logscale,contour=T,col.regions=topo.colors(100),main=titre,names=names(imgsclim))
	 titre<-paste(name," (",unite,") monthly boxplot",sep="")
	 pltbw<-rasterVis::bwplot(crop(imgsclim,extent(mpa)),scales=list(x=list(labels=names(imgsclim))),main=titre)
 	#build terretrial polygon
	 #mappoly<-map("worldHires",fill=T,plot=FALSE,xlim=c(extent(mpa)@xmin-1,extent(mpa)@xmax+1),ylim=c(extent(mpa)@ymin-1,extent(mpa)@ymax+1))
	 #IDs <- sapply(strsplit(mappoly$names, ":"), function(x) x[1])
	 #coast<- map2SpatialPolygons(mappoly, IDs=IDs, proj4string=CRS("+proj=longlat +datum=WGS84"))
 	 #add the coast and the mpa to the maps
 	 #pltall<-pltall+latticeExtra::layer(sp.polygons(coast,fill="grey",col="grey"))+latticeExtra::layer(sp.polygons(mpa))
 	 #pltmean<-pltmean+latticeExtra::layer(sp::sp.polygons(coast,fill="grey",col="grey"))+latticeExtra::layer(sp::sp.polygons(mpa))
 	 #pltclim<-pltclim+layer(sp.polygons(coast,fill="grey",col="grey"))+layer(sp.polygons(mpa))
         #return the plot in a list
         return(list(pltall=pltall,pltmean=pltmean,pltclim=pltclim,pltbw=pltbw))
 }
}

#' This function compute some indicators for a series of map of a EMIS/GMIS parameter for a marine protected area
#'
#' The parameters are monthly data available on EMIS or GMIS.
#' "name" gives the short name of the parameter processed.
#' "timerange" is the time range of the parameter.
#' "pixelnb" gives the number of 4 by 4 km pixels covering the study area: this area is the geographical extent of the MPA polygon (ie a rectangle encompassing the MPA).
#' "fracnbna" is the percentage of missing data in the MPA polygon over the time range. 
#' "fractsna" is the percentage of missing value for the time series of the parameters. This time serie is calculated using the monthly mean of the pixel values INSIDE the MPA polygon.
#' "mean" is the average value of the parameter INSIDE the MPA polygon over the time range.
#' "sd" is the standard deviation of the parameter values INSIDE the MPA polygon over the time range.
#' "min" is the minimum of the parameter values INSIDE the MPA polygon over the time range.
#' "max" is the maximum of the parameter values INSIDE the MPA polygon over the time range.
#' The monthly time series X(t) of the parameter is decomposed as  X(t)=S(t)+T(t)+I(t) (see Vantrepotte and Melin, 2009), where S, T, and I represent, respectively, a seasonal signal, a trend, and an irregular (or residual) component (or remainder). These components are presented in the figure seasonal composition. In this framework:
#' "varseason" is the percentage of variance associated the seasonal component and "vartrend" the percentage of variance associated with the linear component (the trend).
#' A statistical test (Kendall's tau test) is performed on the linear component in order to assess the significance of a monotic trend:
#' "trendtest" gives the p.value associated to this test.
#' "senslope" gives the value of the Sen slope (alternately, Theil slope). This value is the median
#' slope joining all pairs of observations and is expressed by
#' quantity per year.
#
#' References:
#' Vantrepotte, V. & Melin, F. Temporal variability of 10-year global SeaWiFS time-series of phytoplankton chlorophyll a concentration ICES J. Mar. Sci., 2009, 66, 1547-1556

#' @param imgs A rasterstack of the series of map extracted using mpaextract 
#' @param mpa A SpatialPolygonsDataFrame of the marine protected area of interest
#' @param name A character vector of the variable name
#' @param unite A character vector of the unit of variable name
#' 
#' @export
#' @return A list of 2 objects: a dataframe and a plot
#' @keywords EMIS
#' @keywords GMIS
#' @keywords wcs-t
#' @examples
#'      #analysis of the MODIS sea surface temperature at 2 km 2009 and 2012 on the Pantelleira mpa (Italy)
#'	stat<-mpaprocessstat(imgs=pantelleria_sst,mpa=pantelleria_mpa,name="EMIS_T_SST",unite="oC")
#'	#statistics
#'	stat[[1]]
#'	#plot of the time series decomposition of the parameter averaged on the MPA surface
#'	stat[[2]]
#'
mpaprocessstat<-function(imgs,mpa,name,unite){
  #############################################
  # compute indicators"
  # check data availability in the polygon and create a dataframe with the corresponding pixel"
  datinpoly<-raster::extract(imgs,mpa,df=TRUE,cellnumbers=TRUE)
  tabdat<-as.vector(as.matrix(datinpoly[,3:ncol(datinpoly)]))
  pixelnb<-nrow(datinpoly)
  #compute the time series for the given pixels
  tabts<-apply(datinpoly[,3:ncol(datinpoly)],2,mean,na.rm=T)
  #number of missing value
  fracnbna<-sum(is.na(tabts))/length(tabts)
  #if more than 20% of the value are not present do not time series analysis
  if(fracnbna>.2){
  print("not enough data for time series analysis")
   vseason<-NA	
   vtrend<-NA	
   vtot<-NA	
   testtrend<-data.frame(p.value=NA,sen.slope=NA)	
   pltts<-NA
  }else{
   startyear<-as.numeric(substr(names(imgs)[1],2,5))
   startmonth<-as.numeric(substr(names(imgs)[1],7,8))
   tsm<-na.approx(ts(tabts, frequency = 12, start = c(startyear, startmonth)),na.rm=FALSE)
   tsmdec<-stl(tsm, s.window=12,na.action=na.omit)
   tsseason<-tsmdec$time.series[,1]
   tsremind<-tsm-tsseason
   linfit<-lm(tsremind~c(1:length(tsremind)))
   tstrend<-ts(predict(linfit), frequency = 12, start = c(startyear,startmonth))
   tsremind<-tsm-tsseason-tstrend
   #check this
   pipo1<-data.frame(time=as.numeric(time(tsm)),value=as.numeric(tsm),name="Time series",type="Original time series and trend")
   pipo2<-data.frame(time=as.numeric(time(tsseason)),value=as.numeric(tsseason),name="Seasonal component",type="Seasonal and remainder components")
   pipo3<-data.frame(time=as.numeric(time(tstrend)),value=as.numeric(tstrend),name="Trend",type="Original time series and trend")
   pipo4<-data.frame(time=as.numeric(time(tsremind)),value=as.numeric(tsremind),name="Residual",type="Seasonal and residual components")
   pipoall<-rbind(pipo1,pipo2,pipo3,pipo4)
   pltts<-ggplot(pipoall,aes(x=time,y=value,group=name))+geom_line()+facet_grid(name~., scales="free")+ylab(paste(name," (",unite,")",sep=""))
   vseason<- var(tsseason,na.rm=T)
   vtrend<- var(tstrend,na.rm=T)
   vremind<- var(tsremind,na.rm=T)
   vtot<-vseason+vtrend+vremind
   #some stat
   #test linear trend
   testtrend<-mannKen(tstrend)
   #return results
  }
  datstat<-data.frame(name=name,unite=unite,timerange=paste(substr(names(imgs)[1],2,8),substr(names(imgs)[dim(imgs)[3]],2,8),sep="/"),pixelnb=pixelnb,fracnbna=round(sum(is.na(tabdat))/length(tabdat)*100,2),fractsna=fracnbna*100,mean=mean(tabdat,na.rm=T),sd=sd(tabdat,na.rm=T),min=min(tabdat,na.rm=T),max=max(tabdat,na.rm=T),varseason=100*vseason/vtot,vartrend=100*vtrend/vtot,trendtest=testtrend$p.value,sen.slope=testtrend$sen.slope)
  return(list(stat=datstat,pltts=pltts))
}
ldbk/EMISR documentation built on May 20, 2019, 11:28 p.m.