version 1.0.0

DESCRIPTION

@@ -0,0 +1,15 @@
+Package: osc
+Type: Package
+Title: Orthodromic Spatial Clustering
+Version: 1.0.0
+Date: 2016-05-10
+Depends: raster, R (>= 2.14)
+Suggests: testthat
+Author: Steffen Kriewald, Till Fluschnik, Dominik Reusser, Diego Rybski
+Maintainer: Steffen Kriewald <[email protected]>
+Description: Allows distance based spatial clustering of georeferenced data by implementing the City Clustering Algorithm - CCA. Multiple versions allow clustering for matrix, raster and single coordinates on a plain (euclidean distance) or on a sphere (great-circle or orthodromic distance).
+License: GPL
+NeedsCompilation: yes
+Packaged: 2016-05-10 13:05:35 UTC; kriewald
+Repository: CRAN
+Date/Publication: 2016-05-10 19:05:05

MD5

@@ -0,0 +1,34 @@
+f48a7147e02c0724e05b6358e127ab85 *DESCRIPTION
+33268da2b5a0174d693206835bc49e9e *NAMESPACE
+3b7548d94355d3aac0d7a20f96f8bf3f *R/cca.R
+456b71ba9225142206fdc6d41077d4f2 *R/cca.single.R
+911fcd63d0bc1d8c0ae950319404655d *R/code_revised_cca.R
+ebd670ecc096df9ce8e3f53dcc0419cb *R/getPart.R
+94bdc1f82757c353cb1003e56241cdd6 *build/vignette.rds
+6d7bc5a2d4e9a0646afa0c7e770a26bb *data/datalist
+58fd5c862c5b503e50d5245be0be11d4 *data/exampledata.rda
+e4610df58b13625c9b23bd995a55faff *data/landcover.rda
+4c396c53f51248f850075e93a4d3ddfc *data/population.rda
+971ec966bdc8fc42d3f2eab9a45fc424 *inst/doc/paper.R
+b0f652cbee75143ec695dc018de995b2 *inst/doc/paper.pdf
+bda20f004bf359c3fffe7c064e7dfc40 *inst/doc/paper.rnw
+71ed23ec67d3a47015db2b0b8dc06b07 *man/assign.data.Rd
+f784386850049b09d50f411c8cca098f *man/cca.Rd
+a48d5876310a824b38f8bedd1fe49b96 *man/coordinate.list.Rd
+caf4c194ad1c8880a0991f019062f05f *man/exampledata.Rd
+83dddfbd273fc93fc26ca39690fa42b7 *man/landcover.Rd
+0d846cb42d6a64028ee1bf6b573ce27c *man/osc.buffer.Rd
+3fefbda3fcd0894582ba971665e7c3fc *man/population.Rd
+436aff5a71c41276d06bbd6a31f2b4e6 *src/cca-core.c
+6e3121545e3e3dc8caabba7ba9200906 *src/ccaRevolution.c
+251bdb10958e0b5ac1737af4a05d85e2 *tests/testthat.R
+35f54901471879ab9191912c70ffc422 *tests/testthat/test_cca_multi.R
+1173e8ffffc2fc75fee4a0260c93a919 *tests/testthat/test_cca_single.R
+d8266d70da0e7737622bb663ea664944 *tests/testthat/test_coordinate_list.R
+716806c191b146997b4f6aaa4bd3be08 *vignettes/compare.pdf
+bda20f004bf359c3fffe7c064e7dfc40 *vignettes/paper.rnw
+49979365c9b21512171323f8aebf9c87 *vignettes/pics/exdat1.pdf
+1fdeb0ba794228c0c0a232b901148d0c *vignettes/pics/exdat2.pdf
+8183cf69156264257790404ae1a68a15 *vignettes/pics/fig-concordance.tex
+63dcbf8e81d4a3f795ff4ce9bc7ed9be *vignettes/pics/landcover.pdf
+b273fda9b9218cf2c707d4a4a5bd8b7a *vignettes/pics/raster.pdf

NAMESPACE

@@ -0,0 +1,5 @@
+useDynLib(osc)
+
+export(cca, cca.single, assign.data, osc.buffer, coordinate.list)
+import(raster)
+importFrom("utils", "setTxtProgressBar", "txtProgressBar")

R/cca.R

