It appears that the raster package in R doesn't distinguish between positive and negative rotations of GeoTIFFs. I have a feeling that it is because R is ignoring the negative signs in the rotation matrix. I'm not quite savvy enough to dig down into the raster
source code to verify, but I did create a reproducible example that demonstrates the problem:
Read R
logo and save as GeoTIFF.
library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
proj4string(b) <- crs("+init=epsg:32616")
writeRaster(b, "R.tif")
Add rotation to tiff with Python
import sys
from osgeo import gdal
from osgeo import osr
import numpy as np
from math import *
def array2TIFF(inputArray,gdalData,datatype,angle,noData,outputTIFF):
# this script takes a numpy array and saves it to a geotiff
# given a gdal.Dataset object describing the spatial atributes of the data set
# the array datatype (as a gdal object) and the name of the output raster, and rotation angle in degrees
# get the file format driver, in this case the file will be saved as a GeoTIFF
driver = gdal.GetDriverByName("GTIFF")
#set the output raster properties
tiff = driver.Create(outputTIFF,gdalData.RasterXSize,gdalData.RasterYSize,inputArray.shape[0],datatype)
transform = []
originX = gdalData.GetGeoTransform()[0]
cellSizeX = gdalData.GetGeoTransform()[1]
originY = gdalData.GetGeoTransform()[3]
cellSizeY = gdalData.GetGeoTransform()[5]
rotation = np.radians(angle)
transform.append(originX)
transform.append(cos(rotation) * cellSizeX)
transform.append(sin(rotation) * cellSizeX)
transform.append(originY)
transform.append(-sin(rotation) * cellSizeY)
transform.append(cos(rotation) * cellSizeY)
transform = tuple(transform)
#set the geotransofrm values which include corner coordinates and cell size
#once again we can use the original geotransform data because nothing has been changed
tiff.SetGeoTransform(transform)
#next the Projection info is defined using the original data
tiff.SetProjection(gdalData.GetProjection())
#cycle through each band
for band in range(inputArray.shape[0]):
#the data is written to the first raster band in the image
tiff.GetRasterBand(band+1).WriteArray(inputArray[band])
#set no data value
tiff.GetRasterBand(band+1).SetNoDataValue(0)
#the file is written to the disk once the driver variables are deleted
del tiff, driver
inputTif = gdal.Open("R.tif")
inputArray = inputTif.ReadAsArray()
array2TIFF(inputArray,inputTif, gdal.GDT_Float64, -45, 0, "R_neg45.tif")
array2TIFF(inputArray,inputTif, gdal.GDT_Float64, 45, 0, "R_pos45.tif")
Read in the rotated tiffs in R
.
c <- brick("R_neg45.tif")
plotRGB(c,1,2,3)
d <- brick("R_pos45.tif")
plotRGB(d,1,2,3)
> c
class : RasterBrick
rotated : TRUE
dimensions : 77, 101, 7777, 3 (nrow, ncol, ncell, nlayers)
resolution : 0.7071068, 0.7071068 (x, y)
extent : 0, 125.865, 22.55278, 148.4178 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /Users/erker/g/projects/uft/code/R_neg45.tif
names : R_neg45.1, R_neg45.2, R_neg45.3
> d
class : RasterBrick
rotated : TRUE
dimensions : 77, 101, 7777, 3 (nrow, ncol, ncell, nlayers)
resolution : 0.7071068, 0.7071068 (x, y)
extent : 0, 125.865, 22.55278, 148.4178 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /Users/erker/g/projects/uft/code/R_pos45.tif
names : R_pos45.1, R_pos45.2, R_pos45.3
The plots are the same and note the equivalent extents. However, gdalinfo
tells a different story
$ gdalinfo R_neg45.tif
Driver: GTiff/GeoTIFF
Files: R_neg45.tif
Size is 101, 77
Coordinate System is:
PROJCS["WGS 84 / UTM zone 16N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-87],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32616"]]
GeoTransform =
0, 0.7071067811865476, -0.7071067811865475
77, -0.7071067811865475, -0.7071067811865476
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 0.0000000, 77.0000000) ( 91d29'19.48"W, 0d 0' 2.50"N)
Lower Left ( -54.4472222, 22.5527778) ( 91d29'21.23"W, 0d 0' 0.73"N)
Upper Right ( 71.4177849, 5.5822151) ( 91d29'17.17"W, 0d 0' 0.18"N)
Lower Right ( 16.9705627, -48.8650071) ( 91d29'18.93"W, 0d 0' 1.59"S)
Center ( 8.4852814, 14.0674965) ( 91d29'19.20"W, 0d 0' 0.46"N)
Band 1 Block=101x3 Type=Float64, ColorInterp=Gray
NoData Value=0
Band 2 Block=101x3 Type=Float64, ColorInterp=Undefined
NoData Value=0
Band 3 Block=101x3 Type=Float64, ColorInterp=Undefined
NoData Value=0
$ gdalinfo R_pos45.tif
Driver: GTiff/GeoTIFF
Files: R_pos45.tif
Size is 101, 77
Coordinate System is:
PROJCS["WGS 84 / UTM zone 16N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-87],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32616"]]
GeoTransform =
0, 0.7071067811865476, 0.7071067811865475
77, 0.7071067811865475, -0.7071067811865476
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 0.0000000, 77.0000000) ( 91d29'19.48"W, 0d 0' 2.50"N)
Lower Left ( 54.4472222, 22.5527778) ( 91d29'17.72"W, 0d 0' 0.73"N)
Upper Right ( 71.418, 148.418) ( 91d29'17.17"W, 0d 0' 4.82"N)
Lower Right ( 125.865, 93.971) ( 91d29'15.42"W, 0d 0' 3.05"N)
Center ( 62.9325035, 85.4852814) ( 91d29'17.45"W, 0d 0' 2.78"N)
Band 1 Block=101x3 Type=Float64, ColorInterp=Gray
NoData Value=0
Band 2 Block=101x3 Type=Float64, ColorInterp=Undefined
NoData Value=0
Band 3 Block=101x3 Type=Float64, ColorInterp=Undefined
NoData Value=0
Is this a bug, or am I missing something? The raster
package is incredibly powerful and useful, and I'd rather help add more functionality than have to use other software to properly handle these (very annoyingly) rotated tiffs. Thanks! Also here's a R-sig-Geo mailing post related to rotated tiffs.
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With