library(sp)
library(spatstat)
library(raster)
library(rgeos)
library(rgl)
library(rgdal)
library(ggplot2)
library(arules)
library(rLiDAR)
main=readLAS("New folder/main.las")
head(main)
## X Y Z Intensity ReturnNumber
## [1,] 404761.0 281158.3 132.38 5 2
## [2,] 404761.3 281158.6 141.07 33 1
## [3,] 404761.1 281157.7 132.60 26 2
## [4,] 404760.5 281158.4 132.44 28 1
## [5,] 404758.8 281158.3 132.45 40 1
## [6,] 404759.4 281157.6 132.48 26 1
e <- extent(main[,1:2])
e=raster(e,ncol=300,nrow=200)
avg=rasterize(main[,1:2],e,main[,3],fun=mean)
int=rasterize(main[,1:2],e,main[,4],fun=mean)
par(mfrow=c(1,2))
plot(avg)
plot(int)
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(strsplit(strsplit(as.character(datamax[1,1]),"\\[")[[1]][2],",")[[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
}
delground=function(main,r=0.5,step=1){
e <- extent(main[,1:2])
e=raster(e,ncol=300,nrow=200)
avg=rasterize(main[,1:2],e,main[,3],fun=mean)
w=height(avg,avg,avg,r)
if(w$Freq[1]>w$Freq[2]){
ind=1
} else {
q=diff(diff(w$Freq)>=0)<0
indq=which(q==TRUE)
if(length(indq)==0){return("Please Choose Smaller r")}
i=which(w$Freq[diff(diff(w$Freq)>=0)<0]==max(w$Freq[diff(diff(w$Freq)>=0)<0]))
ind=indq[i]+step}
h=w$upper.bound[ind]
main=as.data.frame(main)
main=main[main$Z>h,]
main
}
head(height(avg,avg,avg,r=0.5),20)
## Intervals Freq avgfreq minfreq sum height lower.bound upper.bound
## 1 [126,126.5] 1 1 1 3 126.5 126.0 126.5
## 2 (126.5,127] 0 0 0 0 127.0 126.5 127.0
## 3 (127,127.5] 1 1 1 3 127.5 127.0 127.5
## 4 (127.5,128] 36 36 36 108 128.0 127.5 128.0
## 5 (128,128.5] 433 433 433 1299 128.5 128.0 128.5
## 6 (128.5,129] 370 370 370 1110 129.0 128.5 129.0
## 7 (129,129.5] 515 515 515 1545 129.5 129.0 129.5
## 8 (129.5,130] 1417 1417 1417 4251 130.0 129.5 130.0
## 9 (130,130.5] 572 572 572 1716 130.5 130.0 130.5
## 10 (130.5,131] 718 718 718 2154 131.0 130.5 131.0
## 11 (131,131.5] 454 454 454 1362 131.5 131.0 131.5
## 12 (131.5,132] 2024 2024 2024 6072 132.0 131.5 132.0
## 13 (132,132.5] 5048 5048 5048 15144 132.5 132.0 132.5
## 14 (132.5,133] 7867 7867 7867 23601 133.0 132.5 133.0
## 15 (133,133.5] 1202 1202 1202 3606 133.5 133.0 133.5
## 16 (133.5,134] 702 702 702 2106 134.0 133.5 134.0
## 17 (134,134.5] 918 918 918 2754 134.5 134.0 134.5
## 18 (134.5,135] 813 813 813 2439 135.0 134.5 135.0
## 19 (135,135.5] 737 737 737 2211 135.5 135.0 135.5
## 20 (135.5,136] 740 740 740 2220 136.0 135.5 136.0
main=as.matrix(delground(main))
main1=as.matrix(delground(main))
main2=as.matrix(delground(main,step=2))
main3=as.matrix(delground(main,step=4))
par(mfrow=c(1,3))
e <- extent(main1[,1:2])
e=raster(e,ncol=300,nrow=200)
avg=rasterize(main1[,1:2],e,main1[,3],fun=mean)
plot(avg)
e <- extent(main2[,1:2])
e=raster(e,ncol=300,nrow=200)
avg=rasterize(main2[,1:2],e,main2[,3],fun=mean)
plot(avg)
e <- extent(main3[,1:2])
e=raster(e,ncol=300,nrow=200)
avg=rasterize(main3[,1:2],e,main3[,3],fun=mean)
plot(avg)