@@ -0,0 +1,32 @@
+cca <- function(data, s=1, mode=3, count.cells=FALSE, count.max=ncol(data)*3, res.x=NULL, res.y=NULL, cell.class=1, unit="", compare=""){
+	if(class(data) %in% c("data.frame", "RasterLayer")){
+    if(unit=="m"){
+      ccaRm(data=data, d=s,res.x=res.x, res.y=res.y, cell.class=cell.class, compare)
+    } else {
+      ccaRd(data=data, d=s, cell.class=cell.class, compare)
+    }
+	} else {
+		ccaM(data=data, s=s,mode=mode, count.cells=count.cells, count.max=count.max)
+	}
+
+}
+ccaM <- function(data, s=1, mode=3, count.cells=FALSE, count.max=ncol(data)*3){
+	#do checks
+	stopifnot(is.numeric(data))
+	stopifnot(is.matrix(data))
+	stopifnot(is.numeric(s))
+	stopifnot(mode==1 | mode==2 | mode==3)
+	the.data <- as.integer(t(data))
+	clu <- as.integer(rep(0, ncol(data)*nrow(data)))
+        count.max <- as.integer(count.max)
+	count <- as.integer(rep(0, count.max))
+	if(count.cells==TRUE){
+	  out <- .C("callburn_count",  s=as.integer(s), xmax=nrow(data), ymax=ncol(data), mode=as.integer(mode)[1], data=the.data, clu=clu, count=count, count.max=count.max,CLASSES=c("integer", "integer", "integer", "integer", "integer","integer"))
+	  return(list(clusters=matrix(out$clu, ncol=ncol(data), byrow=TRUE), cluster.count=out$count[1:max(out$clu)]))
+	  }
+	if(count.cells==FALSE){
+	  out <- .C("callburn",  s=as.integer(s), xmax=nrow(data), ymax=ncol(data), mode=as.integer(mode)[1], data=the.data, clu=clu, count=count, count.max=count.max,CLASSES=c("integer", "integer", "integer", "integer", "integer","integer"))
+	  return(clusters=matrix(out$clu, ncol=ncol(data), byrow=TRUE))
+	  }
+}
+

R/cca.single.R

@@ -0,0 +1,29 @@
+cca.single <- function(data, s, x,y, mode=3){
+	#do checks
+	stopifnot(is.numeric(data))
+	stopifnot(is.matrix(data))
+	stopifnot(is.numeric(s))
+	stopifnot(is.numeric(x))
+	stopifnot(is.numeric(y))
+	x <- as.integer(x) - 1 # as C starts counting at 0
+	y <- as.integer(y) - 1 # as C starts counting at 0
+	xmax <- nrow(data) 
+	ymax <- ncol(data) 
+	stopifnot(x>=0 & x<xmax) # starting cell must be inside
+	stopifnot(y>=0 & y<ymax) # starting cell must be inside
+	stopifnot(mode==1 | mode==2 | mode==3)
+	method <- switch(mode, "burnn", "burns", "burnr")
+	the.data <- as.integer(t(data))
+	clu <- as.integer(rep(0, ncol(data)*nrow(data)))
+	count.max <- as.integer(ncol(data)*3)
+	count <- as.integer(rep(0, count.max))
+	if(mode==1){
+		out <- .C(method, x=as.integer(x), y=as.integer(y),  c=as.integer(1), xmax=as.integer(xmax), ymax=as.integer(ymax),  data=the.data, clu=clu, CLASSES=c("integer", "integer", "integer", "integer", "integer", "integer","integer"))
+	} else {
+		out <- .C(method,  data=as.integer(the.data), clu=as.integer(clu), x=as.integer(x), y=as.integer(y),  c=as.integer(1), s=as.integer(s), xmax=as.integer(xmax), ymax=as.integer(ymax))
+	}
+
+	return(matrix(out$clu, ncol=ncol(data), byrow=TRUE))
+}
+
+

R/code_revised_cca.R

