How to make R's 'raster' package distinguish between positive and negative rotation matrices in GeoTIFFs?
EDIT
I suppose that the presented fix below was not accessible to most people here, therefore I have wrapped this up nicely, such that people can hopefully check and comment.
I have taken the sources from the current release (2.6-7
) of the raster
package on CRAN
:
https://cran.r-project.org/web/packages/raster/index.html
and created a new Github repository from there.
Afterwards, I have committed the proposed rotation-fix, a handful of associated tests and rotated tiffs to use in those. Finally I added some onLoad
messages to indicate clearly that this is not an official version of the raster
package.
You can now test by simply running the following commands:
devtools::install_github("miraisolutions/raster")
library(raster)
## modified raster 2.6-7 (2018-02-23)
## you are using an unofficial, non-CRAN version of the raster package
R_Tif <- system.file("external", "R.tif", package = "raster", mustWork = TRUE)
R_Tif_pos45 <- system.file("external", "R_pos45.tif", package = "raster", mustWork = TRUE)
R_Tif_neg45 <- system.file("external", "R_neg45.tif", package = "raster", mustWork = TRUE)
R_Tif_pos100 <- system.file("external", "R_pos100.tif", package = "raster", mustWork = TRUE)
R_Tif_neg100 <- system.file("external", "R_neg100.tif", package = "raster", mustWork = TRUE)
R_Tif_pos315 <- system.file("external", "R_pos315.tif", package = "raster", mustWork = TRUE)
RTif <- brick(R_Tif)
plotRGB(RTif, 1, 2, 3)
pos45Tif <- suppressWarnings(brick(R_Tif_pos45))
plotRGB(pos45Tif, 1, 2, 3)
neg45Tif <- suppressWarnings(brick(R_Tif_neg45))
plotRGB(neg45Tif,1,2,3)
pos100Tif <- suppressWarnings(brick(R_Tif_pos100))
plotRGB(pos100Tif, 1, 2, 3)
neg100Tif <- suppressWarnings(brick(R_Tif_neg100))
plotRGB(neg100Tif, 1, 2, 3)
pos315Tif <- suppressWarnings(brick(R_Tif_pos315))
plotRGB(pos315Tif,1,2,3)
For the example provided I could fix it with the following modifications to raster:::.rasterFromGDAL
(see comments addition 1 and addition 2):
# ... (unmodified initial part of function)
# condition for rotation case
if (gdalinfo["oblique.x"] != 0 | gdalinfo["oblique.y"] != 0) {
rotated <- TRUE
res1 <- attributes(rgdal::readGDAL(filename))$bbox # addition 1
if (warn) {
warning("\n\n This file has a rotation\n Support for such files is limited and results of data processing might be wrong.\n Proceed with caution & consider using the \"rectify\" function\n")
}
rotMat <- matrix(gdalinfo[c("res.x", "oblique.x", "oblique.y", "res.y")], 2)
# addition 2 below
if (all(res1[, "min"] < 0)) {
rotMat[2] <- rotMat[2] * -1
rotMat[3] <- rotMat[3] * -1
}
# ... (original code continues)
I have tested this with the R.tif
and rotations of +45, -45, +315, +100 and -100, which all look like what I would expect.
At the same time, given the warning
in the code, I would expect that there are deeper potential issues with rotated files, so I cannot say how far this might take you.