-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path.Rhistory
226 lines (226 loc) · 10.2 KB
/
.Rhistory
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
# List of packages to update on R-Studio
# Run the following code:
x <- c("raster","rgdal", "ggplot2", "rasterVis", "ps", "RStoolbox", "dplyr","magrittr","sf", "devtools",
"dplyr", "magrittr", "lubridate", "xts", "tmap", "rgdal","spData", "sp", "gdalUtils", "RColorBrewer", "randomForest",
"RCurl", "rgeos", "maptools", "recipes")
ipak <- function(pkg){
newpkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(newpkg))
install.packages(newpkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
ipak(x)
# install getSpatialData
devtools::install_github("16EAGLE/getSpatialData")
# update other packages: select 1
knitr::opts_chunk$set(echo = TRUE, fig.width=12, fig.height=8, fig.align="center", warning=FALSE, message=FALSE)
library(knitr)
# Information for this tutorial: https://www.neonscience.org/raster-data-r
# Load libraries
library(raster)
ipak <- function(pkg){
newpkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(newpkg))
install.packages(newpkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c("raster","rgdal", "ggplot2", "rasterVis", "ps", "RStoolbox", "dplyr","magrittr","sf", "devtools", "stars", "terra")
ipak(packages)
#setwd("W:\\Ulloa\\SHK\\Lehre\\basic_raster")
w <- setwd("G:\\Documents\\_SHK\\M1")
setwd("C:/Users/ulloa-to/git/Advanced_Remote_sensing_HM/week_1/")
#w <- setwd("G:\\Documents\\_SHK\\M1")
readClipboard()
readClipboard()
readClipboard()
knitr::opts_chunk$set(echo = TRUE, fig.width=12, fig.height=8, fig.align="center", warning=FALSE, message=FALSE)
# Check the projection of the vector
crs(aoi)
# import the corrected boa image
l8_boa_br <- brick("data\\wk2_results\\l8_boa_ref.tif")
detach("package:base", unload = TRUE)
library(base, lib.loc = "C:/Program Files/R/R-4.1.1/library")
library(raster, lib.loc = "C:/Program Files/R/R-4.1.1/library")
# import the corrected boa image
l8_boa_br <- brick("data\\wk2_results\\l8_boa_ref.tif")
library(ipak)
install.packages("ipak")
# Load the packages
lapply(PkgToLoad, library, character.only = TRUE)
readClipboard()
# import the corrected boa image
l8_boa_br <- stack("data\\wk2_results\\l8_boa_ref.tif")
l8_boa_br <- stack("data\\wk2_results\\l8_boa_ref.tif")
readClipboard()
# import the corrected boa image
l8_boa_br <- stack("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM\\data\\wk2_results\\l8_boa_ref.tif")
# import the corrected boa image
l8_boa_br <- stack(paste0(wd, "\\data\\wk2_results\\l8_boa_ref.tif"))
getwd()
# Setting up working directory
wd <- setwd("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM/")
getwd()
# import the corrected boa image
l8_boa_br <- brick(paste0(wd, "\\data\\wk2_results\\l8_boa_ref.tif"))
# import the corrected boa image
l8_boa_br <- brick("\\data\\wk2_results\\l8_boa_ref.tif")
# import the corrected boa image
l8_boa_br <- brick("data\\wk2_results\\l8_boa_ref.tif")
# import the corrected boa image
l8_boa_br <- stack("data\\wk2_results\\l8_boa_ref.tif")
paste0(wd, "data\\wk2_results\\l8_boa_ref.tif" )
# import the corrected boa image
l8_boa_br <- stack(paste0(wd, "data\\wk2_results\\l8_boa_ref.tif"))
# import the corrected boa image
l8_boa_br <- stack("C:/Users/ulloa-to/git/Advanced_Remote_sensing_HM/week_2data\\wk2_results\\l8_boa_ref.tif")
# import the corrected boa image
l8_boa_br <- brick("C:/Users/ulloa-to/git/Advanced_Remote_sensing_HM/week_2data\\wk2_results\\l8_boa_ref.tif")
# import the corrected boa image
l8_boa_br <- brick("C:/Users/ulloa-to/git/Advanced_Remote_sensing_HM/week_2/data\\wk2_results\\l8_boa_ref.tif")
# import the corrected boa image
l8_boa_br <- stack("C:/Users/ulloa-to/git/Advanced_Remote_sensing_HM/week_2/data\\wk2_results\\l8_boa_ref.tif")
# import the corrected boa image
l8_boa_br <- stack("C:/Users/ulloa-to/git/Advanced_Remote_sensing_HM/week_2/data/wk2_results/l8_boa_ref.tif")
# import the corrected boa image
l8_boa_br <- raster("C:/Users/ulloa-to/git/Advanced_Remote_sensing_HM/week_2/data/wk2_results/l8_boa_ref.tif")
# Import meta-data and bands based on MTL file
metaData <- readMeta("data\\LC08_L1TP_174037_20180904_20180912_01_T1\\LC08_L1TP_174037_20180904_20180912_01_T1_MTL.txt")
# Load the packages
lapply(PkgToLoad, library, character.only = TRUE)
# Define the multiple mutiple packages in a vector
PkgToLoad <- c("raster", "rgdal", "ggplot2", "devtools", "sf")
# Load the packages
lapply(PkgToLoad, library, character.only = TRUE)
# import the corrected boa image
l8_boa_br <- raster("C:/Users/ulloa-to/git/Advanced_Remote_sensing_HM/week_2/data/wk2_results/l8_boa_ref.tif")
# Import meta-data and bands based on MTL file
metaData <- readMeta("data\\LC08_L1TP_174037_20180904_20180912_01_T1\\LC08_L1TP_174037_20180904_20180912_01_T1_MTL.txt")
library("RStoolbox")
# Import meta-data and bands based on MTL file
metaData <- readMeta("data\\LC08_L1TP_174037_20180904_20180912_01_T1\\LC08_L1TP_174037_20180904_20180912_01_T1_MTL.txt")
readClipboard()
# Import meta-data and bands based on MTL file
metaData <- readMeta("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM\\data\\LC08_L1TP_174037_20180904_20180912_01_T1\\LC08_L1TP_174037_20180904_20180912_01_T1_MTL.txt")
# Stack the metadata bands
lsMeta <- stackMeta(metaData)
# Correct DN to at-surface-reflecatance with DOS (Chavez decay model)
l8_boa_ref <- radCor(lsMeta, metaData, method = "dos")
writeRaster(l8_boa_ref, datatype="FLT4S", filename = "data\\wk2_results\\l8_boa_ref.tif", format = "GTiff", overwrite=TRUE)
# Setting up working directory
wd <- setwd("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM/")
getwd()
# Setting up working directory
wd <- setwd("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM/")
getwd()
# Setting up working directory
wd <- setwd("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM/")
getwd()
l8_boa_cr
knitr::opts_chunk$set(echo = TRUE, fig.width=12, fig.height=8, fig.align="center", warning=FALSE, message=FALSE)
# Setting up working directory
wd <- setwd("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM/")
getwd()
# Define the multiple mutiple packages in a vector
PkgToLoad <- c("raster", "rgdal", "ggplot2", "devtools", "sf")
# Load the packages
lapply(PkgToLoad, library, character.only = TRUE)
library("RStoolbox")
# Import meta-data and bands based on MTL file
metaData <- readMeta("data\\LC08_L1TP_174037_20180904_20180912_01_T1\\LC08_L1TP_174037_20180904_20180912_01_T1_MTL.txt")
# Stack the metadata bands
lsMeta <- stackMeta(metaData)
# Correct DN to at-surface-reflecatance with DOS (Chavez decay model)
l8_boa_ref <- radCor(lsMeta, metaData, method = "dos")
writeRaster(l8_boa_ref, datatype="FLT4S", filename = "data\\wk2_results\\l8_boa_ref.tif", format = "GTiff", overwrite=TRUE)
# import the corrected boa image
l8_boa_br <- brick("data/wk2_results/l8_boa_ref.tif")
# Load vector shapefile
aoi <- readOGR("data\\vector", "aoi_beirut")
# Check the projection of the vector
crs(aoi)
# Check the projection of your raster
crs(l8_boa_br)
# Reproject aoi if needed
aoi <- spTransform(aoi, crs(l8_boa_br))
# Crop all rasters to the extent of the aoi
l8_boa_cr <- crop(l8_boa_br, aoi)
# plot
plotRGB(l8_boa_cr, r=4, g=3, b=2, axes=TRUE, stretch="lin", main =
"Bottom of Atmosphere Reflectance. Landsat 8")
# Export raster to folder
writeRaster(l8_boa_cr, datatype="FLT4S", filename = "data\\wk2_results\\l8_boa_cr.tif", format = "GTiff", overwrite=TRUE)
# Import the boa_reflectance cropped file
l8_boa_cr <- brick("data\\wk2_results/l8_boa_cr.tif")
# plot
plotRGB(l8_boa_cr, r=4, g=3, b=2, axes=TRUE, stretch="lin", main =
"Bottom of Atmosphere Reflectance. Landsat 8")
# Import the boa_reflectance cropped file
l8_boa_cr <- brick("data\\wk2_results/l8_boa_cr.tif")
# Plot
plotRGB(l8_boa_cr, r=4, g=3, b=2, axes=TRUE, stretch="lin", main =
"Bottom of Atmosphere Reflectance. Landsat 8")
# Plot
plotRGB(l8_boa_cr, r=4, g=3, b=2, axes=TRUE, stretch="lin", main =
"Bottom of Atmosphere Reflectance. Landsat 8")
readClipboard()
readClipboard()
readClipboard()
# set my working directory
w <- setwd("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM")
# set my working directory
w <- setwd("C:\\Users\\ulloa-to\\git\\Advanced_Remote_sensing_HM")
getwd()
libToload <- c("rgdal", "raster")
library(libToload)
library("rgdal")
library("raster")
library("sp")
library("RColorBrewer")
library("RStoolbox")
library(sf)
readClipboard()
s2_pre_june <- stack("data\\s2_change_rasters\\juneCrop_32N.tif")
s2_post_july <- stack("data\\s2_change_rasters\\julyCrop_32N.tif")
readClipboard()
aoi <- readOGR("data\\vector\\wk3_aoi.shp")
crs(aoi)
crs(s2_pre_june)
# reproject the aoi, because it has the wrong projection
aoi_32N <- spTransform(aoi, crs(s2_pre_june))
aoi_32N
# let'S mask both sentinel images to the shape of the vector AOI_32N.
s2_pre_june_mask <- mask(s2_pre_june, aoi_32N)
s2_post_july_mask <- mask(s2_post_july, aoi_32N)
x11()
plotRGB(s2_pre_june, r=3, g=2, b=1, axes=FALSE, stretch="lin", main = "raster cropped")
x11()
plotRGB(s2_pre_june_mask, r=3, g=2, b=1, axes=FALSE, stretch="lin", main = "raster masked")
# create a change raster stack
s2_change <- stack(s2_pre_june_mask, s2_post_july_mask)
# create a change raster stack
s2_pre_june_mask <- resample(s2_pre_june_mask, s2_post_july_mask)
s2_change <- stack(s2_pre_june_mask, s2_post_july_mask)
s2_change
readClipboard()
training_change <- readOGR("data\\vector\\wk3_change_data.shp")
flood_change <- superClass(img = s2_change, nSamples = 1000, trainData = training_change,
responseCol = "class_name")
# reproject our training data
training_change_32N <- spTransform(training_change, crs(s2_change))
training_change_32N
flood_change <- superClass(img = s2_change, nSamples = 1000, trainData = training_change,
responseCol = "class_name")
# let's do the classification
flood_change <- superClass(img = s2_change, nSamples = 1000, trainData = training_change_32N,
responseCol = "class_name")
flood_change$model
x11()
plot(flood_change$map, breaks = c(0, 1, 2, 3, 4, 5, 6),
col = c("darkolivegreen", "tomato", "yellowgreen", "tan", "blue2", "cyan"),
main = "Change Detection")
x11()
plot(flood_change$map, breaks = c(0, 1, 2, 3, 4, 5, 6),
col = c("cyan", "cyan", "tomato", "cyan", "tomato", "cyan"),
main = "Change Detection Binary")
writeRaster(s2_pre_june_mask, datatype="FLT4S", filename = "data/wk3_results/s2_june_32N_mask.tif", format = "GTiff", overwrite=TRUE)
writeRaster(s2_post_july_mask, datatype="FLT4S", filename = "data/wk3_results/s2_july_32N_mask.tif", format = "GTiff", overwrite=TRUE)