@@ -0,0 +1,189 @@
+
+#########################################################
+# extract urban coordinates from land-cover data
+#########################################################
+
+coordinate.list <- function(raster, cell.class, compare=""){
+
+  if(compare=="g"){compare <- function(v){return(which(v > cell.class))}}
+  else{ if(compare=="s"){compare <- function(v){return(which(v < cell.class))}}
+        else{compare <- function(v){return(which(v %in% cell.class))}}
+  }
+
+
+  all <- NULL
+  coordinates <- NULL
+
+  bs <- blockSize(raster)
+
+  print("get coordinates:")
+  # create progress bar
+  pb <- txtProgressBar(min = 0, max = bs$n, style = 3)
+
+  for (i in 1:bs$n) {
+      v <- getValues(raster, row=bs$row[i], nrows=bs$nrows[i] )
+      v <- compare(v)
+
+      if(length(v)>0){
+        coordinates <- xyFromCell(object=raster, cell= (raster@ncols * (bs$row[i]-1)) +v)
+      }
+      all <- rbind(all, coordinates)
+      coordinates <- NULL
+      setTxtProgressBar(pb, i)
+    }
+
+  close(pb)
+
+  return(all)
+}
+
+#########################################################
+# cluster urban coordinates
+#########################################################
+
+
+ccaR.order <- function(m) { m <- m[order(-m[,2]),]; m <- m[order(m[,1]),]; return(m); }
+
+ccaRm <- function(data="", d=500, res.x=NULL, res.y=NULL, cell.class, compare)
+{
+
+  # If raster get coordinate list
+  if(class(data)=="RasterLayer"){
+    res.x <- xres(data)
+    res.y <- yres(data)
+    data <- coordinate.list(data, cell.class, compare)
+  }
+
+  # preprocessing
+  n1 <- length(data[,2]);
+  data <- ccaR.order(data);
+  data <- data[order(data[,2]),];
+  print("Sorting... Done");
+  data[,1] <- data[,1]*pi/180;
+  data[,2] <- data[,2]*pi/180;
+
+  # calculate exit condition/distance
+  max_abs <- max(abs(data[,2]));
+  min_dist <- acos(sin(max_abs)*sin(max_abs)+cos(max_abs)*cos(max_abs)*cos(res.y*pi/180))*6371000;
+  min_dist_h <- acos(sin(0)*sin(res.x*pi/180)+cos(0)*cos(res.x*pi/180)*cos(0))*6371000;
+
+  # create place holder for cluster id
+  data <- cbind(data, rep(0,length(data[,1])))
+  data <- as.matrix(data)
+  data <- as.numeric(data[,1:3]);
+
+  step_w <- floor(1+(d/min_dist));
+  step_h <- ceiling(d/min_dist_h);
+  print("Start Clustering...");
+  out <- .C("ccaRevT", m=as.numeric(data), n=as.integer(n1), l=as.numeric(d),
+              step_w=as.integer(step_w), step_h=as.integer(step_h),
+              res_x=as.numeric(res.x*pi/180), res_y=as.numeric(res.y*pi/180),
+              w=as.integer(array(0,n1)));
+  rm(data)
+  print("Clustering... Done");
+  mat1 <- matrix(out$m, ncol=3, byrow=FALSE);
+  m1 <- array(0, max(mat1[,3]));
+  m1 <- .C("ccaSumT", m=as.numeric(mat1), m3=as.integer(mat1[,3]), mm=as.numeric(m1), n=as.integer(n1));
+  print("Summary... Done");
+
+  mat1[,1] <- mat1[,1]*180/pi;
+  mat1[,2] <- mat1[,2]*180/pi;
+  mat1 <- as.data.frame(mat1)
+  colnames(mat1) <- c("long","lat","cluster_id")
+  return(list(cluster=mat1,size=m1$mm));
+}
+
+
+ccaRd <- function(data="", d=500, cell.class, compare)
+{
+
+  # If raster get coordinate list
+  if(class(data)=="RasterLayer"){
+     data <- coordinate.list(data, cell.class, compare)
+  }
+
+  # preprocessing
+  n1 <- length(data[,2]);
+  #data <- ccaR.order(data);
+  data <- data[order(data[,1]),];
+  print("Sorting... Done");
+
+  # create place holder for cluster id
+  data <- cbind(data, rep(0,length(data[,1])))
+  data <- as.matrix(data)
+  data <- as.numeric(data[,1:3]);
+
+  print("Start Clustering...");
+  out <- .C("ccaRev", m=as.numeric(data), n=as.integer(n1), l=as.numeric(d),w=as.integer(array(0,n1)));
+  rm(data)
+  print("Clustering... Done");
+  mat1 <- matrix(out$m, ncol=3, byrow=FALSE);
+  m1 <- array(0, max(mat1[,3]));
+  m1 <- .C("ccaSum", m=as.numeric(mat1), m3=as.integer(mat1[,3]), mm=as.numeric(m1), n=as.integer(n1));
+  print("Summary... Done");
+
+  mat1 <- as.data.frame(mat1)
+  colnames(mat1) <- c("long","lat","cluster_id")
+  return(list(cluster=mat1,size=m1$mm));
+}
+
+#########################################################
+# population assignment to clusters 
+#########################################################
+
+assign.data <- function(cluster, points, dist=1000){
+
+  # conversion from meter to degree
+  wgs84 <- 6371000
+  scale <- 360 / (2*pi*wgs84)
+  y_dist <- dist * scale
+
+  points <- cbind(points, cluster_id=rep(0, length(points[, 1])))
+
+  pop_cluster <- NULL
+
+  for(j in 1:length(points[, 1])){
+    #print(j)
+    x <- points$long[j]
+    y <- points$lat[j]
+    x_dist <- dist * scale * (1 / cos( y*pi / 180 ))
+    temp <- which((x-x_dist) < cluster[, 1] &
+                    cluster[, 1] < (x+x_dist) &
+                    (y-y_dist) < cluster[, 2] &
+                    cluster[, 2] < (y+y_dist))
+    if(length(temp)>0){
+      pdist <- pointDistance(c(x, y), matrix(c(cluster[temp, 1],cluster[temp, 2]),
+                                             ncol=2, byrow=FALSE), longlat=TRUE)
+      if(min(pdist)< dist){
+        index <- which(pdist[]==min(pdist))
+        cluster_id <- cluster[temp[index], 3]
+        points$cluster_id[j] <- cluster_id
+        }
+    }
+  }
+  return(points)
+}
+
+#########################################################
+# buffer
+#########################################################
+
+osc.buffer <- function(input, width)
+{	
+  if(class(input)=="RasterLayer"){
+  m <- matrix(input[], nrow=input@nrows, byrow=TRUE)
+  returnraster <- T}else{
+    m <- input
+    returnraster <- F
+  }
+	m1 <- .C("ccaBuffED", m=as.integer(m), nr=as.integer(dim(m)[1]), nc=as.integer(dim(m)[2]), sz=as.integer(width))
+	m1 <- m1$m
+#	m1[which(m1<0)] <- -1
+	m1 <- matrix(m1, nrow=dim(m)[1])
+	if(returnraster){
+	  input[] <- raster(m1)[]
+	  return(input)
+	  }else{
+	    return(m1)
+	  }
+}

