1. 首頁
  2. 問答
  3. 如何使R的'光柵'包區分GeoTIFF中的正旋矩陣和負旋轉矩陣？

如何使R的'光柵'包區分GeoTIFF中的正旋矩陣和負旋轉矩陣？

2017-02-03 62 views
50

似乎R中的柵格包並沒有區分GeoTIFF的正向旋轉和負向旋轉。我有一種感覺，那是因爲R忽略了旋轉矩陣中的負號。我沒有足夠的精力深入到源代碼raster來驗證，但我確實創建了一個可重現的示例來演示問題：如何使R的'光柵'包區分GeoTIFF中的正旋矩陣和負旋轉矩陣？

閱讀R徽標並保存爲GeoTIFF。

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

writeRaster(b, "R.tif")

添加旋轉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")

閱讀在R旋轉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。謝謝！這裏還有一個R-sig-Geo mailing post與旋轉tiffs相關。

來源

2017-02-03 Tedward

+0

自從你問了差不多一年了，但是對於當前版本的'raster'包，當我看'raster :::. rasterFromGDAL'時，我會看到關於旋轉文件的'warning'。從這個角度來看，我不確定這是否能找到更好的答案。 – RolandASc

回答

0

編輯

我想，下面介紹的修復是不是這裏的大多數人訪問，因此，我裹着這件事很好，這樣人們可以希望檢查和評論。

我承擔了CRAN從raster包的當前版本（2.6-7）來源：
https://cran.r-project.org/web/packages/raster/index.html
，並從那裏創建了一個新的Github上庫。

之後，我犯了提出旋轉固定的相關測試和一把旋轉TIFF格式在那些使用。最後我添加了一些onLoad消息來清楚地表明這不是raster包的正式版本。

現在，您可以通過簡單地運行以下命令進行測試：

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)

對於提供，我可以用下面的修改解決它raster:::.rasterFromGDAL的例子（見註釋另外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)

我已經測試了這個R.tif和旋轉+45，-45，+ 315，+ 100和-100，這些都看起來像我期望的。

與此同時，由於代碼中的warning，我預計旋轉文件存在更深層次的潛在問題，所以我不能說這會對您造成多大的影響。

來源

2018-01-31 11:22:13 RolandASc

相關問題

 相關問題