xyz traslacion transformación rotación rotacion robotica resueltos matriz matrices las identificar ejes ejercicios desplazamiento coordenadas cartesianas r raster image-rotation r-raster geotiff

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 (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.

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.