R/getPart.R

@@ -0,0 +1,38 @@
+cca.row <- function(pop, row){
+	stopifnot(is.numeric(pop))
+	stopifnot(is.matrix(pop))
+	stopifnot(is.numeric(row))
+	stopifnot(row<=NROW(pop))
+	stopifnot(row>=0)
+	the.pop <- as.integer(t(pop))
+	clu <- rep(999, ncol(pop))
+	out <- .C("getrow",  col=as.integer(row-1), xmax=nrow(pop), ymax=ncol(pop),  pop=as.integer(the.pop), ret=as.integer(clu))
+
+	return(out$ret)
+}
+cca.col <- function(pop, col){
+	stopifnot(is.numeric(pop))
+	stopifnot(is.matrix(pop))
+	stopifnot(is.numeric(col))
+	stopifnot(col<=NCOL(pop))
+	stopifnot(col>=1)
+	the.pop <- as.integer(t(pop))
+	clu <- rep(999, nrow(pop))
+	out <- .C("getcol",  col=as.integer(col-1), xmax=nrow(pop), ymax=ncol(pop),  pop=as.integer(the.pop), ret=as.integer(clu))
+
+	return(out$ret)
+}
+#cca.block <- function(pop, row,col, drow,dcol){
+#	stopifnot(is.numeric(pop))
+#	stopifnot(is.matrix(pop))
+#	stopifnot(is.numeric(row))
+#	stopifnot(row+drow<=NROW(pop))
+#	stopifnot(col+dcol<=NROW(pop))
+#	stopifnot(row>=1)
+#	stopifnot(col>=1)
+#	the.pop <- as.integer(t(pop))
+#	clu <- rep(999, dcol*drow)
+#	out <- .C("getblock",  row=as.integer(row-1), col=as.integer(col-1),drow=as.integer(drow),dcol=as.integer(dcol),xmax=nrow(pop), ymax=ncol(pop),  pop=the.pop, ret=as.integer(clu))
+#
+#	return(matrix(out$ret, ncol=dcol))
+#}

data/datalist

@@ -0,0 +1,3 @@
+exampledata
+landcover
+population