[GIS] Opening rotated raster in R

landsatrraster

I have part of a LANDSAT image, this one: http://ge.tt/80uEl0n2

Using R, I would like to open this image, so I did:

B4<-raster(LE07_L1TP_219076_20021029_20170127_01_T1_B4.tif)

But it gave me this warning:

Warning message:
In .rasterFromGDAL(x, band = band, objecttype, …) :

This file has a rotation
Support such files is limited and results of data processing might be wrong.
Proceed with caution & consider using the "rectify" function

Then, using the rectify() function it gave me:

Error in if (value[1] != nrow(x) | value[2] != ncol(x)) { :
missing value where TRUE/FALSE needed
In addition: Warning message:
In dim<-(*tmp*, value = c(nr, nc)) :
NAs introduced by coercion to integer range

How can I solve this to work with my tif image? At least, when I call B4, it gives me:

class : RasterLayer

rotated : TRUE

dimensions : 505, 872, 440360 (nrow, ncol, ncell)

resolution : 4.662376e-11, 1.623406e-10 (x, y)

extent : 259785, 285945, 7476655, 7491805 (xmin, xmax, ymin, ymax)

coord. ref. : +proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

data source : D:\L7Referencia\LE072190762002102901T1-
SC20171002122402\LE07_L1TP_219076_20021029_20170127_01_T1_B4.tif

names : LE07_L1TP_219076_20021029_20170127_01_T1_B4

values : 0, 255 (min, max)

What I would like to do after this is open the other band and proceed calculating a vegetation index, like NDVI. But I can't even open properly them.

Best Answer

gdalinfo on the file reports the following transformation:

GeoTransform =
  259784.9999996533, 30.0000000002095, -4.662376007821301e-11
  7491805.000000283, -1.623406322724229e-10, -29.99999999998893

To six decimal places those numbers are either integers or zero, and would correspond to an unrotated raster with the expected coordinates. So any transformation (rotation/skew) is probably negligible.

When I read your raster in I can't even plot it without an error:

> r = raster("./LE07_L1TP_219076_20021029_20170127_01_T1_B4.tif")
Warning message:
In .rasterFromGDAL(x, band = band, objecttype, ...) : 

 This file has a rotation
 Support such files is limited and results of data processing might be wrong.
 Proceed with caution & consider using the "rectify" function

> plot(r)
Error in if (nc < 1) { : missing value where TRUE/FALSE needed
In addition: Warning message:
In setExtent(out, e, keepres = TRUE) :
  NAs introduced by coercion to integer range

The rotation matrix parameters are here:

> r@rotation@geotrans
         ll.x         res.x     oblique.x          <NA>     oblique.y 
 2.597850e+05  3.000000e+01  4.662376e-11  7.491805e+06  1.623406e-10 
        res.y 
-3.000000e+01

and if you believe these are negligible then mark the raster as unrotated:

> r@rotated=FALSE
> plot(r)
>

and it looks like its in the right place. I don't know if QGIS is actually doing the rotation correctly, or ignoring it, I'd need a raster with a bigger rotation to figure that out! I do get this message from QGIS:

 Warning: Creating Warped VRT.

which I think is what it says when it has to do a transformation.

But for your raster and any more with near-identical transformations, you can set the rotate flag to false.

The reason it fails in R is because res is not getting the resolution (cell size) correctly. Look:

 resolution : 4.662376e-11, 1.623406e-10 (x, y)

and that is because the res function does this:

 if (rotated(x)) {
    return(x@rotation@geotrans[c(3, 5)])
 }

and the geotrans slot is this:

> r@rotation@geotrans
         ll.x         res.x     oblique.x          <NA>     oblique.y 
 2.597850e+05  3.000000e+01  4.662376e-11  7.491805e+06  1.623406e-10 
        res.y 
-3.000000e+01

so res should look at elements 2 and 6 to get the resolution, not 3 and 5.

Related Question