###load libraries###
library(raster)
library(maptools)
library(rgdal)
library(rgeos)

###define data paths###
rm(list=ls())
shp_data <- "./shp/"
light_data <- "R:/desktop//data/dmsp-ols/old/"
lc_data <- "R:/desktop/data/modis/world/"
C<-grep(".shp$",dir(shp_data),value=TRUE)
lf <- grep(".tif$",paste(light_data,dir(light_data),sep=""),value=TRUE)
lc <- grep(".tif$",paste(lc_data,dir(lc_data),sep=""),value=TRUE)
tt <- as.data.frame(matrix(rep("",length=23),ncol=23))

###define data items###
colnames(tt) <- c("year","light","light63","area_light","area_light63","area","lc0","lc1","lc2","lc3","lc4","lc5","lc6","lc7","lc8","lc9","lc10","lc11","lc12","lc13","lc14","lc15","lc16")

###read shape###
shp<-readOGR("shp",gsub(".shp","",C))
shp<- spTransform(shp, CRS("+proj=longlat +datum=WGS84"))
adm <- shp@data

###create data frame###
d_name <- cbind(as.character(shp$P_name_11),as.character(shp$D_name_11))
d_name <- as.data.frame(matrix(d_name,ncol=2))
colnames(d_name) <- c("P_name_11","D_name_11")
d_name <- d_name[order(d_name$P_name_11,d_name$D_name_11),]
d_name <- unique(d_name)
dat <- cbind(d_name,tt)	

###create data file###		
if(!file.exists("light_lc_district_cambodia_rev.csv")){
	csv<-file("light_lc_district_cambodia_rev.csv",open="wt")
	tmp<-colnames(dat)[1]
	for(i in 2:length(colnames(dat))) tmp<-paste(tmp,colnames(dat[i]),sep=",")
	writeLines(tmp,con=csv)
	close(csv)
}

###aggregate###
for(year in 1992:2012){
	##read global light data##
	light_fn<-grep(year,lf,value=TRUE)
	lc_fn<-grep(year,lc,value=TRUE)
	result_light<-try(light<-raster(light_fn))
	
	##read global land cover data##
	if(year>=2001)
		result_lc<-try(landc<-raster(lc_fn))
	
	##loop by district##
	for(ms in 1:length(dat[,1])){
		#cut and dissolve distirct polygon#
		tmp_poly<-shp[shp$P_name_11==as.character(dat$P_name_11[ms]) & shp$D_name_11==as.character(dat$D_name_11[ms]),]
      	id <- tmp_poly$D_name_11
		tmp_poly <- unionSpatialPolygons(tmp_poly,id)

		#aggregate light# 
		if(year<2001){
			#light
			tmp_light <- extract(light,tmp_poly,small=TRUE)
			s <- NA
			s63 <- NA
			for(ii in 1:length(tmp_light)){
				s <- sum(s,sum(tmp_light[[ii]][tmp_light[[ii]]<255],na.rm=TRUE),na.rm=TRUE)
				s63 <- sum(s63,sum(tmp_light[[ii]][tmp_light[[ii]]==63],na.rm=TRUE),na.rm=TRUE)
			}
			t <- dat[ms,]
			t$year<-year
			t$light<-s
			t$light63<-s63
			t$area <- sum(areaPolygon(tmp_poly)/1000000)	
			t$area_light<-NA
			t$area_light63<-NA
			tmp_light_a <- crop(light,tmp_poly)
			tmp_light_a <- rasterize(tmp_poly,tmp_light_a,mask=TRUE)
			tmp_light_a_p <- rasterToPolygons(tmp_light_a,dissolve=TRUE)	
			try(t$area_light<- sum(areaPolygon(tmp_light_a_p[tmp_light_a_p$layer>0 & tmp_light_a_p$layer<64,])/1000000))
			try(t$area_light63<- sum(areaPolygon(tmp_light_a_p[tmp_light_a_p$layer==63,])/1000000))
			print(t)

			#write to data file#
			write.table(t, file = "light_lc_district_cambodia_rev.csv", append=TRUE,col.names = FALSE, row.names = FALSE, sep = ",")
		
		#aggregate light and land cover#
		}else{
			#light#
			tmp_light <- extract(light,tmp_poly,small=TRUE)
			s <- NA
			s63 <- NA
			for(ii in 1:length(tmp_light)){
				s <- sum(s,sum(tmp_light[[ii]][tmp_light[[ii]]<255],na.rm=TRUE),na.rm=TRUE)
				s63 <- sum(s63,sum(tmp_light[[ii]][tmp_light[[ii]]==63],na.rm=TRUE),na.rm=TRUE)
			}
			t <- dat[ms,]
			t$year<-year
			t$light<-s
			t$light63<-s63
			t$area <- sum(areaPolygon(tmp_poly)/1000000)				
			t$area_light<-NA
			t$area_light63<-NA
			tmp_light_a <- crop(light,tmp_poly)
			tmp_light_a <- rasterize(tmp_poly,tmp_light_a,mask=TRUE)
			tmp_light_a_p <- rasterToPolygons(tmp_light_a,dissolve=TRUE)	
			try(t$area_light<- sum(areaPolygon(tmp_light_a_p[tmp_light_a_p$layer>0 & tmp_light_a_p$layer<64,])/1000000))
			try(t$area_light63<- sum(areaPolygon(tmp_light_a_p[tmp_light_a_p$layer==63,])/1000000))

			#land cover#				
			tmp_lc <- crop(landc,tmp_poly)
			tmp_lc <- rasterize(tmp_poly,tmp_lc,mask=TRUE)
			tmp_lc_p <- rasterToPolygons(tmp_lc,dissolve=TRUE)				
			t$lc16<-t$lc15<-t$lc14<-t$lc13<-t$lc12<-t$lc11<-t$lc10<-t$lc9<-t$lc8<-t$lc7<-t$lc6<-t$lc5<-t$lc4<-t$lc3<-t$lc2<-t$lc1<-t$lc0<-NA		
			try(t$lc0 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==0,]),na.rm=TRUE)/1000000)
			try(t$lc1 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==1,]),na.rm=TRUE)/1000000)
			try(t$lc2 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==2,]),na.rm=TRUE)/1000000)
			try(t$lc3 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==3,]),na.rm=TRUE)/1000000)
			try(t$lc4 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==4,]),na.rm=TRUE)/1000000)
			try(t$lc5 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==5,]),na.rm=TRUE)/1000000)
			try(t$lc6 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==6,]),na.rm=TRUE)/1000000)
			try(t$lc7 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==7,]),na.rm=TRUE)/1000000)
			try(t$lc8 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==8,]),na.rm=TRUE)/1000000)
			try(t$lc9 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==9,]),na.rm=TRUE)/1000000)
			try(t$lc10 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==10,]),na.rm=TRUE)/1000000)
			try(t$lc11 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==11,]),na.rm=TRUE)/1000000)
			try(t$lc12 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==12,]),na.rm=TRUE)/1000000)
			try(t$lc13 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==13,]),na.rm=TRUE)/1000000)
			try(t$lc14 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==14,]),na.rm=TRUE)/1000000)
			try(t$lc15 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==15,]),na.rm=TRUE)/1000000)
			try(t$lc16 <- sum(areaPolygon(tmp_lc_p[tmp_lc_p$layer==16,]),na.rm=TRUE)/1000000)
			print(t)

			#write to data file#
			write.table(t, file = "light_lc_district_cambodia_rev.csv", append=TRUE,col.names = FALSE, row.names = FALSE, sep = ",")
		}
	}
}