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

Wed Feb 10 18:59:32 CET 2010

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