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))
}
} |
Partager