[R-sig-Geo] smothing surface/volume in 3D

Samuel Field fieldsh at mail.med.upenn.edu
Fri Aug 22 17:43:47 CEST 2008


List,

I have a set of locations in three dimensional space.  Each location has a continously scaled attribute and I would like to shade or color the volume between these locations using a simple smoother.  At the moment, I am simply using the cloud() function and coloring and sizing the points based on the value of the attribute, but I would like to think of these locations as sampling from a continously scaled "latent" attribute. I would like something simple to use if possible, since what I need it for is simply illustrative. Although I imagine that shading a volume in 3D isn't very striaght forward! There is an example below.  



#creating data
library(spdep)
library(RColorBrewer)


#sample size
n <- 200

#coordinates

x_coord <- runif(n,0,10)
y_coord <- runif(n,0,10)
z_coord <- runif(n,0,10)


## w matrix and row normalized

w_raw <- matrix(nrow=length(x_coord),ncol=length(x_coord),1)

for(i in 1:length(x_coord)){for(j in 1:length(x_coord)) 
{w_raw[i,j]<- 1/(sqrt((x_coord[i]-x_coord[j])^2 + (y_coord[i]-y_coord[j])^2 + (z_coord[i]-z_coord[j])^2 ))^3}}


diag(w_raw) <- rep(0,n)


row_sum <- rep(1,length(x_coord))
for(i in 1:length(x_coord))  
        {row_sum[i] <- sum(w_raw[i,]) }

w <- matrix(nrow=length(x_coord),ncol=length(x_coord),1)

for(i in 1:length(x_coord)){for(j in 1:length(x_coord)) 
        {w[i,j] <- w_raw[i,j]/row_sum[i]}}

x <- rbinom(n,1,.5)
                     
parms <- c(12.4,0)
w.listw <- mat2listw(w)

e <- rnorm(n,0,1)
p <- .6

y <- parms[1] + parms[2]*x + (solve(diag(n)- p*w))%*%e 

sim.data <- as.data.frame(cbind(y,x,x_coord,y_coord,z_coord))
names(sim.data) <- c("y","x","x_coord","y_coord","z_coord")
error.mod <- errorsarlm(y~x,data = sim.data,w.listw)

summary(error.mod)



brks <- quantile(y, probs=seq(0,1,.20),na.rm=TRUE)
cols <- brewer.pal(length(brks)-1, "YlOrRd")
size <- seq(1:length(cols))/3

cloud(z_coord~x_coord*y_coord,col=cols[findInterval(y,brks,all.inside=TRUE)],
pch=8,cex=size[findInterval(y,brks,all.inside=TRUE)])




More information about the R-sig-Geo mailing list