# Author: Robert J. Hijmans
# Date : January 2009
# Version 1.0
# Licence GPL v3
.useproj6 <- function() {
FALSE
}
.rawTransform <- function(projfrom, projto, xy, wkt="") {
v <- terra::vect()
out <- terra::project(xy, projfrom, projto)
matrix(out, ncol=2)
# xy <- terra::vect(xy, crs=projfrom)
# xy <- terra::project(xy, projto)
# terra::crds(xy)
}
projectExtent <- function(object, crs) {
use_proj6 <- .useproj6()
object <- raster(object)
dm <- oldm <- dim(object)
# simple way to avoid a bug with a single column/row reported by
# Jon Olav Skoien
dm[1] <- max(10, dm[1])
dm[2] <- max(10, dm[2])
dim(object) <- dm
pfrom <- .getCRS(object)
pto <- .getCRS(crs)
if (use_proj6) {
projfrom <- wkt(pfrom)
projto <- wkt(pto)
if (is.null(projfrom) || is.null(projto)) {
use_proj6 = FALSE
projfrom <- pfrom
projto <- pto
}
} else {
projfrom <- proj4string(pfrom)
projto <- proj4string(pto)
}
# rs <- res(object)
# xmn <- object@extent@xmin - 0.5 * rs[1]
# xmx <- object@extent@xmax + 0.5 * rs[1]
# ymn <- object@extent@ymin - 0.5 * rs[2]
# ymx <- object@extent@ymax + 0.5 * rs[2]
# xha <- (xmn + xmx) / 2
# yha <- (ymn + ymx) / 2
# xy <- matrix(c(xmn, ymx, xha, ymx, xmx, ymx, xmn, yha, xha, yha, xmx, yha, xmn, ymn, xha, ymn, xmx, ymn), ncol=2, byrow=T)
rows <- unique(c(seq(1,nrow(object), by=max(1, round(nrow(object)/50))), nrow(object)))
cols <- unique(c(seq(1,ncol(object), by=max(1, round(ncol(object)/50))), ncol(object)))
xy1 <- xyFromCell(object, cellFromRowCol(object, rows, 1))
xy1[,1] <- xy1[,1] - 0.5 * xres(object)
xy1[1,2] <- xy1[1,2] + 0.5 * yres(object)
xy1[nrow(xy1),2] <- xy1[nrow(xy1),2] + 0.5 * yres(object)
xy2 <- xyFromCell(object, cellFromRowCol(object, rows, ncol(object)))
xy2[,1] <- xy2[,1] + 0.5 * xres(object)
xy2[1,2] <- xy2[1,2] + 0.5 * yres(object)
xy2[nrow(xy2),2] <- xy2[nrow(xy2),2] + 0.5 * yres(object)
xy3 <- xyFromCell(object, cellFromRowCol(object, 1, cols))
xy3[,2] <- xy3[,2] + 0.5 * yres(object)
xy3[1,1] <- xy3[1,1] - 0.5 * xres(object)
xy3[ncol(xy3),1] <- xy3[ncol(xy3),1] + 0.5 * xres(object)
xy4 <- xyFromCell(object, cellFromRowCol(object, nrow(object), cols))
xy4[,2] <- xy4[,2] - 0.5 * yres(object)
xy4[1,1] <- xy4[1,1] - 0.5 * xres(object)
xy4[ncol(xy4),1] <- xy4[ncol(xy4),1] + 0.5 * xres(object)
# added for circumpolar data:
if (nrow(object) > 75 & ncol(object) > 75) {
xy5 <- sampleRegular(object, 500, xy=TRUE)
# rows <- c(seq(min(nrow(object), 25), nrow(object), by=50))
# cols <- c(seq(min(ncol(object), 25), ncol(object), by=50))
# xy5 <- xyFromCell(object, cellFromRowColCombine(object, rows, cols))
xy <- rbind(xy1, xy2, xy3, xy4, xy5)
} else {
xy <- rbind(xy1, xy2, xy3, xy4)
}
xy <- .rawTransform( projfrom, projto, xy)
xy <- subset(xy, !(is.infinite(xy[,1]) | is.infinite(xy[,2]) | is.na(xy[,2]) ))
x <- xy[,1]
y <- xy[,2]
if (length(y) == 0 | length(y) ==0) { stop("cannot do this transformation") }
minx <- min(x)
maxx <- max(x)
if (maxx == minx) {
maxx <- maxx + 0.5
minx <- minx - 0.5
}
miny <- min(y)
maxy <- max(y)
if (maxy == miny) {
maxy <- maxy + 0.5
miny <- miny - 0.5
}
obj <- raster(extent(minx, maxx, miny, maxy), nrows=oldm[1], ncols=oldm[2], crs=crs)
return(obj)
}
.computeRes <- function(obj, crs) {
x <- xmin(obj) + 0.5 * (xmax(obj) - xmin(obj))
y <- ymin(obj) + 0.5 * (ymax(obj) - ymin(obj))
res <- res(obj)
x1 <- x - 0.5 * res[1]
x2 <- x + 0.5 * res[1]
y1 <- y - 0.5 * res[2]
y2 <- y + 0.5 * res[2]
xy <- cbind(c(x1, x2, x, x), c(y, y, y1, y2))
#fromcrs <- .getCRS(obj)
fromsrs <- .getSRS(obj)
pXY <- .rawTransform(fromsrs, crs, xy)
# out <- c((pXY[2,1] - pXY[1,1]), (pXY[4,2] - pXY[3,2]))
outex <- extent(pXY)
out <- c(xmax(outex) - xmin(outex), ymax(outex) - ymin(outex))
if (any(is.na(out))) {
if (isLonLat(obj)) {
out <- pointDistance(cbind(x1, y1), cbind(x2, y2), lonlat=TRUE)
out <- c(out, out)
} else {
out <- res
}
}
# abs should not be necessary, but who knows what a projection might do?
abs( signif(out, digits=3) )
}
.getAlignedRaster <- function(x,y) {
x <- raster(x)
y <- raster(y)
p <- projectRaster(x, crs=.getCRS(y))
m <- merge(extent(y), extent(p))
rx <- extend(y, m)
crop(rx, p)
}
projectRaster <- function(from, to, res, crs, method="bilinear", alignOnly=FALSE, over=FALSE, filename="", ...) {
#projfrom <- .getCRS(from)
projfrom <- .getSRS(from)
if (is.na(projfrom)) {
stop("input projection is NA")
}
lonlat <- isLonLat(projfrom)
if (missing(to)) {
if (missing(crs)) {
stop("both 'to' and 'crs' arguments are missing.")
}
#projto <- .getCRS(crs)
projto <- .getSRS(crs)
#compareCRS(projfrom, projto)
to <- projectExtent(from, projto)
# to@crs <- projto
to@srs <- projto
if (missing(res)) {
res <- .computeRes(from, projto)
}
res(to) <- res
# add some cells to capture curvature
e <- extent(to)
add <- min(5, min(dim(to)[1:2])/10) * max(res)
e@ymin <- e@ymin - add
e@ymax <- e@ymax + add
e@xmin <- e@xmin - add
e@xmax <- e@xmax + add
if (!is.character(projto)) projto <- projto@projargs
if (substr(projto, 1, 13) == "+proj=longlat") {
e@xmin <- max(-180, e@xmin)
e@xmax <- min(180, e@xmax)
e@ymin <- max(-90, e@ymin)
e@ymax <- min(90, e@ymax)
}
to <- extend(to, e)
} else {
# projto <-.getCRS(to)
projto <- .getSRS(to)
if (is.na(projto) || (projto == "")) {
stop("output projection is empty")
}
e <- extent( projectExtent(from, projto) )
add <- min(10, min(dim(to)[1:2])/10) * max(raster::res(to))
e@ymin <- e@ymin - add
e@ymax <- e@ymax + add
e@xmin <- e@xmin - add
e@xmax <- e@xmax + add
if (isLonLat(projto)) {
e@xmin <- max(-180, e@xmin)
e@xmax <- min(180, e@xmax)
e@ymin <- max(-90, e@ymin)
e@ymax <- min(90, e@ymax)
}
}
methods::validObject(to)
# methods::validObject(.getCRS((to)))
#if (identical(projfrom, projto)) {
# warning('projections of "from" and "to" are the same')
#}
# if ((!use_proj6) & lonlat & over) {
# projto_int <- paste(projto, "+over")
# } else {
projto_int <- projto
# }
if (alignOnly) {
to <- .getAlignedRaster(from, to)
return (to)
}
# pbb <- projectExtent(to,.getCRS(from))
# bb <- intersect(extent(pbb), extent(from))
# methods::validObject(bb)
if (!method %in% c('bilinear', 'ngb')) {
stop('invalid method')
}
nl <- nlayers(from)
if ( nl == 1) {
to <- raster(to)
if (method=="ngb") {
colortable(to) <- colortable(from)
}
} else {
to <- brick(to, values=FALSE, nl=nl)
}
if (method=='ngb') {
method <- 'simple' # for extract (.xyValues)
}
names(to) <- names(from)
if ( ! hasValues(from) ) {
#warning("'from' has no cell values")
return(to)
}
if (canProcessInMemory(to, n=nl*4)) {
inMemory <- TRUE
} else {
inMemory <- FALSE
}
# this seems to need smaller chunks
#cz <- max(5, 0.1 * .chunk() / nlayers(to))
if (inMemory) {
xy <- coordinates(to)
xy <- subset(xy, xy[,1] > e@xmin & xy[,1] < e@xmax)
cells <- cellFromXY(to, xy)
xy <- .rawTransform( projto_int, projfrom, xy)
to[cells] <- .xyValues(from, xy, method=method)
if (filename != '') {
to <- writeRaster(to, filename, ...)
}
return(to)
} else {
tr <- blockSize(to, n=nlayers(to)*4)
pb <- pbCreate(tr$n, label='projectRaster', ...)
to <- writeStart(to, filename=filename, ...)
for (i in 1:tr$n) {
cells <- cellFromRowCol(to, tr$row[i], 1):cellFromRowCol(to, tr$row[i]+tr$nrows[i]-1, ncol(to))
xy <- xyFromCell(to, cells )
xy <- subset(xy, xy[,1] > e@xmin & xy[,1] < e@xmax)
if (nrow(xy) > 0) {
ci <- match(cellFromXY(to, xy), cells)
xy <- .rawTransform( projto_int, projfrom, xy)
v <- matrix(nrow=length(cells), ncol=nl)
v[ci, ] <- .xyValues(from, xy, method=method)
to <- writeValues(to, v, tr$row[i])
}
pbStep(pb)
}
pbClose(pb)
to <- writeStop(to)
return(to)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.