SIMLANDER R script

APoLUS

To download APoLUS, a multiple land use version of SIMLANDER, go here

SIMLANDER

Download latest user manual in pdf format (March 2017) here.

DOWNLOAD SIMLANDER LATEST VERSION

Download Simlander latest version v.1.0.5 (Rscript and sample data)

DOWNLOAD ALL VERSIONS

NB: These previous versions are downloadable from WordPress. WordPress does not permit .zip extensions, so after downloading (if you are on a windows system make sure your file extensions are not hidden – if they are change this in Folder Options) you will need to change the file extension (which is .doc) back to .zip. Then you can unzip the file which contains the data and scripts.

Download Simlander latest version v.1.0.4 (Rscript and sample data)

Simlander version 1.0.3 (Rscript and sample data)

Simlander version 1.0.2 (Rscript and sample data)

Simlander version 1.0.1 (Rscript and sample data)

—————————————————————————————-

SIMLANDER R script version 1.0.5

In Version 1.05, the main changes are in the way the demand is calculated, and a lot of general tidying up. The script should be alot easier to follow and contains much less redundant code

Also added the following code to remove the over-allocation issue

#test for duplicates that inflate the number of cells allocated
difftrans 0) {
result2 <- head(result,-difftrans) #remove the duplicates from the end of the file (the weakest candidate cells)
result <- result2
}

Richard Hewitt 28-March-2017

#begin script
#set workspace and libraries
#setwd("/media/hd_sdb1/richard/SIMLANDeR/R_CA_model/R-working-dir")
require(raster)
#begin import data
lu1 <- raster("lu56.asc")
lu2 <- raster("lu99.asc")
lu3 <- raster("lu03.asc")
lu4 <- raster("lu07.asc")
distriv <- raster("distriv.asc")
distroad <- raster("distroad.asc")
elev <- raster("elevation.asc")
rainfall <- raster("precip_71_00.asc")
slope <- raster("slope.asc")
temper <- raster("tempmedan.asc")
#end import data
#———————–
#begin tidy up a bit – calculate a change map, set null values for features and zoned areas
change56_99 <- lu2-lu1
fun <- function(x) { x[x<0] <- NA; return(x) }
ch56_99<-calc(change56_99, fun)
distroad_null <- mask(distroad,ch56_99)
distriv_null <- mask(distriv,ch56_99)
elev_null <- mask(elev,ch56_99)
rain_null <- mask(rainfall,ch56_99)
slope_null <- mask(slope,ch56_99)
temp_null <- mask(temper,ch56_99)
#end tidy up
#———————–
#begin calculate accessibility
roadsco <- reclassify(distroad_null,c(-Inf,100,1,100,300,0.9,300,800,0.8,800,1500,0.7,1500,2500,0.6,2500,5000,0.5,5000,8000,0.4,8000,10000,0.3,10000,12000,0.2,12000,Inf,0.1)) #setting the coefficients
accessibility <- 1+((distroad_null/roadsco)^-1)
funacc 2] <- 2; return(x) } #to make sure we do not have infite values at distance 0
acc<-calc(accessibility, funacc)
#end calculate accessibility
#———————–
#begin calculate suitability
slope_suit <- reclassify(slope_null,c(-Inf,0,1,0,5,0.9,5,10,0.8,10,15,0.5,15,20,0.3,20,Inf,0.1))
model_suitability <- slope_suit
#end calculate suitability
#———————–
#set weights matrix for neighbourhood
w <- matrix(c(0,0,50,0,0,0,50,50,50,0,50,50,500,50,50,0,50,50,50,0,0,0,50,0,0), nr=5,nc=5)
#———————–
#preparations for simulation
T1 <- lu1 #set T1 to first map lu1. Change this if we want to start from a different map.
print (paste("simulation starting at ", Sys.time(), sep=""))
startyear <- 1956
endyear <- 1999
total <- (endyear)-(startyear)
#——————————-calculate demand from lu1 and lu2 (must have 2 rows and 2 columns with urban land in row 2, column 2)
dflu1 <- as.data.frame(freq(lu1))
dflu2 <- as.data.frame(freq(lu2))
urbdemand <- as.numeric(dflu1[2,2])
finaldemand <- as.numeric(dflu2[2,2])
andem <- (finaldemand – urbdemand)/total
andem <- round(andem,0)
#——————————-end calculate demand from lu1 and lu2
print("set demand, starting timer..")
ptm <- proc.time()
#end preparations
#———————–
pb <- txtProgressBar(min = 0, max = total, style = 3)
for (i in 1:total){ #begin running simulation
urbdemand1 <- urbdemand + andem*i
#begin neighbourhood block
n <- focal(T1, w=w)
nhood <- cover(n, T1)
model_nhood <- nhood
#end neighbourhood block
#———————–
#begin random block
x <-runif(ncell(T1)) #corrected to calculate directly from number of cells in the map
weibull <- 1+(-log(1-(x)))*exp(1/2)
funselect <- function(x) { x[x!=0] <- NA; return(x) } #extract only vacant areas (value 0) so that existing functional land uses are not randomized. Set everything else to NA
vacants <- calc(T1, funselect) #BE CAREFUL! NEED TO TAKE VACANTS FROM T1, NOT LU56
random <- lu1 #just copying lu1 as a skeleton map
values(random)<-weibull
model_random <- mask(random, vacants) #generating a mask to apply NAs to the random layer.
model_random <- cover(model_random,T1) #cover to fill the NAs back in with the original values from T1.
#end random block
#———————–
#begin transition potential calculation
model_TTP <- (acc)*(model_suitability+1)*(model_nhood+1)*(model_random)
fun <- function(x) { x[is.na(x)] <- 0; return(x)}
TTP <- calc(model_TTP, fun)
#Select the highest 164 values from the TTP map
#r1demand <- 164 #this is the number of cells we want to allocate from r1
x <- as.matrix(TTP)
n <- urbdemand1 #the demand calculated as above
x2 <- sort(-x, partial = n)
x2h <- -sort(x2[1:n])
ix <- which(x %in% x2h)
rowsT1 <- nrow(T1)
colsT1 <- ncol(T1)
ro <- ix %% rowsT1
ro <- ifelse(ro == 0L, rowsT1, ro)
co <- 1 + ix %/% rowsT1
x3 <- x[ix]
d <- data.frame(row = ro, col = co, x = x3)
result <- d[rev(order(d$x)), ]
#———————–
#test for duplicates that inflate the number of cells allocated
difftrans 0) {
result2 <- head(result,-difftrans) #remove the duplicates from the end of the file (the weakest candidate cells)
result <- result2
}
#turn selected n values in the dataframe into a matrix and then to a raster
x.mat <- matrix(0, rowsT1, colsT1) #create a matrix with the right number of rows and columns and fill with 0 values
#x.mat[cbind(result$row,result$col)] <-result$x #works just fine.
x.mat[cbind(result$row,result$col)] <-1 #but actually we want values to be 1 (urban land), not the TP value.
#which(!is.na(x.mat), arr.ind =TRUE) #pick out the non-NA values. Useful for testing that this has actually worked.
r <- raster(x.mat)
extent(r) <- extent(T1)
newdata <- r
newdata <- mask(newdata, lu1) #generating a mask to apply NAs to the final map layer.
T1 <- newdata #directly allocated all the cells at their most favourable locations according to TTP

