paper.R
In osc: Orthodromic Spatial Clustering

### R code from vignette source 'paper.rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: paper.rnw:39-44
###################################################
library(osc)
data(exampledata)
str(exampledata)
pop.list <- cca(exampledata[,1:2],s=1)
str(pop.list)


###################################################
### code chunk number 2: paper.rnw:49-56
###################################################
  palette(rainbow(12))
  pdf(file="pics/exdat1.pdf", paper="special", width=8, height=4)
  par(mfrow=c(1,2))
  plot(exampledata$x,exampledata$y,col="darkblue", pch=15, xlab="", ylab="", cex=1.2)
  plot(pop.list$cluster$long,pop.list$cluster$lat,col=pop.list$cluster$cluster_id, pch=15, xlab="", ylab="", cex=1.2)
    dev.off()
  palette("default")


###################################################
### code chunk number 3: paper.rnw:69-76
###################################################
#initiate empty matrix
exampledata.pop <- matrix(0, nrow=max(exampledata$x), ncol=max(exampledata$y))

#restructure data
for(i in 1:NROW(exampledata)){
  exampledata.pop[exampledata$x[i],exampledata$y[i]] <- exampledata$z[i]
}


###################################################
### code chunk number 4: paper.rnw:81-83
###################################################
example.result <- cca(exampledata.pop, s=1)
str(example.result)


###################################################
### code chunk number 5: paper.rnw:88-93
###################################################
  pdf(file="pics/exdat2.pdf", paper="special", width=8, height=4)
  par(mfrow=c(1,2))
  image(x=1:20, y=1:20,exampledata.pop, xlab="", ylab="")
  image(x=1:20, y=1:20,example.result, col=c("white",rep(rainbow(12),2)), xlab="", ylab="")
  dev.off()


###################################################
### code chunk number 6: paper.rnw:112-118
###################################################
# create a raster and set the projection information
raster <- raster(extent(0,5,0,5),nrow=5,ncol=5)
raster[c(1,2,3,5,6,10,17,18,22,23,24)] <- 1
proj4string(raster) <- CRS("+proj=longlat")
# get a feeling for the dimensions
summary(distance(raster)[])


###################################################
### code chunk number 7: paper.rnw:121-126
###################################################
# cluster all cells with value 1
# for various cluster distances of 150, 240 and 111 km
cluster <- cca(raster, cell.class=1, s=1.5e+05, unit="m")
cluster2 <- cca(raster, cell.class=1, s=2.4e+05, unit="m")
cluster3 <- cca(raster, cell.class=1, s=1.11e+05, unit="m")


###################################################
### code chunk number 8: paper.rnw:133-134
###################################################
pixel <- cca(raster, cell.class=1, s=1)


###################################################
### code chunk number 9: paper.rnw:136-137
###################################################
str(pixel)


###################################################
### code chunk number 10: paper.rnw:142-151
###################################################
pdf("pics/raster.pdf", width=9, height=3)
par(mfrow=c(1,3))
raster[cellFromXY(raster,cluster$cluster[,1:2])] <- cluster$cluster[,3]
plot(raster,col=rainbow(11),legend=FALSE, main="s = 150 km")
raster[cellFromXY(raster,cluster2$cluster[,1:2])] <- cluster2$cluster[,3]
plot(raster,col=rainbow(11),legend=FALSE, main="s = 240 km")
raster[cellFromXY(raster,cluster3$cluster[,1:2])] <- cluster3$cluster[,3]
plot(raster,col=rainbow(11),legend=FALSE, main="s = 111 km")
dev.off()


###################################################
### code chunk number 11: paper.rnw:164-166
###################################################
data("landcover")
cities <- cca(landcover, cell.class=1, s=2000, unit="m")


###################################################
### code chunk number 12: paper.rnw:171-174
###################################################
str(cities)
result <- landcover*NA
result[cellFromXY(result, cities$cluster[,1:2])] <- cities$cluster[,3]


###################################################
### code chunk number 13: paper.rnw:179-185
###################################################
pdf("pics/landcover.pdf", width=8, height=6)
par(mfrow=c(1,2))
cols <- c("lightblue", "darkred","yellow","orange","green","darkgreen")
plot(landcover, col=cols)
plot(result, col=rainbow(9))
dev.off()


###################################################
### code chunk number 14: paper.rnw:197-200
###################################################
library(maps)
head(world.cities)
str(world.cities)


###################################################
### code chunk number 15: paper.rnw:203-205
###################################################
city.cluster <- cca(data = world.cities[,c(5,4,1:3,6)], s = 0.5)



###################################################
### code chunk number 16: paper.rnw:209-212
###################################################
city.cluster.m <- cca(data = world.cities[,c(5,4,1:3,6)], s = 56*1e+03,
                      res.x = 0.05, res.y = 0.05, unit = "m")



###################################################
### code chunk number 17: paper.rnw:215-228
###################################################
pdf("pics/point_data.pdf", width=8, height=4)
par(mfrow=c(1,2))
plot(x = city.cluster[[1]]$long,
     y = city.cluster[[1]]$lat,
     col=city.cluster[[1]]$cluster_id,
     main = "d = 0.5°",
     pch=19, cex=.1, xlab="long", ylab="lat")
plot(x = city.cluster.m[[1]]$long,
     y = city.cluster.m[[1]]$lat,
     col=city.cluster.m[[1]]$cluster_id,
     main = "d = 56 km",
     pch=19, cex=.1, xlab="long", ylab="lat")
dev.off()


###################################################
### code chunk number 18: paper.rnw:238-240
###################################################
  coordinates(city.cluster[[1]]) <- ~long+lat
  proj4string(city.cluster[[1]]) <- CRS("+init=epsg:4326")


###################################################
### code chunk number 19: paper.rnw:252-281 (eval = FALSE)
###################################################
## library(osc)
## data("population")
## 
## tl <- rep(NA,100)
## tm <- rep(NA,100)
## 
## for(i in 1:100){
##   print(i)
##   tm[i] <- system.time(mat <- cca(population, s=i, mode=3))[[1]]
## }
## 
## index <- which(population[]>0)
## long <- ceiling(index/nrow(population))
## lat <- index - ((long-1)*nrow(population))
## popl <- data.frame(long, lat, cluster_id=rep(0, length(lat)))
## 
## for(i in 1:100){
##   print(i)
##   tl[i] <- system.time(list <- cca(popl,s=i) )[[1]]
## }
## 
## saveRDS(tm, "tm.rds")
## saveRDS(tm, "tl.rds")
## 
## pdf("compare.pdf", width=5, height=4)
## plot(tm, xlab="cluster distance", ylab="time in s", col="darkblue", log="", type="l")
## lines(tl, col="darkred")
## legend("topleft", legend=c("matrix","list"), col=c("darkblue","darkred"), pch="-")
## dev.off()