Wie unterscheidet das R-Paket von R zwischen positiven und negativen Rotationsmatrizen in GeoTIFFs?

Question

Es scheint, dass das Rasterpaket in R nicht zwischen positiven und negativen Rotationen von GeoTIFFs unterscheidet. Ich habe das Gefühl, dass dies daran liegt, dass R die negativen Vorzeichen in der Rotationsmatrix ignoriert. Ich bin nicht ganz clever genug, um in den raster Quellcode zu graben, um zu überprüfen, aber ich habe ein reproduzierbares Beispiel erstellt, das das Problem demonstriert:

R lesen und als GeoTIFF speichern.

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

writeRaster(b, "R.tif") 

Add Rotation zu tiff mit 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") 

liest in der gedrehten 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 

Die Plots sind die gleichen und beachten Sie die entsprechenden Extents. Allerdings sagt gdalinfo eine andere Geschichte

$ 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 

Ist das ein Bug, oder bin ich etwas fehlt? Das raster Paket ist unglaublich leistungsfähig und nützlich, und ich würde eher helfen, mehr Funktionalität hinzuzufügen, als andere Software zu verwenden, um diese (sehr ärgerlichen) gedrehten Tiffs richtig zu handhaben. Vielen Dank! Auch hier ist ein R-Sig-Geo mailing post bezogen auf rotierte Tiffs.

Answer 1

EDIT

Ich nehme an, dass die im Folgenden dargestellt fix hier, um die die meisten Menschen nicht zugänglich war, deshalb habe ich das oben schön verpackt, so dass die Menschen überprüfen und kommentieren hoffentlich können.
https://cran.r-project.org/web/packages/raster/index.html
und erstellt eine neue Github-Repository von dort:

Ich habe die Quellen aus der aktuellen Version (2.6-7) des raster Pakets auf CRAN genommen.

Danach habe ich begangen die vorgeschlagene Rotations-fix, eine Handvoll zugehörigen Prüfungen und tiffs in solche zu verwenden, gedreht. Schließlich fügte ich einige onLoad Nachrichten hinzu, um deutlich anzuzeigen, dass dies keine offizielle Version des raster Pakets ist.

Sie können nun testen, indem Sie einfach die folgenden Befehle ausführen:

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) 

---

Für das Beispiel vorgesehen ich es mit den folgenden Modifikationen raster:::.rasterFromGDAL beheben konnte (siehe Kommentare zusätzlich 1 und zusätzlich 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) 

ich habe dies mit demgetestetund Rotationen von +45, -45, +315, +100 und -100, die alle so aussehen, wie ich es erwarten würde.

Zur gleichen Zeit, angesichts der warning im Code, würde ich erwarten, dass es tiefere potentielle Probleme mit gedrehten Dateien gibt, so dass ich nicht sagen kann, wie weit dies Sie bringen könnte.