Probleme mit NA-Werten beim Lesen der DEM-Datei mit dem R-Raster-Paket in Windows

Ich versuche, eine Rasterdatei im .DEM-Format unter Windows mit dem 'Raster'-Paket in R zu lesen.

Ich bekomme Probleme mit NA-Werten, wenn ich die Daten in R in Windows 7 lade, aber ich habe das Problem auf einem Mac mit OSX Lion nicht. Unter Windows scheinen die NA-Werte nicht richtig gelesen zu werden. Die Frage ist, warum das passiert?

Die verwendete Rasterdatei wurde von USGS mit dem folgenden R-Code heruntergeladen:

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

Dann habe ich das Raster mit dem 'Raster'-Paket in R eingelesen. In OSX Lion und R64 Version 2.13.1 werden die NA-Werte erkannt:

> 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

Unter Windows 7 (64 Bit, dieselbe R-Version) werden die Zellenwerte, die NAs sein sollten, in Zahlen konvertiert:

> 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

Warum enthält das Raster keine NA-Werte, wenn ich es unter Windows lese? Wie könnte ich das umgehen? Ich vermute, es hat mit der Art und Weise zu tun, wie Zahlen gespeichert werden. Viele der NA-Werte werden in 55540 konvertiert.

Info von Windows (nach dem Laden des Rasters):

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

Infos von OSX (nach dem Laden des Rasters):

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

Eine Problemumgehung besteht darin, nur die Rohdaten zu erfassen, da dies ein sehr einfaches Dateiformat ist.

Nicht für jedermann, aber es kann aufschlussreich sein, zu sehen, was passiert.

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

An diesem Punkt können Sie die verschiedenen Optionen für Ganzzahlzeichen und Endianness direkt ausprobieren. Wenn Sie auf diese Weise lesen, erreichen Sie, was Robert mit der > 32767Transformation macht, nachdem die Datei gelesen wurde.

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)

Stellen Sie abschließend die Projektion so ein, wie sie vom Raster gelesen wird (und dies würde das gleiche Seitenverhältnis in der Darstellung ergeben, das beim Lesen auf diese Weise angezeigt wird).

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

EDIT: Hoppla, hatte vergessen, von oben zu subtrahieren, jetzt behoben - es gibt immer noch ein Problem mit halben Zellen, dem ich auch nicht auf den Grund gegangen bin.

Es gibt einige Probleme mit dieser Datei oder mit GDAL. Ich benutze Windows 7

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

und

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

Beachten Sie, dass der NoDataValue mit dem ungeraden Bmin-Wert (-9999) identisch ist. Schlimmer ist, dass GDType UInt16 ist - vorzeichenlose 2-Byte-Ganzzahlen - was bedeutet, dass Sie keine Werte unter Null haben können. Dies ist wahrscheinlich ein Fehler, der in GDAL 1.8.0 behoben wurde

Das Problem wird dargestellt, wenn Sie dies tun

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

Ich denke, der schnellste Weg, dies zu beheben, ist:

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

Das Problem scheint durch ein Problem verursacht zu werden, das die Tatsache erkennt, dass die Daten im vorzeichenbehafteten 2-Byte-Ganzzahlformat vorliegen. Es wird fälschlicherweise als vorzeichenloses 2-Byte-Ganzzahlformat interpretiert. Daher lautet Ihr Knotenwert von -9999: 2 Byte = 256 * 256 -9999 = 55537

Was ich seltsam finde, ist, dass der Min-Wert: -9999 und der Max-Wert: 5483 für Windows und Mac gleich sind. Es scheint, dass in beiden Fällen beim Erstellen der Header keine Daten korrekt identifiziert wurden, aber bei der tatsächlichen Verwendung für die Werte ist ein Fehler aufgetreten.

Problemumgehung:

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

Um tiefer zu graben: Es scheint, dass Raster rgdal aufruft, das wiederum gdal selbst aufruft. Höchstwahrscheinlich haben Sie eine andere Version von gdal auf Ihrem System. Überprüfen Sie beim Laden von rgdal zB:

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

Ich habe gerade eine kurze Überprüfung unter Linux durchgeführt: gdal 1.8 lädt die Datei einwandfrei, aber gdal 1.6 schlägt fehl. Es scheint also von gdal verursacht zu werden.

Obwohl ich mir über Ihre Anforderung nicht sicher bin, können Sie konvertieren. DEM-Dateien in .GRID-Dateien. Der Arcgis-Geoprozessor oder R erkennt dann automatisch .GRIDs mit N / A-Werten während der Manipulation des Rasterrasters.