如何使 R 的“光栅”包区分 GeoTIFF 中的正旋转矩阵和负旋转矩阵？
Posted
技术标签:
【中文标题】如何使 R 的“光栅”包区分 GeoTIFF 中的正旋转矩阵和负旋转矩阵？【英文标题】：How to make R's 'raster' package distinguish between positive and negative rotation matrices in GeoTIFFs? 【发布时间】：2017-06-21 03:56:53 【问题描述】：R 中的光栅包似乎无法区分 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
这是一个错误，还是我遗漏了什么？ raster 软件包非常强大和有用，我宁愿帮助添加更多功能，也不愿使用其他软件来正确处理这些（非常烦人的）旋转 tiff。谢谢！这里还有一个与旋转 tiff 相关的 R-sig-Geo mailing post。
【问题讨论】：
你问已经快一年了，但是当我查看raster:::.rasterFromGDAL 时，使用raster 软件包的当前版本，我确实看到了一个关于旋转文件的warning。据此判断，我不确定这是否是获得更好答案的正确地方。
【参考方案1】：
编辑
我想这里的大多数人都无法访问下面提供的修复程序，因此我已经很好地总结了这一点，以便人们可以检查和评论。
我已从CRAN 上的raster 软件包的当前版本（2.6-7）获取源代码：https://cran.r-project.org/web/packages/raster/index.html
并从那里创建了一个新的 Github 存储库。
之后，我提交了建议的旋转修复、一些相关测试和旋转的tiff以用于这些。最后我添加了一些onLoad 消息，明确表明这不是raster 包的正式版本。
您现在只需运行以下命令即可进行测试：
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)
对于提供的示例，我可以通过对raster:::.rasterFromGDAL 进行以下修改来修复它（参见 cmets addition 1 和 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)
我已经用R.tif 和 +45、-45、+315、+100 和 -100 的旋转对此进行了测试，这一切看起来都符合我的预期。
同时，鉴于代码中的warning，我预计旋转文件存在更深层次的潜在问题，所以我不能说这可能会让你走多远。
【讨论】：
以上是关于如何使 R 的“光栅”包区分 GeoTIFF 中的正旋转矩阵和负旋转矩阵？的主要内容，如果未能解决你的问题，请参考以下文章
使用 Polygon 裁剪光栅文件并使用相同的文件名写入输出