Skip to content

Reversed raster after readGDAL()

9 messages · Agustin Lobo, Roger Bivand, Sarah Goslee +1 more

Original
Roger Bivand
#
On Wed, 10 Feb 2010, Agustin Lobo wrote:

            

      
      
      
    

        
Collected from:

http://sites.google.com/site/eospansite/dummy

Yes, it is flipped on the Y axis. What software made it? Could you please 
show the images with axes, as I think that the GTiff has a wrong sign on 
its y-step? What are the coordinates of the NW and SW corners of the image 
in QGIS? GDAL expects images to be read by row N to S, so a y-step with 
the wrong sign (+ instead of -) flips the image northward.

Roger

PS:

gdalinfo externally has lower left north of upper left:

            

      
      
      
    

        
Driver: GTiff/GeoTIFF
Files: test2.tif
Size is 256, 256
Coordinate System is:
PROJCS["ED50 / UTM zone 31N",
     GEOGCS["ED50",
         DATUM["European_Datum_1950",
             SPHEROID["International 1924",6378388,297.0000000000014,
                 AUTHORITY["EPSG","7022"]],
             AUTHORITY["EPSG","6230"]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4230"]],
     PROJECTION["Transverse_Mercator"],
     PARAMETER["latitude_of_origin",0],
     PARAMETER["central_meridian",3],
     PARAMETER["scale_factor",0.9996],
     PARAMETER["false_easting",500000],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]],
     AUTHORITY["EPSG","23031"]]
Origin = (424389.000000000000000,4635822.000000000000000)
Pixel Size = (345.007812500000000,194.019531250000000)
Metadata:
   AREA_OR_POINT=Area
Image Structure Metadata:
   INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  424389.000, 4635822.000) (  2d 5'20.10"E, 41d52'11.88"N)
Lower Left  (  424389.000, 4685491.000) (  2d 4'56.97"E, 42d19'2.06"N)
Upper Right (  512711.000, 4635822.000) (  3d 9'11.42"E, 41d52'24.52"N)
Lower Right (  512711.000, 4685491.000) (  3d 9'15.31"E, 42d19'14.90"N)
Center      (  468550.000, 4660656.500) (  2d37'10.89"E, 42d 5'47.83"N)
Band 1 Block=256x4 Type=Float64, ColorInterp=Gray

and:

            

      
      
      
    

        
rows        256
columns     256
bands       1
origin.x        424389
origin.y        4636016
res.x       345.0078
res.y       194.0195
oblique.x   0
oblique.y   0
driver      GTiff
projection  +proj=utm +zone=31 +ellps=intl +units=m +no_defs
file        test2.tif
apparent band summary:
    GDType        Bmin       Bmax
1 Float64 -4294967295 4294967295
Metadata:
AREA_OR_POINT=Area

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
I think I tripped on this same stone, but cannot find anything in my
records, perhaps it's on the lost computer.

The file test2.tif was made by plugin Interpolation in QGIS out of a
vector layer of points with nitrate concentration values, just
a simple inverse distance method. I think that the author f that
pluginis Marco Hugentobler (Marco, sorry if it's not you) to whom I'm
forwarding.

Relevant info from the metadata in QGIS for this raster:

Layer Spatial Reference System:
+proj=utm +zone=31 +ellps=intl +units=m +no_defs

Origin
424389,4.68549e+06

Pixel Size:
279.887,-279.887

NW corner in QGIS (interactive):424388.820,4685490.931
SW corner in QGIS (interactive):424389.122,4635951.740

Image with grid: http://sites.google.com/site/eospansite/dummy/test2gridqgis.jpg

The weird thing is that QGIS uses gdal to read the geotif files.


Agus

2010/2/10 Roger Bivand <Roger.Bivand at nhh.no>

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Sarah Goslee
    

    

      
      #
    

  


  

      
This may well be too simplistic an answer, but it might just be image().

            

      
      
      
    

        
Notice that ?image? interprets the ?z? matrix as a table of
     ?f(x[i], y[j])? values, so that the x axis corresponds to row
     number and the y axis to column number, with column 1 at the
     bottom, i.e. a 90 degree counter-clockwise rotation of the
     conventional printed layout of a matrix.

x <- matrix(rep(1:10, each=10), nrow=10)
x
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1    2    3    4    5    6    7    8    9    10
 [2,]    1    2    3    4    5    6    7    8    9    10
 [3,]    1    2    3    4    5    6    7    8    9    10
 [4,]    1    2    3    4    5    6    7    8    9    10
 [5,]    1    2    3    4    5    6    7    8    9    10
 [6,]    1    2    3    4    5    6    7    8    9    10
 [7,]    1    2    3    4    5    6    7    8    9    10
 [8,]    1    2    3    4    5    6    7    8    9    10
 [9,]    1    2    3    4    5    6    7    8    9    10
[10,]    1    2    3    4    5    6    7    8    9    10
image(x) # not what you might expect
image(t(x)) # this view matches the matrix representation

What happens if you use plot(mygeotiff) instead of image(mygeotiff)?

Sarah

        
On Wed, Feb 10, 2010 at 11:17 AM, Agustin Lobo <alobolistas at gmail.com> wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Paul Hiemstra
    

    

      
      #
    

  


  

      
Hi,

What if you make an spplot, still the same problem?

cheers,
Paul

        
Agustin Lobo wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
Thanks Sarah. But I think that package sp has an specific method of
image() for objects SpatialGridDataFrame, which is the class returned
by readGDAL(), and this specific method takes care of this problem. Check
image.SpatialGridDataFrame in package {sp}

Agus

2010/2/10 Sarah Goslee <sarah.goslee at gmail.com>:

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
Yes, same problem
Agus

2010/2/10 Paul Hiemstra <p.hiemstra at geo.uu.nl>:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 10 Feb 2010, Agustin Lobo wrote:

        

            

      
      
      
    

        
Sorry, can't find such a posting - maybe my search keys were wrong - you 
could try the nabble archive.

            

      
      
      
    

        
No, it is 345.0078, 194.0195.

            

      
      
      
    

        
SW in R is 424389, 4635822
or SW cell centre 424561.5, 4635919.0 which is the SW corner plus 0.5* the 
resolution.

Try:

GDALinfo(system.file("pictures/cea.tif", package = "rgdal")[1])
system(paste("gdalinfo", system.file("pictures/cea.tif", package =
   "rgdal")[1]))

to see the oddness of gdalinfo on your raster. The SW (LL) corner is north 
of the NW (UL) one. The provided raster is OK, so I guess one might think 
that yours isn't?

            

      
      
      
    

        
How is the plugin writing the GTiff?

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Agustin Lobo
    

    

      
      #
    

  


  

      
I forward to qgis devel and Marco Hugentobler.
I really cannot understand, Qgis uses gdal for writing and reading geotif files
as far as I know.
I think the point is that QGIS reports the origin in the NW corner and then uses
negative resolution for Y. But do not understand the discrepancy in
the absolute values
of the pixel size that you report.

I'm going to check if this is just for this plugin or for any raster
written by QGIS
as geotiff.

Agus

2010/2/10 Roger Bivand <Roger.Bivand at nhh.no>: