Problèmes avec les valeurs NA lors de la lecture du fichier .DEM avec le package R 'raster' dans Windows

J'essaie de lire un fichier raster au format .DEM sur Windows en utilisant le package «raster» dans R.

J'ai des problèmes avec les valeurs NA, lors du chargement des données dans R dans Windows 7, mais je n'ai pas le problème sur un Mac avec OSX Lion. Sous Windows, les valeurs NA ne semblent pas être lues correctement. La question est pourquoi cela se produit-il?

Le fichier raster utilisé a été téléchargé depuis USGS avec le code R suivant:

download.file('http://edcftp.cr.usgs.gov/pub/data/gtopo30/global/e020n90.tar.gz', 'e020n90.tar.gz')
untar('e020n90.tar.gz')

Ensuite, j'ai lu le raster dans R en utilisant le package «raster». Dans OSX Lion et R64 version 2.13.1, les valeurs NA sont reconnues:

> onMac <- raster('E020N90.DEM')
> onMac
class       : RasterLayer 
dimensions  : 6000, 4800, 28800000  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 20, 60, 40, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
values      : /Users/Tam/Desktop/E020N90.DEM 
min value   : -9999 
max value   : 5483 

> summary(values(onMac))
Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
-137       85      148      213      213     5483 13046160

Mais sur Windows 7 (64Bit, même version R), il convertit les valeurs de cellule qui devraient être NA en nombres:

> onWindows <- raster('E020N90.DEM')
> onWindows
class       : RasterLayer 
dimensions  : 6000, 4800, 28800000  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 20, 60, 40, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
values      : E:/WorldDegreeDays/gsoddata/gtopo/E020N90.DEM 
min value   : -9999 
max value   : 5483 

> summary(values(onWindows))
Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1     150     946   27190   55540   65540

Pourquoi n'y a-t-il pas de valeurs NA dans le raster lorsque je le lis sur Windows? Comment pourrais-je contourner ce problème? Je suppose que cela a à voir avec la façon dont les nombres sont stockés, une grande partie des valeurs NA sont converties en 55540.

Informations de Windows (après le chargement du raster):

SessionInfo()
R version 2.13.1 (2011-07-08)
Platform: x86_64-pc-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rgdal_0.7-1   raster_1.9-12 sp_0.9-88    

loaded via a namespace (and not attached):
[1] grid_2.13.1     lattice_0.19-30

Informations d'OSX (après le chargement du raster):

R version 2.13.1 (2011-07-08)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] rgdal_0.6-33  raster_1.9-12 sp_0.9-88    

loaded via a namespace (and not attached):
[1] grid_2.13.1     lattice_0.19-33

Une solution consiste simplement à opter pour les données brutes, car il s'agit d'un format de fichier très simple.

Pas pour tout le monde, mais il peut être éclairant de voir ce qui se passe.

## all these details are in the .HDR file
NROWS   <-      6000
NCOLS   <-      4800

À ce stade, vous pouvez essayer les différentes options pour le signe entier et l'endianité directement, et en lisant de cette façon, nous réalisons ce que Robert fait avec la > 32767transformation après la lecture du fichier.

x1 <- readBin("E020N90.DEM", "integer", size = 2, signed = TRUE, n = NROWS * NCOLS, endian = "big")

range(x1)
[1] -9999  5483

x1[x1 < -9998] <- NA

## now for the simple georeferencing, also in the HDR file

ULXMAP   <-     20.00416666666667
ULYMAP   <-     89.99583333333334
XDIM     <-     0.00833333333333
YDIM     <-     0.00833333333333

## now generate x/y coordinates, and the data matrix (flip on Y)
x <- list(x = seq(ULXMAP, by = XDIM, length = NCOLS),
       y = seq(ULYMAP - NROWS * YDIM, by = YDIM, length = NROWS),
      z = matrix(x1, nrow = NCOLS)[ , NROWS:1])

library(sp)

x <- image2Grid(x)

library(raster)
r <- raster(x)

plot(r)

Enfin, définissez la projection telle qu'elle est lue par raster (et cela donnerait le même rapport hauteur / largeur dans l'intrigue qui est vu lorsqu'il est lu de cette façon).

projection(r) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

EDIT: Oups, j'avais oublié de soustraire du haut, maintenant corrigé - il y a toujours un problème de demi-cellule que je ne suis pas allé aussi loin.

Il y a quelques problèmes avec ce fichier ou avec GDAL. J'utilise Windows 7

R version 2.13.1 (2011-07-08)
Platform: x86_64-pc-mingw32/x64 (64-bit)

et

> getGDALVersionInfo()
[1] "GDAL 1.7.2, released 2010/04/23"


> GDALinfo('E020N90.DEM')
rows        6000 
columns     4800 
bands       1 
origin.x        20 
origin.y        40 
res.x       0.008333333 
res.y       0.008333333 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      EHdr 
projection  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs 
file        E020N90.DEM 
apparent band summary:
 GDType  Bmin Bmax   Bmean    Bsd hasNoDataValue NoDataValue
1 UInt16 -9999 5483 -4412.9 5088.6           TRUE       -9999
> 

Notez que le NoDataValue est identique à la valeur Bmin (-9999), ce qui est impair. Le pire est que GDType est UInt16 - Entiers non signés sur 2 octets - ce qui signifie que vous ne pouvez pas avoir de valeurs inférieures à zéro. C'est probablement un bug qui a été corrigé dans gdal 1.8.0

Le problème est illustré lorsque vous faites

r <- 'E020N90.DEM'
plot(r)

Je pense que la façon la plus rapide de résoudre ce problème est la suivante:

r <- raster('E020N90.DEM')
fun <- function(x){ x[x > 32767] <- x[x > 32767] - 65536; x[x == -9999] <- NA; x}
r[] <- fun(values(r))

plot(r)
r <- writeRaster(r, 'E020N90.TIF')

Le problème semble provenir d'un problème de reconnaissance du fait que les données sont au format entier signé de 2 octets. Il est interprété à tort comme un format entier de 2 octets non signé. Par conséquent, votre valeur nodale de -9999 devient: 2 octets = 256 * 256 -9999 = 55537

Ce que je trouve étrange, c'est que la valeur min: -9999 et la valeur max: 5483 sont les mêmes pour Windows et Mac. Il semble que dans les deux cas, aucune donnée n'ait été correctement identifiée lors de la création des en-têtes, mais lors de leur utilisation réelle pour les valeurs, une erreur s'est produite.

solution de contournement:

values(onWindows)[values(onWindows)>128*256]<-values(onWindows)[values(onWindows)>128*256]-256*256
values(onWindows)[values(onWindows)==-9999]<-NA

Pour creuser plus profondément: Il semble que raster appelle rgdal, qui à son tour appelle gdal lui-même. Vous avez probablement une version différente de gdal sur votre système. Vérifiez lors du chargement de rgdal, par exemple:

Loaded GDAL runtime: GDAL 1.8.0, released 2011/01/12

Je viens de faire une vérification rapide sur linux: gdal 1.8 charge bien le fichier, mais gdal 1.6 échoue. Donc, cela semble être causé par gdal.

Bien que je ne sois pas sûr de vos besoins, vous pouvez convertir. Fichiers DEM en fichiers .GRID. Ensuite, le géoprocesseur arcgis ou R reconnaîtra automatiquement les .GRID avec des valeurs N / A lors de la manipulation du raster de la grille.