forked from OpenDroneMap/FIELDimageR
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfct_fieldObject.R
153 lines (153 loc) · 7.37 KB
/
fct_fieldObject.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
#' fieldObject
#'
#' @title Evaluate object area, "x" distance, "y" distance, number, extent, etc.
#'
#' @description Calculating the object dimensions (e.g., area, "x" distance, "y" distance, number, extent, etc.)
#' in the entire mosaic or per plot using the \code{\link{fieldShape}} file.
#' @param mosaic object mask of class stack from the function \code{\link{fieldMask}}.
#' @param fieldShape evaluate the object per plot using the fieldShape as reference.
#' If fieldShape=NULL, the object will be evaluated directly for the entire original image.
#' @param minArea used to set the minimum size percentage of plant canopy (to remove weeds and more).
#' @param areaValue referent value of object area in the image.
#' @param watershed if TRUE the "watershed" algorithm will be used to differentiate objects that are touching each other.
#' @param dissolve if TRUE, polygons with the same attribute value will be dissolved into multi-polygon regions.
#' This option requires the rgeos package.
#' @param perimeter If it is TRUE the object perimeter and the perimeter of the box around the object will be calculated.
#' @param n.rem number of objects that should be removed by decreasing size (n.rem=1 is the background).
#' @param na.rm logical. Should missing values (including NaN) be removed?.
#' @param plot if TRUE the crop image and fieldShape will be plotted.
#'
#' @importFrom raster raster projection values
#' @importFrom methods slot
#' @importFrom stats dist
#' @importFrom spatialEco polyPerimeter
#'
#'@return A list with elements:
#' \itemize{
#' \item \code{mosaic} (cropped by plot)
#' \item \code{Dimension} (area, x.dist, y.dist)
#' \item \code{numObjects} (number of objects)
#' \item \code{Objects} (all objects polygon shape)
#' \item \code{Polygons} (all extent polygons shape)
#' \item \code{single.obj} (single object polygon shape)
#' \item \code{obj.extent} (each object extent)
#' \item \code{x.position} (coordinates of "x" length per object)
#' \item \code{y.position} (coordinates of "y" length per object).
#' }
#'
#'
#' @export
fieldObject <- function(mosaic, fieldShape = NULL, minArea = 0.01, areaValue = 0, watershed = FALSE,
dissolve = TRUE, perimeter = FALSE, n.rem = 1, na.rm = FALSE, plot = TRUE) {
mosaic <- stack(mosaic)
num.band <- length(mosaic@layers)
print(paste(num.band, " layer available", sep = ""))
if (num.band > 1) {
stop("Only mask mosaic with values of 1 and 0 can be evaluated, please use the mask output from fieldMask()")
}
if (!areaValue %in% c(1, 0)) {
stop("The value must be 1 or 0 to represent the objects in the mask mosaic, please use the mask output from fieldMask()")
}
if (na.rm) {
mosaic[is.na(mosaic)] <- c(0, 1)[c(0, 1) != areaValue]
}
print("Identifying objects ... ")
if (is.null(fieldShape)){
r <- stack(mosaic)
fieldShape <- as(extent(r), "SpatialPolygons")
fieldShape <- SpatialPolygonsDataFrame(fieldShape, data.frame(z = 1))
raster::projection(fieldShape) <- raster::projection(r)}
if (plot) {
par(mfrow=c(1,1))
raster::plot(mosaic, col = grey(1:100/100), axes = FALSE,
box = FALSE, legend=FALSE)
sp::plot(fieldShape, add = T)
}
if(watershed){
names(mosaic) <- "mask"
mask <- raster::as.matrix(mosaic$mask) == areaValue
dd <- distmap(mask)
mosaic$watershed <- watershed(dd)
n.obj<-unique(values(mosaic$watershed))[-1]
if (plot) {
par(mfrow=c(1,1))
raster::plot(mosaic$watershed, col = grey((length(n.obj)+10):1/(length(n.obj)+10)), axes = FALSE,
box = FALSE, add=T, legend=FALSE)
sp::plot(fieldShape, add = T)
}}
Out<-list()
numObjects<-NULL
for(s1 in 1:dim(fieldShape)[1]){
single <- fieldShape[s1, ]
CropPlot <- crop(x = mosaic, y = single)
if(watershed){
n.obj <- as.numeric(names(table(values(CropPlot$watershed))[order(table(values(CropPlot$watershed)),decreasing = T)]))[-c(1:n.rem)]
SP <-rasterToPolygons(clump(CropPlot$watershed==n.obj[1]), dissolve = dissolve)
if(length(n.obj)>1){
for(m1 in 2:length(n.obj)){
SP <- rbind(SP,rasterToPolygons(clump(CropPlot$watershed==n.obj[m1]), dissolve = dissolve))
}}
SP_df<- as.data.frame(sapply(slot(SP, "polygons"), function(x) slot(x, "ID")))
row.names(SP_df) <- sapply(slot(SP, "polygons"), function(x) slot(x, "ID"))
colnames(SP_df) <- "ID"
SP <- SpatialPolygonsDataFrame(SP, data =SP_df)
}
if(!watershed){
SP <- rasterToPolygons(clump(CropPlot==areaValue), dissolve = dissolve)
}
sps2 <- lapply(SP@polygons, function(x) SpatialPolygons(list(x)))
print("Taking measurements...")
obj.extent<-list()
x.position<-list()
y.position<-list()
single.obj<-list()
Dimension<-NULL
Polygons<-NULL
for(f1 in 1:length(sps2)){
sps3<-sps2[[f1]]
raster::projection(sps3)<-raster::projection(mosaic)
P<-rasterToPolygons(raster(extent(sps3)),dissolve=TRUE)
area<-raster::area(sps3)
if(area>minArea){
xy1<-cbind(x=c(extent(sps3)[c(1)],
extent(sps3)[c(2)]),
y=c(sum(extent(sps3)[c(3,4)])/2,
sum(extent(sps3)[c(3,4)])/2))
xy2<-cbind(x=c(sum(extent(sps3)[c(1,2)])/2,
sum(extent(sps3)[c(1,2)])/2),
y=c(extent(sps3)[c(3)],
extent(sps3)[c(4)]))
obj.extent[[f1]]<-extent(sps3)
x.position[[f1]]<-xy1
y.position[[f1]]<-xy2
Dimension<-rbind(Dimension,c(area=area,x.dist=dist(xy1),y.dist=dist(xy2)))
single.obj[[f1]]<-sps3
if(!is.null(Polygons)){Polygons<-rbind(Polygons,P);Objects<-rbind(Objects,sps3)}
if(is.null(Polygons)){Polygons<-P;Objects<-sps3}
if (plot) {
raster::plot(P, add=T, border="blue")
lines(xy1,col="red",lty=2)
lines(xy2,col="red",lty=2)}
}
}
Objects_df<- as.data.frame(sapply(slot(Objects, "polygons"), function(x) slot(x, "ID")))
row.names(Objects_df) <- sapply(slot(Objects, "polygons"), function(x) slot(x, "ID"))
colnames(Objects_df) <- "ID"
Objects <- SpatialPolygonsDataFrame(Objects, data =Objects_df)
Polygons_df<- as.data.frame(sapply(slot(Polygons, "polygons"), function(x) slot(x, "ID")))
row.names(Polygons_df) <- sapply(slot(Polygons, "polygons"), function(x) slot(x, "ID"))
colnames(Polygons_df) <- "ID"
Polygons <- SpatialPolygonsDataFrame(Polygons, data =Polygons_df)
print(paste("number of objects on plot ",s1,": ",length(Objects),sep = ""))
if(dim(fieldShape)[1]==1){Out<-list(mosaic=mosaic,Dimension=as.data.frame(Dimension),numObjects=c(numObjects,length(Objects)),Objects=Objects,Polygons=Polygons,single.obj=single.obj,obj.extent=obj.extent,x.position=x.position,y.position=y.position)}
if(dim(fieldShape)[1]!=1){Out[[s1]]<-list(mosaic=CropPlot,Dimension=as.data.frame(Dimension),numObjects=c(numObjects,length(Objects)),Objects=Objects,Polygons=Polygons,single.obj=single.obj,obj.extent=obj.extent,x.position=x.position,y.position=y.position)}
}
if(dim(fieldShape)[1]!=1){names(Out)<-fieldShape$fieldID}
if(dim(fieldShape)[1]!=1){if("PlotName"%in%as.character(names(fieldShape))){names(Out)<-fieldShape$PlotName}}
if(perimeter){
perimeter1 <- polyPerimeter(Out$Objects)
box <- polyPerimeter(Out$Polygons)
Out$Dimension <- cbind.data.frame(Out$Dimension,perimeter=perimeter1,box=box)
}
return(Out)
}