Reversed raster after readGDAL()

classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

Reversed raster after readGDAL()

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

Agus

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Reversed raster after readGDAL()

Roger Bivand
Administrator
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
> [hidden email]
> 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: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Fwd: Reversed raster after readGDAL()

Agustin Lobo
In reply to this post by Agustin Lobo
---------- Forwarded message ----------
From: Agustin Lobo <[hidden email]>
Date: 2010/2/10
Subject: Re: [R-sig-Geo] Reversed raster after readGDAL()
To: [hidden email]


It is, but the link does not work (don't know why)
Write http://sites.google.com/site/eospansite/dummy
in your browser and you'll see the file, at the bottom of the page,
the third one in reverse order.

Agus

2010/2/10 Roger Bivand <[hidden email]>

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
>>
>
> The second file is not available.
>
> Roger
>
>
>> Agus
>>
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> 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: [hidden email]
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Reversed raster after readGDAL()

Agustin Lobo
In reply to this post by Roger Bivand
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 <[hidden email]>

>
> 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
>> [hidden email]
>> 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: [hidden email]
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Reversed raster after readGDAL()

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

>From the help for 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 <[hidden email]> wrote:

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


--
Sarah Goslee
http://www.functionaldiversity.org

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Reversed raster after readGDAL()

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 <[hidden email]>:

> This may well be too simplistic an answer, but it might just be image().
>
> From the help for 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 <[hidden email]> wrote:
>> 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
>>
>
>
> --
> Sarah Goslee
> http://www.functionaldiversity.org
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Reversed raster after readGDAL()

Paul Hiemstra
Hi,

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

cheers,
Paul

Agustin Lobo wrote:

> 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 <[hidden email]>:
>  
>> This may well be too simplistic an answer, but it might just be image().
>>
>> From the help for 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 <[hidden email]> wrote:
>>    
>>> 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
>>>
>>>      
>> --
>> Sarah Goslee
>> http://www.functionaldiversity.org
>>
>>    
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>  


--
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Reversed raster after readGDAL()

Agustin Lobo
Yes, same problem
Agus

2010/2/10 Paul Hiemstra <[hidden email]>:

> Hi,
>
> What if you make an spplot, still the same problem?
>
> cheers,
> Paul
>
> Agustin Lobo wrote:
>>
>> 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 <[hidden email]>:
>>
>>>
>>> This may well be too simplistic an answer, but it might just be image().
>>>
>>> From the help for 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 <[hidden email]>
>>> wrote:
>>>
>>>>
>>>> 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
>>>>
>>>>
>>>
>>> --
>>> Sarah Goslee
>>> http://www.functionaldiversity.org
>>>
>>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
> --
> Drs. Paul Hiemstra
> Department of Physical Geography
> Faculty of Geosciences
> University of Utrecht
> Heidelberglaan 2
> P.O. Box 80.115
> 3508 TC Utrecht
> Phone:  +3130 274 3113 Mon-Tue
> Phone:  +3130 253 5773 Wed-Fri
> http://intamap.geo.uu.nl/~paul
>
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Reversed raster after readGDAL()

Roger Bivand
Administrator
In reply to this post by Agustin Lobo
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 <[hidden email]>
>>
>> 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
>>> [hidden email]
>>> 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: [hidden email]
>>
>
--
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: [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: Reversed raster after readGDAL()

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 <[hidden email]>:

> 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 <[hidden email]>
>>>
>>> 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
>>>> [hidden email]
>>>> 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: [hidden email]
>>>
>>
>
> --
> 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: [hidden email]
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo