https://github.com/cran/plotKML
Tip revision: 7e79f0094a9ff898febe51089ecf40eb542c523e authored by Tomislav Hengl on 09 March 2020, 21:40:02 UTC
version 0.6-1
version 0.6-1
Tip revision: 7e79f00
reproject.R
reproject.SpatialPoints <- function(obj, CRS = get("ref_CRS", envir = plotKML.opts), ...) {
message(paste("Reprojecting to", CRS, "..."))
res <- spTransform(x = obj, CRSobj = CRS(CRS))
return(res)
}
reproject.RasterLayer <- function(obj, CRS = get("ref_CRS", envir = plotKML.opts), program = "raster", tmp.file = TRUE, NAflag = get("NAflag", envir = plotKML.opts), show.output.on.console = FALSE, method, ...) {
if(raster::is.factor(obj)){
method <- "ngb"
} else {
if(missing(method)){ method <- "bilinear" }
}
if(program=="raster"){
message(paste("Reprojecting to", CRS, "..."))
res <- raster::projectRaster(obj, crs = CRS, method = method, progress='text')
names(res) <- names(obj)
} else {
if(program=="GDAL"){
gdalwarp <- get("gdalwarp", envir = plotKML.opts)
# look for GDAL path:
if(nchar(gdalwarp)==0){
plotKML.env(show.env = FALSE)
gdalwarp <- get("gdalwarp", envir = plotKML.opts)
}
if(tmp.file==TRUE){
tf <- tempfile()
} else {
tf <- normalizeFilename(deparse(substitute(obj, env = parent.frame())))
}
if(!nchar(gdalwarp)==0){
if(method == "ngb") { method <- "near" }
writeRaster(obj, paste(tf, ".tif", sep=""), overwrite=TRUE, NAflag=NAflag)
# resample to WGS84 system:
message(paste("Using gdalwarp function:", gdalwarp))
message(paste("Reprojecting to", CRS, "..."))
system(paste(gdalwarp, " ", tf, ".tif", " -t_srs \"", CRS, "\" ", tf, "_ll.tif -dstnodata \"", NAflag, "\" ", " -r ", method, sep=""), show.output.on.console = show.output.on.console)
res <- raster(paste(tf, "_ll.tif", sep=""), silent = TRUE)
names(res) <- names(obj)
} else {
stop("Could not locate GDAL. See 'plotKML.env()' for more info.") }
}
}
return(res)
}
reproject.SpatialGrid <- function(obj, CRS = get("ref_CRS", envir = plotKML.opts), tmp.file = TRUE, program = "raster", NAflag = get("NAflag", envir = plotKML.opts), show.output.on.console = FALSE, ...) {
## convert all character vectors to factors:
for(j in 1:ncol(obj)){
if(is.character(obj@data[,j])){ obj@data[,j] <- as.factor(obj@data[,j]) }
}
if(program=="raster"){
# if multiple layers:
if(ncol(obj) > 1) {
r <- raster::stack(obj)
res <- list(NULL)
for(j in 1:ncol(obj)){
if(is.factor(obj@data[,j])){
r[[j]] <- raster::as.factor(r[[j]])
}
res[[j]] <- reproject(r[[j]], CRS = CRS)
}
res <- as(raster::stack(res), "SpatialGridDataFrame")
## TH: time consuming but would be preferred:
# res <- as(res, "SpatialPixelsDataFrame")
names(res) <- names(obj)
}
# single layer:
else {
r <- raster(obj)
if(is.factor(obj@data[,1])){
r <- raster::as.factor(r)
message(paste("Reprojecting to", CRS, "..."))
res <- as(raster::projectRaster(r, crs = CRS, method = "ngb"), "SpatialGridDataFrame")
} else {
res <- as(reproject(r, CRS = CRS), "SpatialGridDataFrame")
}
names(res) <- names(obj)
}
# try to fix factor-type objects:
for(j in 1:ncol(obj)){
if(is.factor(obj@data[,j])){
# copy levels:
res@data[,j] <- as.factor(res@data[,j])
try( levels(res@data[,j]) <- levels(obj@data[,j]), silent = TRUE )
}
}
}
if(program=="GDAL"){
gdalwarp <- get("gdalwarp", envir = plotKML.opts)
# look for GDAL path if missing:
if(nchar(gdalwarp)==0){
plotKML.env()
gdalwarp <- get("gdalwarp", envir = plotKML.opts)
}
message(paste("Using gdalwarp function:", gdalwarp))
if(!nchar(gdalwarp)==0){
for(j in 1:ncol(obj)){
if(tmp.file==TRUE){
tf <- tempfile()
}
else {
tf <- paste(normalizeFilename(deparse(substitute(obj, env = parent.frame()))), names(obj)[j], sep="_")
}
# write SPDF to a file:
if(is.factor(obj@data[,j])){
x <- obj[j]
x@data[,1] <- as.integer(x@data[,1])
if(max(x@data[,1], na.rm=TRUE)<254){
writeGDAL(x, paste(tf, ".tif", sep=""), "GTiff", mvFlag = 255, type="Byte")
} else {
writeGDAL(x, paste(tf, ".tif", sep=""), "GTiff")
}
}
else {
writeGDAL(obj[j], paste(tf, ".tif", sep=""), "GTiff")
}
message(paste("Reprojecting to", CRS, "..."))
# resample to WGS84 system:
if(is.factor(obj@data[,j])){
system(paste(gdalwarp, ' ', tf, '.tif', ' -t_srs \"', CRS, '\" ', tf, '_ll.tif -dstnodata \"', NAflag, '\" -r near', sep=""), show.output.on.console = show.output.on.console)
}
else {
system(paste(gdalwarp, ' ', tf, '.tif', ' -t_srs \"', CRS, '\" ', tf, '_ll.tif -dstnodata \"', NAflag, '\" -r bilinear', sep=""), show.output.on.console = show.output.on.console)
}
## read images back to R:
if(j==1){
res <- readGDAL(paste(tf, "_ll.tif", sep=""), silent = TRUE)
names(res) <- names(obj)[j]
}
else{
res@data[names(obj)[j]] <- readGDAL(paste(tf, "_ll.tif", sep=""), silent = TRUE)$band1
}
## reformat to the original factors:
if(is.factor(obj@data[,j])){
res@data[,j] <- as.factor(res@data[,j])
try( levels(res@data[,j]) <- levels(obj@data[,j]) , silent = TRUE)
}
unlink(paste(tf, ".tif", sep=""))
unlink(paste(tf, "_ll.tif", sep=""))
}
}
else {
stop("Could not locate GDAL. See 'plotKML.env()' for more info.")
}
}
## TH: time consuming but preferred:
# res <- as(res, "SpatialPixelsDataFrame")
return(res)
}
reproject.RasterStack <- function(obj, CRS = get("ref_CRS", envir = plotKML.opts)) {
rs <- list(NULL)
for(j in 1:nlayers(obj)){
rs[[j]] <- reproject(obj[[j]], CRS = CRS)
}
return(stack(rs))
}
reproject.RasterBrick <- function(obj, CRS = get("ref_CRS", envir = plotKML.opts)) {
r <- stack(obj)
rs <- reproject.RasterStack(r, CRS = CRS)
return(brick(rs))
}
# connect all methods and classes:
setMethod("reproject", "SpatialPoints", reproject.SpatialPoints)
setMethod("reproject", "SpatialPolygons", reproject.SpatialPoints)
setMethod("reproject", "SpatialLines", reproject.SpatialPoints)
setMethod("reproject", "RasterLayer", reproject.RasterLayer)
setMethod("reproject", "SpatialGridDataFrame", reproject.SpatialGrid)
setMethod("reproject", "SpatialPixelsDataFrame", reproject.SpatialGrid)
setMethod("reproject", "RasterStack", reproject.RasterStack)
setMethod("reproject", "RasterBrick", reproject.RasterBrick)
# end of script;