在 R 中，似乎 raster 包不区分 GeoTIFF 的正旋转和负旋转。我有一种感觉，那就是 R 忽略了旋转矩阵中的负号。我不太懂得深入挖掘 raster 源代码来验证，但我确实创建了一个可重现的示例来演示这个问题:
读取 R 徽标并保存为 GeoTIFF。
library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
proj4string(b) <- crs("+init=epsg:32616")
writeRaster(b, "R.tif")
使用Python为tiff文件添加旋转功能
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")
使用R读取旋转后的tiff文件。
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
图形相同，注意等价范围。然而，gdalinfo 告诉了一个不同的故事。
$ 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
这是一个bug吗，还是我漏掉了什么？raster包非常强大和有用，我宁愿帮助添加更多功能，而不是使用其他软件来正确处理这些（非常烦人的）旋转tiff文件。谢谢！另外，这里有一篇与旋转tiff相关的R-sig-Geo邮件帖子。
raster包时，当我查看raster:::.rasterFromGDAL时确实会出现有关旋转文件的warning。根据这个判断，我不确定这是否是获取更好答案的正确位置。 - RolandASc