Rの「ラスタ」パッケージを作成すると、GeoTIFFの正と負の回転行列が区別されますか？

Question

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") 

Rに回転TIFFファイルには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") 

読むと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を適切に処理するために他のソフトウェアを使用するよりも多くの機能を追加するのに役立ちます。ありがとう！回転Tiffsに関連するR-sig-Geo mailing postもここにあります。

Answer 1

EDIT

私は、以下のとおり修正はここでほとんどの人がアクセス可能ではなかったと仮定ので、私がうまくこの包まれている、人々はうまくいけばチェックして、コメントできるようになっています。

私はCRANにrasterパッケージの現在のリリース（2.6-7）からソースを取ってきた：
https://cran.r-project.org/web/packages/raster/index.html
をし、そこから新しいGithubのリポジトリを作成しました。その後

、Iは関連試験との一握りは、TIFFを回転もので使用する、提案を回転固定コミットしています。最後に、rasterパッケージの公式バージョンではないことを明確に示すために、onLoadメッセージを追加しました。

あなたは今、単純に次のコマンドを実行してテストすることができます

：私はraster:::.rasterFromGDALに以下のように変更して、それを修正することができてたとえば

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) 

---

（コメントさらに1と加え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) 

私はこれをでテストしましたと+45、-45、+155、+100、-100の回転です。これらはすべて私が期待するように見えます。

同時にコード内にwarningがあると、回転したファイルの潜在的な問題が深刻になることが予想されるため、これでどれくらいの時間がかかるかはわかりません。