Nothing
#########################################################
# 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)
}
print(paste("use as X-coordinates column '",colnames(data)[1],"'"))
print(paste("use as Y-coordinates column '",colnames(data)[2],"'"))
# 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;
# separate additional data for later
add.data <- NULL
if(ncol(data)>2){
add.data <- data[,3:ncol(data)]
data <- data[,1:2]
}
# 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")
# re-add additional data
if(!is.null(add.data)){mat1 <- cbind(mat1, add.data)}
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)
}
print(paste("use as X-coordinates column '",colnames(data)[1],"'"))
print(paste("use as Y-coordinates column '",colnames(data)[2],"'"))
# preprocessing
n1 <- length(data[,2]);
#data <- ccaR.order(data);
data <- data[order(data[,1]),];
print("Sorting... Done");
# separate additional data
add.data <- NULL
if(ncol(data)>2){
add.data <- data[,3:ncol(data)]
data <- data[,1:2]
}
# 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")
# re-add additional data
if(!is.null(add.data)){mat1 <- cbind(mat1, add.data)}
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=max(dim(input)), return.width=F, complete=F)
{
if(class(input)=="RasterLayer"){
m <- matrix(input[], nrow=input@nrows, byrow=TRUE)
returnraster <- T}else{
m <- input
returnraster <- F
}
if(return.width){
# print(sum(m==0))
if(complete){
# m1 <- .C("ccaBuffEDszS",
# m=as.integer(m),
# nr=as.integer(dim(m)[1]),
# nc=as.integer(dim(m)[2]),
# sz=as.integer(width),
# nz=as.integer(sum(m==0)))
m1 <- .C("ccaBuffEDszNN",
m=as.integer(m),
nr=as.integer(dim(m)[1]),
nc=as.integer(dim(m)[2]),
sz=as.integer(width))
}else{
m1 <- .C("ccaBuffEDszN",
m=as.integer(m),
nr=as.integer(dim(m)[1]),
nc=as.integer(dim(m)[2]),
sz=as.integer(width),
nz=as.integer(sum(m==0)))
}
}else{
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)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.