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
| library(sp)
library(rgdal)
library(raster) ; rasterOptions(tmpdir = 'D://TEMP')
library(data.table)
library(colorspace) # HLS
library(rgl)
library(png)
library(sm) # Density
library(FactoMineR) # ACP
library(MASS)
library(cluster)
library(tree)
library(randomForest)
library(e1071) # SVM
library(caret)
library(leaps)
rasterize.gdal <- function(x, rast, attr) {
grid.extent <- extent(rast)
resolution <- res(rast)
command.base <- "gdal_rasterize -ot Byte -a_nodata 255"
command.attr <- paste("-a", attr)
command.extent <- paste("-te", grid.extent@xmin, grid.extent@ymin, grid.extent@xmax, grid.extent@ymax)
command.resolution <- paste("-tr", res(rast)[1], res(rast)[2])
command.src <- x
command.dst <- paste(x, ".tif", sep = "")
command <- paste(command.base, command.attr, command.extent, command.resolution, command.src, command.dst)
system(command)
raster(command.dst)
}
extract.pixels <- function(x, repository) {
file.name <- basename(x)
# Build search pattern from shapefile name
layername <- sub("\\.[^.]+$", "", file.name)
file.basename <- sub(".CLASS.+$", "-TOA", file.name)
file.pattern <- paste(file.basename, ".*\\.TIF$", sep = "")
# Look for matching rasters
file.images <- list.files(repository, full.names = T, recursive = T, pattern = file.pattern) ; print(file.images)
# Shapefile class attribute name
classes.attr <- "class_code"
# Load rasters
# beginCluster()
layers <- stack(file.images)
classes.vector <- readOGR(f, layername)
classes.raster <- rasterize(f, layers, classes.attr) # Affecte la classe à chaque pixel
# Shorten variable names
names(layers) <- sub(gsub("-", ".", file.basename), "", names(layers))
names(layers) <- sub("^.", "", names(layers))
names(layers) <- sub("\\.", "_", names(layers))
names(layers) <- sub("^S_", "s", names(layers))
names(layers) <- sub("^W_", "w", names(layers))
names(layers) <- sub("^V_", "v", names(layers))
names(layers) <- sub("^T_", "t", names(layers))
names(layers) <- sub("HLS_1", "H", names(layers))
names(layers) <- sub("HLS_2", "L", names(layers))
names(layers) <- sub("HLS_3", "S", names(layers))
names(layers) <- sub("^X[^0-9]", "X", names(layers)) #; print(names(layers))
# Scale reflectance to [0, 1]
layers[["X1"]] <- layers[["X1"]]/1000.0
layers[["X2"]] <- layers[["X2"]]/1000.0
layers[["X3"]] <- layers[["X3"]]/1000.0
layers[["X4"]] <- layers[["X4"]]/1000.0 #; print("Scaled reflectance.")
# Pixel centers
point <- rasterToPoints(classes.raster, spatial = TRUE, progress = "text") #; print("Computed pixel centers.")
names(point@data)[[1]] <- classes.attr
# Load class pixel variable values as point attributes
print("Extracting variable values...")
point <- extract(layers, point, sp = TRUE, progress = "text")
# endCluster()
print("Done.")
rm(layers)
# Retrieve attribute table as a dataframe
point.data <- point@data ; print(paste("Exported dataframe. Length:", nrow(point.data)))
return(point.data)
}
repo <- "/SEAS/02-Images" # Répertoire
pattern <- "-CLASS-RD\\.shp$" # R&D
files <- list.files(repo, full.names = T, recursive = T, pattern = pattern) ; print(files) #, ignore.case = T) # files <- rev(files)
rm(pattern)
rm(dataset)
for (f in files) {
print(f)
if (exists("dataset")){ #print("Extracting new data...")
dataset.new <- extract.pixels(f, repo) #; print("Extracted new data.")
dataset.new <- cbind(dataset.new, basename(f)) ; print(colnames(dataset)) ; print(colnames(dataset.new))
dataset <- rbind(dataset, dataset.new) #; print("Appending new data.")
rm(dataset.new) #; print("Removing temp new data.")
} else { #print("Creating data...")
dataset <- extract.pixels(f, repo) #; print("Created data.")
dataset <- cbind(dataset, basename(f))
}
} |
