library(sp)
library(spatstat)
library(raster)
library(rgeos)
library(rgl)
library(rgdal)
library(ggplot2)
library(arules)
library(rLiDAR)
min=raster("raster-min.tif")
avg=raster("raster-avr.tif")
max=raster("raster-max.tif")height=function(minras,maxras,avgras,r,groundheight=0.1){
avgv=getValues(avgras)
maxv=getValues(maxras)
minv=getValues(minras)
avgv=avgv[avgv>groundheight]
minv=minv[minv>groundheight]
maxv=maxv[maxv>groundheight]
#h=round(r/2,2)
datamax=as.data.frame(table(cut_interval(maxv,length = r)))
datamin=as.data.frame(table(cut_interval(minv,length = r)))
dataavg=as.data.frame(table(cut_interval(avgv,length = r)))
m=max(dim(datamin)[1],dim(datamax)[1],dim(dataavg)[1])
maxfreq=c(datamax$Freq,rep(NA,m-length(datamax$Freq)))
minfreq=c(datamin$Freq,rep(NA,m-length(datamin$Freq)))
avgfreq=c(dataavg$Freq,rep(NA,m-length(dataavg$Freq)))
w=cbind(datamax,avgfreq,minfreq)
colnames(w)[1]="Intervals"
row.names(w)=1:dim(w)[1]
lll=as.numeric(unlist(regmatches(as.character(datamax[1,1]),gregexpr("[0-9]",as.character(datamax[1,1])))))[1]
lowerb=seq(lll,by=r,length.out = dim(w)[1])
upperb=seq(lll+r,by=r,length.out = dim(w)[1])
#newh=seq(h,by=r,length.out = dim(w)[1]+1)
#newh=newh[-1]
w$sum=w$Freq+w$avgfreq+w$minfreq
w$height=upperb#newh
w$lower.bound=lowerb
w$upper.bound=upperb
w[is.na(w)]=0
w
}
w=height(min,max,avg,r=0.5)
head(w,10)## Intervals Freq avgfreq minfreq sum height lower.bound upper.bound
## 1 [7,7.5] 1 1 6 8 7.5 7.0 7.5
## 2 (7.5,8] 21 26 57 104 8.0 7.5 8.0
## 3 (8,8.5] 22 23 50 95 8.5 8.0 8.5
## 4 (8.5,9] 359 370 578 1307 9.0 8.5 9.0
## 5 (9,9.5] 37 37 46 120 9.5 9.0 9.5
## 6 (9.5,10] 28 46 41 115 10.0 9.5 10.0
## 7 (10,10.5] 33 33 46 112 10.5 10.0 10.5
## 8 (10.5,11] 26 44 34 104 11.0 10.5 11.0
## 9 (11,11.5] 39 52 48 139 11.5 11.0 11.5
## 10 (11.5,12] 78 92 92 262 12.0 11.5 12.0buildthewallscoord=function(min,avg,max,wallthreshold=3,roofthreshold=1.5,wallmanypoints=100,roofmanypoints=10,r=0.5){
w=height(min,max,avg,r)
k=which.max(w$sum[!is.na(w$h)])
wr=w[(k+1):nrow(w),]
w=w[1:k,]
w=w[w$sum>=wallmanypoints,]
w=w[complete.cases(w),]
w=w[order(-w$sum),]
k=dim(w)[1]
tw=trunc(wallthreshold/k*100,0)/100
c=data.frame()
ind=0
nn=0
for(i in 1:k){
pol1=rasterToPolygons(min,fun=function(x){x>w$lower.bound[i] & x<=w$upper.bound[i]})
pol2=rasterToPolygons(avg,fun=function(x){x>w$lower.bound[i] & x<=w$upper.bound[i]})
pol3=rasterToPolygons(max,fun=function(x){x>w$lower.bound[i] & x<=w$upper.bound[i]})
poly1=gUnaryUnion(pol1)
poly2=gUnaryUnion(pol2)
poly3=gUnaryUnion(pol3)
poly=union(poly1,poly2)
poly=union(poly,poly3)
poly=gUnaryUnion(poly)
poly=gBoundary(poly)
poly=gSimplify(poly,tol=wallthreshold)
wallthreshold=wallthreshold-tw
c=rbind(c,pcoordinate(poly,w$height[i],indict = ind))
ind=c[dim(c)[1],dim(c)[2]]+1
}
WallRoof=rep("W",dim(c)[1])
wr=wr[wr$sum>=roofmanypoints,]
wr=wr[complete.cases(wr),]
wr=wr[order(-wr$sum),]
k=dim(wr)[1]
tr=trunc(roofthreshold/k*100,0)/100
if(k!=0){
for(i in 1:k){
pol1=rasterToPolygons(min,fun=function(x){x>wr$lower.bound[i] & x<=wr$upper.bound[i]})
pol2=rasterToPolygons(avg,fun=function(x){x>wr$lower.bound[i] & x<=wr$upper.bound[i]})
pol3=rasterToPolygons(max,fun=function(x){x>wr$lower.bound[i] & x<=wr$upper.bound[i]})
poly1=gUnaryUnion(pol1)
poly2=gUnaryUnion(pol2)
poly3=gUnaryUnion(pol3)
poly=union(poly1,poly2)
poly=union(poly,poly3)
poly=gUnaryUnion(poly)
poly=gBoundary(poly)
poly=gSimplify(poly,tol=roofthreshold)
roofthreshold=roofthreshold-tw
c=rbind(c,pcoordinate(poly,wr$height[i],indict = ind))
ind=c[dim(c)[1],dim(c)[2]]+1
}}
WallRoof=c(WallRoof,rep("R",(dim(c)[1]-length(WallRoof))))
c$WallRoof=WallRoof
c
}
pcoordinate=function(x,height.p=0,indict=0){
c=coordinates(x)
l=length(c[[1]])
cc=data.frame()
p=0
m=3
ll=c()
for(j in 1:l){
if(dim(c[[1]][[j]])[1]>m){
cc=rbind(cc,c[[1]][[j]])
p=p+dim(c[[1]][[j]])[1]
ww=rep(indict,dim(c[[1]][[j]])[1])
indict=indict+1
ll=c(ll,ww)
}
}
cc$z=rep(height.p,p)
cc$indicator=ll
cc
}
c=buildthewallscoord(min,avg,max,r=1)
head(c)## x y z indicator WallRoof
## 1 404744 281061.8 27 0 W
## 2 404743 281069.8 27 0 W
## 3 404748 281074.8 27 0 W
## 4 404744 281061.8 27 0 W
## 5 404794 281034.8 27 1 W
## 6 404783 281034.8 27 1 Wtail(c)## x y z indicator WallRoof
## 1556 404772 281032.8 31 241 R
## 1557 404772 281034.8 31 241 R
## 1558 404778 281035.8 31 241 R
## 1559 404778 281034.8 31 241 R
## 1560 404777 281034.8 31 241 R
## 1561 404777 281030.8 31 241 Rbuildthewall3D=function(c,t=10){
l=dim(table(c$indicator))
d=rainbow(l+4)
q=list()
a=c()
s=c()
p=0:(l-1)
while(length(p)>=1 && t>=1){
for(i in p){
a[i+1]=tryCatch({ q[[i+1]]=extrude3d(c[c$indicator==i,1],c[c$indicator==i,2],thickness=c[c$indicator==i,3][1])},error=function(expr) {
return(i)
})
if(is.integer(a[[i+1]])){s=c(s,a[[i+1]])}
}
p=s
s=c()
t=t-1
}
q=q[!sapply(q,is.null)]
z=shapelist3d(q,col=d,plot=F)
z
}
z=buildthewall3D(c)
class(z)## [1] "shapelist3d" "shape3d"shade3d(z)You must enable Javascript to view this page properly.
z=list()
for(i in 1:20){
c=buildthewallscoord(min,avg,max,r = i)
z[[i]]=buildthewall3D(c,t=3)
}mfrow3d(2,3)
for(i in 1:6){
shade3d(z[[i]])
if(i!=6){ next3d() }
}You must enable Javascript to view this page properly.
mfrow3d(2,3)
for(i in 7:12){
shade3d(z[[i]])
if(i!=12){ next3d() }
}You must enable Javascript to view this page properly.
mfrow3d(2,3)
for(i in 12:17){
shade3d(z[[i]])
if(i!=17){ next3d() }
}You must enable Javascript to view this page properly.
mfrow3d(1,3)
for(i in 18:20){
shade3d(z[[i]])
if(i!=20){ next3d() }
}You must enable Javascript to view this page properly.