Problemy z wartościami NA podczas odczytu pliku .DEM z pakietem R „raster” w systemie Windows

Próbuję odczytać plik rastrowy w formacie .DEM w systemie Windows za pomocą pakietu „raster” w języku R.

Występują problemy z wartościami NA podczas ładowania danych do R w systemie Windows 7, ale nie mam problemu na komputerze Mac z systemem OSX Lion. W systemie Windows wartości NA nie są poprawnie odczytywane. Pytanie brzmi, dlaczego tak się dzieje?

Używany plik rastrowy został pobrany z USGS z następującym kodem R:

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

Następnie czytam raster do R za pomocą pakietu „raster”. W OSX Lion i R64 wersja 2.13.1 rozpoznawane są wartości NA:

> 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

Ale w Windows 7 (64Bit, ta sama wersja R) konwertuje wartości komórek, które powinny być NA, na liczby:

> 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

Dlaczego w rastrze nie ma żadnych wartości NA, gdy czytam je w systemie Windows? Jak mogłem to obejść? Domyślam się, że ma to związek ze sposobem przechowywania liczb, wiele wartości NA jest konwertowanych na 55540.

Informacje z systemu Windows (po załadowaniu rastra):

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

Informacje z OSX (po załadowaniu rastra):

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

Jednym obejściem jest po prostu użycie surowych danych, ponieważ jest to bardzo prosty format pliku.

Nie dla wszystkich, ale może być pouczające, aby zobaczyć, co się dzieje.

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

W tym momencie możesz wypróbować różne opcje znaku liczby całkowitej i endianizmu, a czytając w ten sposób osiągamy to, co Robert robi z > 32767transformacją po odczytaniu pliku.

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)

Na koniec ustaw projekcję tak, jak jest odczytywana przez raster (a to dałoby ten sam współczynnik kształtu na wykresie, który jest widziany podczas odczytu w ten sposób).

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

EDYCJA: Ups, zapomniałem odjąć od góry, teraz naprawione - wciąż jest problem z półkomórką, do którego też nie dotarłem.

Istnieją pewne problemy z tym plikiem lub GDAL. Korzystam z systemu Windows 7

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

i

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

Zauważ, że wartość NoDataValue jest taka sama jak wartość Bmin (-9999), co jest nieparzyste. Co gorsza, GDType to UInt16 - 2-bajtowe liczby całkowite bez znaku - co oznacza, że ​​nie możesz mieć wartości niższych niż zero. Jest to prawdopodobnie błąd, który został naprawiony w gdal 1.8.0

Problem jest ilustrowany, kiedy to robisz

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

Myślę, że najszybszym sposobem na rozwiązanie tego jest:

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

Problem wydaje się być spowodowany problemem polegającym na rozpoznaniu faktu, że dane są w formacie 2-bajtowej liczby całkowitej. Jest źle interpretowany jako 2-bajtowy format liczb całkowitych bez znaku. Dlatego twoja wartość nodata -9999 wynosi: 2 bajty = 256 * 256 -9999 = 55537

Dziwne jest to, że wartość minimalna: -9999 i wartość maksymalna: 5483 są takie same dla systemu Windows i Mac. Wydaje się, że w obu przypadkach żadne dane nie zostały poprawnie zidentyfikowane podczas budowania nagłówków, ale podczas faktycznego używania ich do wartości wystąpił błąd.

obejście:

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

Aby kopać głębiej: Wygląda na to, że raster wywołuje rgdal, co z kolei wywołuje samą gdal. Najprawdopodobniej masz w systemie inną wersję gdal. Sprawdź podczas ładowania rgdal np .:

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

Właśnie zrobiłem szybkie sprawdzenie linuksa: gdal 1.8 ładuje plik dobrze, ale gdal 1.6 nie działa. Więc wydaje się, że jest to spowodowane przez gdal.

Chociaż nie jestem pewien co do twoich wymagań, możesz dokonać konwersji. Pliki DEM do plików .GRID. Następnie geoprocesor arcgis lub R automatycznie rozpoznają .GRID z wartościami N / A podczas manipulacji rastrem siatki.