r - ضرب - طريقة حل المحددات بالالة الحاسبة



كيفية جعل حزمة R 'raster' تميز بين مصفوفات الدوران الموجب والسالب في GeoTIFFs؟ (1)

يبدو أن الحزمة النقطية في R لا تميز بين التدوير الإيجابي والسلبي في GeoTIFF. لدي شعور بأن السبب هو أن R يتجاهل العلامات السلبية في مصفوفة الدوران. لست دقيقًا بما فيه الكفاية للدخول إلى شفرة المصدر 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")

قراءة في tiffs استدارة في 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 قوية ومفيدة بشكل لا يصدق ، وأنا أفضل المساعدة في إضافة المزيد من الوظائف بدلاً من استخدام برامج أخرى للتعامل مع هذه المشاجرات المدورة (المزعجة جدًا) بشكل صحيح. شكر! أيضا هنا R-sig-Geo البريدية المشاركة المتعلقة tiffs استدارة.


تعديل

أفترض أن الإصلاح المقدم أدناه لم يكن متاحًا لمعظم الأشخاص هنا ، لذا فقد انتهيت من هذا الأمر بشكل جيد ، بحيث يمكن للأشخاص التحقق والتعليق.

لقد أخذت المصادر من الإصدار الحالي ( 2.6-7 ) من الحزمة raster على CRAN :
https://cran.r-project.org/web/packages/raster/index.html
وإنشاء مستودع جيثب جديد من هناك.

بعد ذلك ، لقد ارتكبت إصلاح التناوب المقترح ، وحفنة من الاختبارات المرتبطة بها والمشاجرات الدوارة لاستخدامها في تلك. وأخيرًا أضفت بعض رسائل 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 الوارد في الكود ، أتوقع أن تكون هناك مشكلات محتملة أعمق مع الملفات المدورة ، لذلك لا يمكنني تحديد المدى الذي قد يستغرقه هذا الأمر.