[R-sig-Geo] Rgdal and unsigned integers

[Roger Bivand]

On Thu, 20 Aug 2009, Robert J. Hijmans wrote:

> Jonathan,
>
> I used the file you send me. This shows that the problem is already
> present in rgdal (i.e. before it gets passed to raster)


The next release of rgdal will contain some tools to permit the handling 
of signed/unsigned integer conversion (for 8 and 16 bit data so far), 
where GDAL drivers are unhelpful (submitted to CRAN, available for source 
check-out from the sourceforge rgdal repository.

In this case, the LAN/GIS driver assumes 16 bit signed integers when the 
data are 16 bit, and there are no "authorities" to say that signed or 
unsigned are correct. Here, the new helper function toUnSigned(x, 16) (for 
x as the band of interest) will do the conversion.

GDALinfo() now reports band-wise the GDAL Data Type (GDT), and the band 
minimum and maximum values known to the driver (either based on metadata, 
which varies with driver - maybe in an auxiliary file, or on the 
characteristics of the GDT, so just the range of the data representation, 
not the data in the band) It will not read the band to find the values.

The toSigned() helper function may be used for example when the data are 
read as 8 bit unsigned (the only option for GDT Byte), but ought to be 
signed.

Does anyone need conversion for 32 bit integers? This will overflow if 
not coerced to numeric 8-byte representation with storage.mode "double", 
as R integers are signed.

In summary, if the data look de-ranged, use GDALinfo() to see what GDT the 
driver assigned, and if appropriate, convert signed/unsigned to suit.

Hope this helps,

Roger

PS. Don't be thrown if the summary() of the data shows rounded maximum or 
minimum values - summary() takes a digits= argument with a tight default, 
because the output is usually read for its most significant (left) digits, 
not the rightmost digits. If in doubt, use range().

> library(rgdal)
> GDALinfo("stand_bm.gis")
rows        839
columns     1758
bands       1
origin.x        33773.26
origin.y        861578
res.x       100
res.y       100
oblique.x   0
oblique.y   0
driver      LAN
projection  NA
file        stand_bm.gis
apparent band summary:
   GDType   Bmin  Bmax
1  Int16 -32768 32767
> xx <- readGDAL("stand_bm.gis")
stand_bm.gis has GDAL driver LAN
and has 839 rows and 1758 columns
> summary(xx$band1)
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  -32140       0       0    8316   19180   29620
> summary(toUnSigned(xx$band1, 16))
    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
       0       0   14250   12830   22230   48500
> summary(toUnSigned(xx$band1, 16), digits=7)
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
     0.00     0.00 14247.00 12833.83 22226.00 48499.00
> range(toUnSigned(xx$band1, 16))
[1]     0 48499


>
>> x <- readGDAL("stand_bm.gis")
> stand_bm.gis has GDAL driver LAN
> and has 839 rows and 1758 columns
>> summary(x)
> (...)
> Data attributes:
>   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
> -32140       0       0    8316   19180   29620
>>
>
> I opened the file in ArcMap, saved it to a esri grid and read it into R
>
>> stndesri <- raster("esri_stand", values = T)
>
>> stnd <- raster("stand_bm.gis", values = T)
>
> I adjusted the adjustment I send before
>> stnd[stnd < 0] <- stnd[stnd < 0] + 65536
>
> and these look good:
>
>> stnd
> class       : RasterLayer
> filename    :
> nrow        : 839
> ncol        : 1758
> ncells      : 1474962
> min value   : 0
> max value   : 48499   <--- same is gdal
> (...)
>
> the values of the two files are the same:
>> dif <- stnd - stndesri
>
>> summary(dif)
> Values
> Min.    0
> 1st Qu. 0
> Median  0
> Mean    0
> 3rd Qu. 0
> Max.    0
>
>
> The values are still very high in some places, but they should be,
> that's what caused the problem to begin with. Whether these high
> values are meaningful is another matter; but that could be caused by
> LANDIS if that's what you are using.
>
> Apparently something needs to be fixed in rgdal as Barry showed that
> gdal gets the right values; it seems somewhere along the line the
> INT16 is interpreted as signed, whereas it should be unsigned.
>
> Best,
> Robert
>
> On Thu, Aug 20, 2009 at 2:44 PM, Barry
> Rowlingson<b.rowlingson at lancaster.ac.uk> wrote:
>> On Thu, Aug 20, 2009 at 9:57 PM, Jonathan
>> Thompson<jthomps at fas.harvard.edu> wrote:
>>
>>> gdalinfo:
>>> Band 1 Block=1758x1 Type=Int16, ColorInterp=Undefined
>>>  Min=0.000 Max=48499.000
>>
>>>> Arcgis:
>>>> Pixel Type  = unsigned integer
>>>> Pixel Depth = 16 bit
>>
>> Arcgis is saying 'unsigned integer' but gdalinfo says 'Int16'. But
>> then gdalinfo says the max is 48499 - which shouldn't be possible with
>> (signed) Int16 values (they go from -32k to +32k).
>>
>>  So something is very odd...
>>
>> Barry
>>
>> _______________________________________________
>> 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