50
R의 래스터 패키지는 GeoTIFF의 양수 회전과 음수 회전을 구분하지 않습니다. 나는 R이 회전 행렬에서 음의 부호를 무시하고 있기 때문에 그것이라고 느낀다. 나는 확인하기 위해 raster 소스 코드를 파고 들기에 충분하지는 않지만 문제를 나타내는 재현 가능한 예제를 만들었습니다.R 's '래스터'패키지로 GeoTIFF의 양수 및 음수 회전 행렬을 구별하는 방법은 무엇입니까?
GeoTIFF로 저장하고 GeoTIFF로 저장하십시오. R의 회전 tiffs 파이썬
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")
읽기와 티파니에
library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
proj4string(b) <- crs("+init=epsg:32616")
writeRaster(b, "R.tif")
추가 회전.
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를 적절하게 처리하는 것보다 더 많은 기능을 추가하는 것을 돕고 싶습니다. 감사! 회전 된 tiff와 관련된 R-sig-Geo mailing post도 있습니다.
당신이 물어 본 지 거의 1 년이 지났지 만,'raster :::. rasterFromGDAL'을 보면'raster' 패키지의 현재 버전에서 회전 된 파일에 대한'경고'를 볼 수 있습니다. 이것으로 판단 할 때 이것이 더 나은 대답을 얻을 수있는 적절한 곳인지 확신 할 수 없습니다. – RolandASc