[R-sig-Geo] Reversed raster after readGDAL()

[Roger Bivand]

On Wed, 10 Feb 2010, Agustin Lobo wrote:

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

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

>
> 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
>

No, it is 345.0078, 194.0195.

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

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?

> 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.

How is the plugin writing the GTiff?

Roger



>
>
> Agus
>
> 2010/2/10 Roger Bivand <Roger.Bivand at nhh.no>
>>
>> On Wed, 10 Feb 2010, Agustin Lobo wrote:
>>
>>> Hi!
>>>
>>> After
>>>
>>>> r3 <- readGDAL("test2.tif")
>>>
>>> test2.tif has GDAL driver GTiff
>>> and has 256 rows and 256 columns
>>>
>>>> image(r3,col=rainbow(128))
>>>
>>> I get the image N-S reversed.
>>>
>>> I've put an screenshot and the actual geotif image here:
>>>
>>> http://sites.google.com/site/eospansite/dummy/test2_screenshot2.jpeg
>>> http://sites.google.com/site/eospansite/dummy/test2.tif
>>
>> 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:
>>
>>> system("gdalinfo test2.tif")
>>
>> 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:
>>
>>> GDALinfo("test2.tif")
>>
>> 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
>>
>>
>>>
>>> Agus
>>>
>>>        [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> --
>> Roger Bivand
>> Economic Geography Section, Department of Economics, Norwegian School of
>> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
>> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
>> e-mail: Roger.Bivand at nhh.no
>>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no