Jak sprawić, by paczka R "rastrowa" odróżniała macierze rotacji dodatniej i ujemnej w GeoTIFF?

Question

Wygląda na to, że pakiet rastrowy w R nie rozróżnia dodatnich i ujemnych obrotów GeoTIFF. Mam wrażenie, że to dlatego, że R ignoruje negatywne znaki w macierzy obrotu. Nie jestem na tyle sprytny, aby zagłębić się w kod źródłowy raster w celu zweryfikowania, ale utworzyłem powtarzalny przykład, który demonstruje problem:

Przeczytaj R logo i zapisz jako GeoTIFF.

library(raster) 
b <- brick(system.file("external/rlogo.grd", package="raster")) 
proj4string(b) <- crs("+init=epsg:32616") 

writeRaster(b, "R.tif") 

dodawania rotacji do tiff z Pythonem

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

przeczytać w obróconych TIFF w 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 

Wykresy są takie same i należy zwrócić uwagę na równoważny zakres. Jednak gdalinfo opowiada inną historię

$ 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 

Jest to błąd, albo ja czegoś brakuje? Pakiet raster jest niesamowicie potężny i użyteczny, a ja wolałbym pomóc w dodaniu większej funkcjonalności niż użycie innego oprogramowania do prawidłowego radzenia sobie z tymi (bardzo irytująco) obróconymi tiffami. Dzięki! Również tutaj jest R-sig-Geo mailing post związane z obróconych tiffs.

Answer 1

EDIT

Przypuszczam, że przedstawiony poniżej poprawka nie była dostępna dla większości ludzi tutaj, więc ja to ładnie opakowane tak, że ludzie mogą z nadzieją sprawdzić i komentarz.

Brałem źródła z bieżącego wydania (2.6-7) opakowania raster na CRAN:
https://cran.r-project.org/web/packages/raster/index.html
i stworzył nowe repozytorium GitHub stamtąd.

Potem mam popełnione proponowany obrót-fix, garść związanych testów i obracany TIFF do wykorzystania w nich. Na koniec dodałem kilka komunikatów onLoad, aby wyraźnie zaznaczyć, że nie jest to oficjalna wersja pakietu raster.

Teraz można przetestować po prostu uruchamiając następujące polecenia:

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) 

---

Dla przykładu przewidzianego mogłem go naprawić z następującymi modyfikacjami do raster:::.rasterFromGDAL (patrz uwagi ponadto 1 i dodatek 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) 

Testowałem to z R.tif i obroty o wartości +45, -45, +315, +100 i -100, które wyglądają tak, jakbym oczekiwał.

W tym samym czasie, biorąc pod uwagę warning w kodzie, spodziewałbym się, że istnieją głębsze potencjalne problemy z obracanymi plikami, więc nie mogę powiedzieć, jak daleko to zajmie.