traslacion - ¿Cómo hacer que el paquete ''raster'' de R distinga entre matrices de rotación positivas y negativas en GeoTIFFs?

rotacion de coordenadas cartesianas

Question

Parece que el paquete ráster en R no distingue entre rotaciones positivas y negativas de GeoTIFF. Tengo la sensación de que es porque R está ignorando los signos negativos en la matriz de rotación. No soy lo suficientemente inteligente como para profundizar en el código fuente raster para verificar, pero creé un ejemplo reproducible que demuestra el problema:

Lea el logotipo R y guárdelo como GeoTIFF.

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

Agregue rotación a tiff con 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")

Lea en los tiffs rotados en 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 

Las parcelas son iguales y observan las extensiones equivalentes. Sin embargo, gdalinfo cuenta una historia diferente

$ 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

¿Es esto un error o me falta algo? El paquete raster es increíblemente poderoso y útil, y prefiero ayudar a agregar más funcionalidad que tener que usar otro software para manejar adecuadamente estos (muy molestos) tiffs rotados. ¡Gracias! También aquí hay una publicación de correo R-sig-Geo relacionada con tiffs rotados.

Answer 1

EDITAR

Supongo que la solución presentada a continuación no era accesible para la mayoría de las personas aquí, por lo tanto, he concluido esto muy bien, de modo que la gente pueda verificar y comentar.

Tomé las fuentes de la versión actual ( 2.6-7 ) del paquete raster en CRAN :
https://cran.r-project.org/web/packages/raster/index.html
y creó un nuevo repositorio de Github a partir de ahí.

Después, he comprometido la corrección de rotación propuesta , un puñado de pruebas asociadas y tiffs rotados para usar en esos. Finalmente agregué algunos mensajes onLoad para indicar claramente que esta no es una versión oficial del paquete raster .

Ahora puede probar simplemente ejecutando los siguientes comandos:

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)

Para el ejemplo proporcionado, podría solucionarlo con las siguientes modificaciones a raster:::.rasterFromGDAL (ver comentarios adición 1 y adición 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)

He probado esto con el R.tif y rotaciones de +45, -45, +315, +100 y -100, que se parecen a lo que esperaría.

Al mismo tiempo, dada la warning en el código, esperaría que haya problemas potenciales más profundos con los archivos rotados, por lo que no puedo decir qué tan lejos podría llevarlo.