-
Notifications
You must be signed in to change notification settings - Fork 0
/
index.Rmd
338 lines (257 loc) · 17.3 KB
/
index.Rmd
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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
---
title: "Monitoring vegetation change for over 30 years"
author: "Ke Li"
output: html_document
editor_options:
chunk_output_type: inline
---
```{r setup, include=FALSE,message=FALSE}
library(knitr)
knitr::opts_chunk$set(cache = TRUE)
```
# Introduction
Knowing change of vegetation plays a significant role in city management and environment issue. Image classification combining with remote sensing satellite data provides an promosing way to display mapping of vegetation change. However, assigning one class in one pixel which has 30 meter sptial resolution always accompany with some classification mistake, citing the reason that one pixel always mixed with different classes. In this project, I intend to use multiple endmember spectral mixture analysis to acquire proportion of each class for per-pixel level. Moreover, to understand change of vegetation, supervised classification method will be used in unmixing image to get the trend of vegetation change.
# Materials and methods
Since vegetation change is focus on change between 30 years, Landsat 8 OLI and Landsat 5 August images will be used in this project to collecting enough pure pixels of each class and finding training samples of vegetation increase, vegetation decrease and vegetation no change class.There are three main steps in this project:
1. The pure pixel of vegetation, urban area, soil and water was collected in 1988, 1998, 2008 and 2018 using ROI method in ENVI. The location of each sample were recorded and export as .shp file. The mean value of each extract satellite value of each class is considered as the spectrum of pure pixel.
2. Acquiring fraction of each class within one-pixel level after performing multiple endmemer spectral mixture analysis, the change images between 1988 and 1998, 1998 and 2008, 2008 and 2018 were acquired using subtraction method. Therefore, each pixel contains change information about vegetation, urban, soil and water.
3. ROI method in ENVI in each change year was used in collecting change pixel about vegetation increase, vegetation decrease, vegetation no change and other change. Change points of vegetation increase, vegetation decrease, vegetation no change are collected about each 300 points, the number of point of other change points are 150 points. I use 70% samples to train maximum likelihood classification model, and 30% samples to test classification method.
## Loading package
The packages used in this project are listed below:
```{r, message=FALSE,'load'}
library(googledrive)
library(raster)
library(rgdal)
library(sf)
library(foreach)
library(doParallel)
library(tidyverse)
library(RStoolbox)
```
```{r, echo=FALSE,message=F,include=FALSE,'import'}
img2018=brick("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/ImageBuffalo2017-2019.tif")
img2018_water=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/2018water.shp")
img2018_veg=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/2018veg.shp")
img2018_urban=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/2018urban.shp")
img2018_soil=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/2018soil.shp")
img2008=brick("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/ImageBuffalo2007-2009.tif")
img1998=brick("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/ImageBuffalo1997-1999.tif")
img1988=brick("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/ImageBuffalo1987-1989.tif")
img1988_water=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/1988water.shp")
img1988_veg=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/1988veg.shp")
img1988_urban=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/1988urban.shp")
img1988_soil=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/1988soil.shp")
img1998_water=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/1998water.shp")
img1998_veg=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/1998veg.shp")
img1998_urban=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/1998urban.shp")
img1998_soil=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/1998soil.shp")
img2008_water=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/2008water.shp")
img2008_veg=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/2008veg.shp")
img2008_urban=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/2008urban.shp")
img2008_soil=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/2008soil.shp")
buffalo<-st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/buffalo/buffalo_shp.shp")
repro_buffalo=st_transform(buffalo,crs=crs(img2018))
buffalo_lonlat<-as.data.frame(st_coordinates(repro_buffalo$geometry))[,c(1,2)]
id=data.frame(name="buffalo")
poly=Polygons(list(Polygon(buffalo_lonlat)),1)
buffalo_poly<-SpatialPolygons(list(poly))
buffalo_frame<-SpatialPolygonsDataFrame(buffalo_poly,id)
crs(buffalo_frame)<-crs(img2018)
buffalo_img=brick("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/buffalo_img.png")
```
## Study area
The post image is Landsat8OLI image I used in 2018.
```{r, echo=FALSE,message=FALSE}
plotRGB(buffalo_img,1,2,3)
```
## Extract sample value
After extract pure pixel satellite value, the value of endmember of each class were averaged as the value of each pure class.
```{r, message=FALSE,'extract'}
extract_water_2018=raster::extract(img2018,img2018_water)
extract_veg_2018=raster::extract(img2018,img2018_veg)
extract_urban_2018=raster::extract(img2018,img2018_urban)
extract_soil_2018=raster::extract(img2018,img2018_soil)
water_2018<-colMeans(foreach(i=1:NROW(extract_water_2018),.combine = "rbind")%do%
unlist(extract_water_2018[[i]]))
soil_2018<-colMeans(foreach(i=1:NROW(extract_soil_2018),.combine = "rbind")%do%
unlist(extract_soil_2018[[i]]))
urban_2018<-colMeans(foreach(i=1:NROW(extract_urban_2018),.combine = "rbind")%do%
unlist(extract_urban_2018[[i]]))
veg_2018<-colMeans(foreach(i=1:NROW(extract_veg_2018),.combine = "rbind")%do%
unlist(extract_veg_2018[[i]]))
end_2018<-rbind(water_2018,soil_2018,urban_2018,veg_2018)
rownames(end_2018)<-c('water',"soil",'urban','veg')
end_2018
```
```{r,echo=FALSE,message=FALSE,include=FALSE,'extract1'}
extract_water_1988=raster::extract(img1988,img1988_water)
extract_veg_1988=raster::extract(img1988,img1988_veg)
extract_urban_1988=raster::extract(img1988,img1988_urban)
extract_soil_1988=raster::extract(img1988,img1988_soil)
extract_water_1998=raster::extract(img1988,img1998_water)
extract_veg_1998=raster::extract(img1988,img1998_veg)
extract_urban_1998=raster::extract(img1988,img1998_urban)
extract_soil_1998=raster::extract(img1988,img1998_soil)
extract_water_2008=raster::extract(img2008,img2008_water)
extract_veg_2008=raster::extract(img2008,img2008_veg)
extract_urban_2008=raster::extract(img2008,img2008_urban)
extract_soil_2008=raster::extract(img2008,img2008_soil)
water_1988<-colMeans(foreach(i=1:NROW(extract_water_1988),.combine = "rbind")%do%
unlist(extract_water_1988[[i]]))
soil_1988<-colMeans(foreach(i=1:NROW(extract_soil_1988),.combine = "rbind")%do%
unlist(extract_soil_1988[[i]]))
urban_1988<-colMeans(foreach(i=1:NROW(extract_urban_1988),.combine = "rbind")%do%
unlist(extract_urban_1988[[i]]))
veg_1988<-colMeans(foreach(i=1:NROW(extract_veg_1988),.combine = "rbind")%do%
unlist(extract_veg_1988[[i]]))
water_1998<-colMeans(foreach(i=1:NROW(extract_water_1998),.combine = "rbind")%do%
unlist(extract_water_1998[[i]]))
soil_1998<-colMeans(foreach(i=1:NROW(extract_soil_1998),.combine = "rbind")%do%
unlist(extract_soil_1998[[i]]))
urban_1998<-colMeans(foreach(i=1:NROW(extract_urban_1998),.combine = "rbind")%do%
unlist(extract_urban_1998[[i]]))
veg_1998<-colMeans(foreach(i=1:NROW(extract_veg_1998),.combine = "rbind")%do%
unlist(extract_veg_1998[[i]]))
water_2008<-colMeans(foreach(i=1:NROW(extract_water_2008),.combine = "rbind")%do%
unlist(extract_water_2008[[i]]))
soil_2008<-colMeans(foreach(i=1:NROW(extract_soil_2008),.combine = "rbind")%do%
unlist(extract_soil_2008[[i]]))
urban_2008<-colMeans(foreach(i=1:NROW(extract_urban_2008),.combine = "rbind")%do%
unlist(extract_urban_2008[[i]]))
veg_2008<-colMeans(foreach(i=1:NROW(extract_veg_2008),.combine = "rbind")%do%
unlist(extract_veg_2008[[i]]))
end_1988<-rbind(water_1988,soil_1988,urban_1988,veg_1988)
rownames(end_1988)<-c('water',"soil",'urban','veg')
end_1998<-rbind(water_1998,soil_1998,urban_1998,veg_1998)
rownames(end_1998)<-c('water',"soil",'urban','veg')
end_2008<-rbind(water_2008,soil_2008,urban_2008,veg_2008)
rownames(end_2008)<-c('water',"soil",'urban','veg')
vegincrease2018=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/2018change/2018vegincrease.shp")
vegdecrease2018=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/2018change/2018vegdecrease.shp")
other2018=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/2018change/2018other.shp")
vegnochange2018=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/2018change/2018vegnochange.shp")
```
## Unmixing image using Multiple Endmember Spectral Mixture Analysis
Multiple Endmember Spectral Mixture Anslysis was used in acquiring fraction of each class. By allowing various endmember for per-pixel level, each pixel can be accounted for fraction of each class. In this study, MESMA acquired five imgae: soil fraction, vegetation fraction, urban fracion and water fraction and RMSE. Considering the change fraction of four class, the first four layer were extracted.
```{r}
unmix2018 <- mesma(img2018,end_2018, method = "NNLS")
unmix2008 <- mesma(img2008,end_2008, method = "NNLS")
change_2018<- (unmix2018-unmix2008)[[1:4]]
```
```{r,echo=FALSE,message=FALSE,include=FALSE}
unmix1988 <- mesma(img1988,end_1988, method = "NNLS")
unmix1998 <- mesma(img1998,end_1998, method = "NNLS")
change_1998<- (unmix1998-unmix1988)[[1:4]]
change_2008<- (unmix2008-unmix1998)[[1:4]]
vegincrease1988=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/1988vegincrease.shp")
vegdecrease1988=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/1988vegdecrease.shp")
other1988=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/1988other.shp")
vegnochange1988=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/1988vegnochange.shp")
vegincrease2008=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/2008change/2008vegincrease.shp")
vegdecrease2008=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/2008change/2008vegdecrease.shp")
other2008=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/2008change/2008other.shp")
vegnochange2008=st_read("C:/Users/like9/Documents/GitHub/2019-geo511-project-kli57/data/change/2008change/2008vegnochange.shp")
```
## Monitor vegetation change using Maximum likelihood classification
The collected shapefile data is not spatialpoint data frame, so I create data frame using imported data point and connect to the longitude and latitude. And Maximum likelihood to classify the four class: vegetation increase, vegetation decrease, vegetation nochange and other.
```{r,message=FALSE}
class_2018=data.frame(class=c(rep('vegincrease2018',length(vegincrease2018$ID)),
rep('vegdecrease2018',length(vegdecrease2018$ID)),
rep('vegnochange2018',length(vegnochange2018$ID)),
rep('other2018',length(other2018$ID))))
vegin_lonlat2018<-as.data.frame(st_coordinates(vegincrease2018$geometry))
vegde_lonlat2018<-as.data.frame(st_coordinates(vegdecrease2018$geometry))
vegno_lonlat2018<-as.data.frame(st_coordinates(vegnochange2018$geometry))
other_lonlat2018<-as.data.frame(st_coordinates(other2018$geometry))
lon_lat_2018<-rbind(vegin_lonlat2018,vegde_lonlat2018,vegno_lonlat2018,other_lonlat2018)
training2018<-SpatialPointsDataFrame(lon_lat_2018,class_2018)
crs(training2018)<-crs(change_2018)
mlh_2018<-superClass(change_2018, trainData = training2018, responseCol = "class",
model="mlc", tuneLength = 1, trainPartition = 0.7)
```
```{r,echo=FALSE,include=FALSE}
class_1988=data.frame(class=c(rep('vegincrease1988',length(vegincrease1988$ID)),
rep('vegdecrease1988',length(vegdecrease1988$ID)),
rep('vegnochange1988',length(vegnochange1988$ID)),
rep('other1988',length(other1988$ID))))
class_2008=data.frame(class=c(rep('vegincrease2008',length(vegincrease2008$ID)),
rep('vegdecrease2008',length(vegdecrease2008$ID)),
rep('vegnochange2008',length(vegnochange2008$ID)),
rep('other2008',length(other2008$ID))))
vegin_lonlat<-as.data.frame(st_coordinates(vegincrease1988$geometry))
vegde_lonlat<-as.data.frame(st_coordinates(vegdecrease1988$geometry))
vegno_lonlat<-as.data.frame(st_coordinates(vegnochange1988$geometry))
other_lonlat<-as.data.frame(st_coordinates(other1988$geometry))
vegin_lonlat2008<-as.data.frame(st_coordinates(vegincrease2008$geometry))
vegde_lonlat2008<-as.data.frame(st_coordinates(vegdecrease2008$geometry))
vegno_lonlat2008<-as.data.frame(st_coordinates(vegnochange2008$geometry))
other_lonlat2008<-as.data.frame(st_coordinates(other2008$geometry))
lon_lat_2008<-rbind(vegin_lonlat2008,vegde_lonlat2008,vegno_lonlat2008,other_lonlat2008)
training2008<-SpatialPointsDataFrame(lon_lat_2008,class_2008)
crs(training2008)<-crs(change_2008)
mlh_2008<-superClass(change_2008, trainData = training2008, responseCol = "class",
model="mlc", tuneLength = 1, trainPartition = 0.7)
lon_lat_1988<-rbind(vegin_lonlat,vegde_lonlat,vegno_lonlat,other_lonlat)
training1988<-SpatialPointsDataFrame(lon_lat_1988,class_1988)
crs(training1988)<-crs(change_1998)
mlh_1988<-superClass(change_1998, trainData = training1988, responseCol = "class",
model="mlc", tuneLength = 1, trainPartition = 0.7)
mlh1998_mask<-mask(mlh_1988$map,buffalo_frame)
mlh2008_mask<-mask(mlh_2008$map,buffalo_frame)
mlh2018_mask<-mask(mlh_2018$map,buffalo_frame)
```
# Results
```{r,echo=FALSE,message=FALSE}
#change classification between 1988 and 1998
mlh_1988$validation$performance$table
mlh_1988$validation$performance$overall
#change classification between 2008 and 1998
mlh_2008$validation$performance$table
mlh_2008$validation$performance$overall
#change classification between 2018 and 2008
mlh_2018$validation$performance$table
mlh_2018$validation$performance$overall
```
According to the confusion matrix of three change images, we can see that the classification accuracy for maximum likelihood method is higher.
```{r,echo=FALSE,fig.width=9, fig.height=6,message=FALSE}
ext=extent(buffalo_frame)
par(mfrow=c(1,2))
par(xpd=FALSE)
plot(mlh1998_mask,ext=extent(buffalo_frame),
legend=F,
col=c('white','red','green','grey'),
main="change image (1998 and 1988)")
par(xpd=T)
legend('bottomleft',legend = c('other','vegdecrease','vegincrease','vegnochange'),
fill = c('white','red','green','grey'),
bty='n',cex = 0.8)
par(xpd=FALSE)
plot(mlh2008_mask,ext=extent(buffalo_frame),
legend=FALSE,col=c('white','red','green','grey'),
main='change image (2008 and 1998)')
par(xpd=T)
legend('bottomleft',legend = c('other','vegdecrease','vegincrease','vegnochange'),
fill = c('white','red','green','grey'),
bty='n',cex = 0.8)
par(xpd=FALSE)
plot(mlh2018_mask,ext=extent(buffalo_frame),
legend=FALSE,col=c('white','red','green','grey'),
main='change image (2018 and 2008)')
par(xpd=T)
legend('bottomleft',legend = c('other','vegdecrease','vegincrease','vegnochange'),
fill = c('white','red','green','grey'),
bty='n',cex = 0.8)
```
```{r,echo=FALSE,message=FALSE,include=FALSE,'calculate'}
fre_bind=cbind(freq(mlh1998_mask)[2:4,],freq(mlh2008_mask)[2:4,],freq(mlh2018_mask)[2:4,])[,c(2,4,6)]
as.data.frame(fre_bind)
row.names(fre_bind)=c('vegdecrease','vegincrease','vegnochange')
colnames(fre_bind)=c('1998','2008','2018')
```
```{r}
fre_bind
```
# Conclusions
This study acquires fraction of each class for per-pixel level in 1988, 1998, 2008 and 2018, and use maximum likelihood classifier to monitor the change of vegetation over 30 years. The classification accuracy of change image between 1998 and 1988, 2008 and 1998,2008 and 2018 have higher accuracy from the testing data. From the vegetation trend for three chang image, the vegetation increase is higher than the vegetation decrease. The area of vegetation no change in 2018 is higher than 2008 and 1998. All in all, the coverage of vegetation is increasing and expand over years.
# References
Rogan, John, Janet Franklin, and Dar A. Roberts. "A comparison of methods for monitoring multitemporal vegetation change using Thematic Mapper imagery." Remote Sensing of Environment 80.1 (2002): 143-156.