आर के 'रास्टर' पैकेज को GeoTIFFs में सकारात्मक और नकारात्मक रोटेशन मैट्रिस के बीच अंतर कैसे बनाया जाए?

Question

ऐसा प्रतीत होता है कि आर में रास्टर पैकेज GeoTIFFs के सकारात्मक और नकारात्मक घूर्णन के बीच अंतर नहीं करता है। मुझे एहसास है कि ऐसा इसलिए है क्योंकि आर रोटेशन मैट्रिक्स में नकारात्मक संकेतों को अनदेखा कर रहा है। मैं काफी पर्याप्त समझ रखने वाले को सत्यापित करने के raster स्रोत कोड में नीचे खुदाई करने के लिए नहीं कर रहा हूँ, लेकिन मैं एक प्रतिलिपि प्रस्तुत करने योग्य उदाहरण है कि समस्या को दर्शाता है बनाने के लिए क्या किया:

पढ़ें R लोगो और GeoTiff रूप में सहेजें।

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

writeRaster(b, "R.tif") 

जोड़ें टिफ को रोटेशन अजगर

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 में घुमाया झगड़े में।

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 पैकेज अविश्वसनीय रूप से शक्तिशाली और उपयोगी है, और मैं इन सॉफ़्टवेयर का उपयोग करने के बजाय अन्य सॉफ़्टवेयर का उपयोग करने की अपेक्षा अधिक कार्यक्षमता जोड़ने में मदद करता हूं (बहुत परेशान) घुमावदार टिफ्स। धन्यवाद! यहां घूर्णन वाले टिफ से संबंधित एक आर-सिग-जिओ mailing post भी है।

Answer 1

संपादित

मुझे लगता है कि नीचे प्रस्तुत ठीक यहाँ ज्यादातर लोगों के लिए सुलभ नहीं था, इसलिए मैं अच्छी तरह से इस की समाप्ति है, लोगों को उम्मीद है कि जाँच करें और टिप्पणी कर सकते हैं जैसे कि।

मैं CRAN पर raster पैकेज की मौजूदा रिलीज (2.6-7) से स्रोतों ले लिया है:
https://cran.r-project.org/web/packages/raster/index.html
और वहां से एक नया Github भंडार बनाया।

बाद में, मैं प्रतिबद्ध है प्रस्तावित रोटेशन-ठीक, जुड़े परीक्षण और के एक मुट्ठी भर झगड़े घुमाया उन में उपयोग करने के लिए। आखिरकार मैंने कुछ onLoad संदेशों को स्पष्ट रूप से इंगित करने के लिए जोड़ा कि यह raster पैकेज का आधिकारिक संस्करण नहीं है।

अब आप बस निम्न कमांड चलाकर परीक्षण कर सकते हैं: मैं raster:::.rasterFromGDAL के लिए निम्न संशोधनों के साथ इसे ठीक कर सकता है, बशर्ते उदाहरण के लिए

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) 

---

(टिप्पणियों को देखने के इसके अलावा 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) 

मैंके साथ इस परीक्षण किया हैऔर +45, -45, +315, +100 और -100 के घूर्णन, जो सभी मेरी अपेक्षा की तरह दिखते हैं।

उसी समय, कोड में warning दिए गए, मुझे उम्मीद है कि घूर्णन वाली फ़ाइलों के साथ गहरे संभावित समस्याएं हैं, इसलिए मैं यह नहीं कह सकता कि यह आपको कितना दूर ले सकता है।