DirPath <- "D:/RScripts/migration_map/"
setwd(DirPath)
#=== read in nc data file
inputfile <- c(paste(DirPath,"1904.nc",sep=""),
paste(DirPath,"1949.nc",sep=""),
paste(DirPath,"1956.nc",sep=""),
paste(DirPath,"1994.nc",sep=""),
paste(DirPath,"2000.nc",sep=""),
paste(DirPath,"2005.nc",sep=""),
paste(DirPath,"2010.nc",sep=""),
paste(DirPath,"2015.nc",sep="") )
library(maptools)
## Warning: package 'maptools' was built under R version 3.4.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.4.3
## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()
library(raster)
## Warning: package 'raster' was built under R version 3.4.3
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.4.3
## rgdal: version: 1.2-16, (SVN revision 701)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
## Path to GDAL shared files: C:/Users/ychen/Documents/R/win-library/3.4/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/ychen/Documents/R/win-library/3.4/rgdal/proj
## Linking to sp version: 1.2-7
#get different type pixcels
lu_table <- data.frame(matrix(ncol = 6, nrow = 0))
colnames(lu_table) = c("Year","Forest","Agri","Grass","W","Built")
id <- c(1, 4, 2, 13, 8)
dd <- NA
for (i in 1:8 ) {
LU_P <- as.matrix(raster(inputfile[i], varname="LU_TYPE"))
dd[1] <- as.numeric(substr(inputfile[i], start=nchar(inputfile[i])-6, stop=nchar(inputfile[i])-3 ))
for (j in 1:5 ){
dd[j+1] <- as.numeric(length(LU_P[which(LU_P==id[j])]))
}
dd<- as.numeric(dd)
dd_all <- sum(as.numeric(dd[2:6]) )
dd[2:6] <- dd[2:6]/dd_all
new_row <- data.frame(Year=dd[1],Forest=dd[2],
Agri=dd[3],Grass=dd[4],
Water=dd[5], Built=dd[6])
lu_table <- rbind(lu_table,new_row)
}
## Loading required namespace: ncdf4
layout(matrix(data=c(1),nrow=1, ncol=1, byrow=TRUE),
widths=c(1,1,1), heights=c(1,1,1))
par( oma=c(0,0,0,0), mgp=c(2.5,1,0), mar=c(1.0,1.0,0.0,0.0),plt = c(0.12,0.9,0.12,0.9))
par(cex.lab=2.0, cex.axis=1.5,cex=1.1)
par(xpd=TRUE)
plot(Forest~Year,data=lu_table,type="l",col="#006400",
ylim=c(0.3,1),ylab="Cover Fration",
xlim=c(1900,2020))
#plot()
x<-NA
y<-NA
for (i in 1:7) {
#plot forest
x[1]<-lu_table$Year[i]
x[2]<-lu_table$Year[i]
x[3]<-lu_table$Year[i+1]
x[4]<-lu_table$Year[i+1]
y[1]<-0.3
y[2]<-lu_table$Forest[i]
y[3]<-lu_table$Forest[i+1]
y[4]<-0.3
polygon(x = x, y = y,
col = "#006400", border = "gray",lty="dashed")
#plot Agri
x[1]<-lu_table$Year[i]
x[2]<-lu_table$Year[i]
x[3]<-lu_table$Year[i+1]
x[4]<-lu_table$Year[i+1]
y[1]<-lu_table$Forest[i]
y[2]<-lu_table$Forest[i]+lu_table$Agri[i]
y[3]<-lu_table$Forest[i+1]+lu_table$Agri[i+1]
y[4]<-lu_table$Forest[i+1]
polygon(x = x, y = y,
col = "#FFB90F", border = "gray",lty="dashed")
#plot Grass
x[1]<-lu_table$Year[i]
x[2]<-lu_table$Year[i]
x[3]<-lu_table$Year[i+1]
x[4]<-lu_table$Year[i+1]
y[1]<-lu_table$Forest[i]+lu_table$Agri[i]
y[2]<-lu_table$Forest[i]+lu_table$Agri[i]+lu_table$Grass[i]
y[3]<-lu_table$Forest[i+1]+lu_table$Agri[i+1]+lu_table$Grass[i+1]
y[4]<-lu_table$Forest[i+1]+lu_table$Agri[i+1]
polygon(x = x, y = y,
col = "#A2CD5A", border = "gray",lty="dashed")
#plot Water
x[1]<-lu_table$Year[i]
x[2]<-lu_table$Year[i]
x[3]<-lu_table$Year[i+1]
x[4]<-lu_table$Year[i+1]
y[1]<-lu_table$Forest[i]+lu_table$Agri[i]+lu_table$Grass[i]
y[2]<-lu_table$Forest[i]+lu_table$Agri[i]+lu_table$Grass[i]+lu_table$Water[i]
y[3]<-lu_table$Forest[i+1]+lu_table$Agri[i+1]+lu_table$Grass[i+1]+lu_table$Water[i+1]
y[4]<-lu_table$Forest[i+1]+lu_table$Agri[i+1]+lu_table$Grass[i+1]
polygon(x = x, y = y,
col = "#0000FF", border = "gray",lty="dashed")
#plot Soil
x[1]<-lu_table$Year[i]
x[2]<-lu_table$Year[i]
x[3]<-lu_table$Year[i+1]
x[4]<-lu_table$Year[i+1]
y[1]<-lu_table$Forest[i]+lu_table$Agri[i]+lu_table$Grass[i]++lu_table$Water[i]
y[2]<-lu_table$Forest[i]+lu_table$Agri[i]+lu_table$Grass[i]+lu_table$Water[i]+lu_table$Built[i]
y[3]<-lu_table$Forest[i+1]+lu_table$Agri[i+1]+lu_table$Grass[i+1]+lu_table$Water[i+1]+lu_table$Built[i+1]
y[4]<-lu_table$Forest[i+1]+lu_table$Agri[i+1]+lu_table$Grass[i+1]+lu_table$Water[i+1]
polygon(x = x, y = y,
col = "#FF1493", border = "gray",lty="dashed")
}
#add blacklines
for (i in 1:8) {
x[1]<-lu_table$Year[i]
y[1]<-lu_table$Forest[i]
y[2]<-lu_table$Forest[i]+lu_table$Agri[i]
y[3]<-lu_table$Forest[i]+lu_table$Agri[i]+lu_table$Grass[i]
y[4]<-lu_table$Forest[i]+lu_table$Agri[i]+lu_table$Grass[i]+lu_table$Water[i]
off_x <- -0.0
off_y <- +0.01
cexset <- 1.0
#add text for forest
pp_text<-paste(format(x=lu_table$Forest[i]*100., digits=3,width=4,format="f"), "",sep="")
text(x=x[1]+off_x,y=y[1]+off_y,pp_text,cex=cexset)
#add text for agri
pp_text<-paste(format(x=lu_table$Agri[i]*100., digits=3,width=4,format="f"), "",sep="")
text(x=x[1]+off_x,y=y[2]+off_y,pp_text,cex=cexset)
#add text for grass
pp_text<-paste(format(x=lu_table$Grass[i]*100., digits=3,width=4,format="f"), "",sep="")
text(x=x[1]+off_x,y=y[3]+off_y,pp_text,cex=cexset)
#add text for water
pp_text<-paste(format(x=lu_table$Water[i]*100., digits=3,width=4,format="f"), "",sep="")
text(x=x[1]+off_x,y=y[4]+off_y,pp_text,cex=cexset)
#add text for Built
pp_text<-paste(format(x=lu_table$Built[i]*100., digits=3,width=4,format="f"), "",sep="")
text(x=x[1]+off_x,y=1.0+off_y,pp_text,cex=cexset)
}
lines(Forest~Year,data=lu_table,type="l",col="black")
lines(Forest+Agri~Year,data=lu_table,type="l",col="black")
lines(Forest+Agri+Grass~Year,data=lu_table,type="l",col="black")
lines(Forest+Agri+Grass+Water~Year,data=lu_table,type="l",col="black")
lines(Forest+Agri+Grass+Water+Built~Year,data=lu_table,type="l",col="black")
# so turn off clipping:
par(xpd=TRUE)
MyCol=c("#006400","#FFB90F","#A2CD5A","#0000FF","#FF1493")
MyLab=c(c("Forest","Agri","Grass","Water","Built/Soil"))
op <- par(cex = 1.2)
legend(x=1900,y=1.11,
MyLab,fill=MyCol,
horiz=T,
bty = "n")
dev.copy2pdf(file = "Net_change_plot.pdf", height=10, width=15)
## png
## 2