filen <- paste("lu", (startyear + i), ".png", sep="")
#filenasc <- paste("lu", (startyear + i), ".asc", sep="")
#LUout <- writeRaster(T1, filename=(filenasc), format="ascii", overwrite=TRUE)
#filen <- paste(filen,".png", sep="")
png(filename = paste(filen),width = 1200, height = 1200, bg="white") #Plot each land use map to be able to make an animation.
#plot(T1)
#plot(T1,breaks=breakpoints,col=colors,main=(startyear+ii))
plot(T1,main=(startyear+i))
dev.off()

removeTmpFiles(h=5)
Sys.sleep(0.1)
setTxtProgressBar(pb, i)
print (paste("FINISHED: landuse simulation for ", (startyear + i), sep=""))
print(proc.time() – ptm)
}
close(pb)
png(filename = paste("lu",startyear,".png",sep=""),width = 1200, height = 1200, bg="white") #Plot initial map also to animations file
plot(lu1, main = paste("initial land use map – ",startyear,sep="")) #lu1 for calibration, lu2 for simulation phase
dev.off()
#eval_anim <- paste("system ('convert -delay 200 -loop 0 *.png lu",startyear,"_",endyear,".gif')",sep="") #creates an animated gif using imagemagick, comment back in if you are on a linux system and have imagemagick installed
#eval(parse(text=eval_anim)) #comment back in if you are on a linux system and have imagemagick installed
#simlu2011 <- T1
sim99 <- T1
sss <- stack(lu2,sim99)
#name1 <- paste("initial land use ",startyear,sep="")
name1 <- paste("real land use ",endyear,sep="")
name2 <- paste("simulated land use ",endyear,sep="")
names(sss) <- c(name1,name2)
print(plot(sss)) #without the print command, the output will not appear on the screen
print (paste("simulation complete at ", Sys.time(), sep=""))
#print(plot(sss,breaks=breakpoints,col=colors)) #same as above, but just in case you have colors and breakpoints defined
#END OF SCRIPT–

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s