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

10

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
yellowcap
la source
version raster 1.9-12 sur les deux systèmes
yellowcap
Pouvez-vous inclure votre sessionInfo()dans votre message?
Roman Luštrik
J'ai obtenu des valeurs différentes sur raster_1.8-12 (mais identiques aux vôtres sur 1.9-12) sur winXP.
Roman Luštrik
Cela a-t-il fonctionné correctement avec raster_1.8-12, ou était-ce juste différent?
yellowcap

Réponses:

11

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)

entrez la description de l'image ici

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.

mdsumner
la source
En fait, vous pouvez combiner les deux méthodes (cette réponse et mes réponses / roberts): r <- raster('E020N90.DEM')puis exécutez values(r)<-readBin("E020N90.DEM", "integer", size = 2, signed = TRUE, n = nrows(r) * ncols(r), endian = "big")puis values(r)[values(r)==-9999]<-NA
johanvdw
Ha oui mais c'est de l'hérésie
mdsumner
6

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')
RobertH
la source
1
Ce correctif est meilleur que le mien car les points de données dans la mer caspienne sont également convertis (ces points sont également négatifs). Agréable!
johanvdw
6

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.

johanvdw
la source
Runtime GDAL chargé: GDAL 1.7.2, publié le 23/04/2010
Roman Luštrik
Sur Windows, ma version GDAL est également celle citée ci-dessus (1.7.2.), Sur OSX j'ai 1.8.0. Mais pourquoi ne puis-je pas lire le fichier DEM en utilisant 1.7.2.? Y a-t-il une solution de contournement?
yellowcap
J'ai obtenu des résultats différents dans différentes versions de raster (voir mes commentaires ci-dessus), donc je ne suis pas totalement convaincu que c'est GDAL en soi .
Roman Luštrik
Pouvez-vous décrire comment rgdaltrouver une gdalinstallation mise à jour sur Win7? J'ai téléchargé et installé les gdalbinaires les plus récents (32 et 64). Ceux-ci ont été installés à l'emplacement par défaut mais rgdalutilisent toujours la 1.7.2, même après la mise à jour.
yellowcap
La mise à jour de rgdal n'est pas évidente et nécessitera une recompilation de rgdal. Plus d'infos ici .
johanvdw
0

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.

EvilInside
la source
L'utilisation d'un autre logiciel pour convertir le fichier en premier est possible mais pas ce que je voulais. L'idée était de n'utiliser R que pour télécharger, lire et analyser le fichier.
yellowcap
en principe, vous pouvez exécuter gdaltranslate via R en utilisant system2.
johanvdw