1

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.

  • gdalinfo reports nothing strange, and the file loads into QGIS at the right spot. – AndreJ Oct 23 '17 at 6:22
3

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.

  • Thank you very much! I asked my teacher how the image was created and it was a LANDSAT 7 image which was cropped using Python with gdal_translate. Then, using Python and gdalmove it was changed from '-s_srs', 'EPSG:32623' (North) to '-t_srs', 'EPSG:32723' (South). Maybe this caused the trouble? – Kol Rocket Oct 23 '17 at 13:13
  • The trouble is in the raster package's handling of the transformation matrix. I've raised an issue in the raster package issue tracker. – Spacedman Oct 23 '17 at 13:16
  • 1
    This was fixed in version 2.6-1 (forthcoming) – Robert Hijmans Oct 31 '17 at 19:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.