Masalah dengan nilai NA saat membaca file .DEM dengan paket R 'raster' di Windows

Saya mencoba membaca file raster dalam format .DEM di windows menggunakan paket 'raster' di R.

Saya mendapatkan masalah dengan nilai-nilai NA, ketika memuat data ke R di Windows 7, tapi saya tidak punya masalah di Mac dengan OSX Lion. Di windows, nilai-nilai NA tampaknya tidak dibaca dengan benar. Pertanyaannya adalah mengapa ini terjadi?

File raster yang digunakan diunduh dari USGS dengan kode R berikut:

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

Lalu saya membaca raster ke R menggunakan paket 'raster'. Dalam OSX Lion dan R64 versi 2.13.1, nilai-nilai NA diakui:

> 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

Tetapi pada Windows 7 (64Bit, versi R yang sama) itu mengubah nilai sel yang harus menjadi angka-angka NA:

> 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

Mengapa tidak ada nilai NA di raster ketika saya membacanya di Windows? Bagaimana saya bisa mengatasinya? Dugaan saya ada hubungannya dengan cara nomor disimpan, banyak nilai NA dikonversi ke 55540.

Info dari Windows (setelah memuat 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

Info dari OSX (setelah memuat 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

Satu solusi hanya untuk mencari data mentah, karena ini adalah format file yang sangat sederhana.

Tidak untuk semua orang, tetapi bisa menjadi cahaya untuk melihat apa yang terjadi.

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

Pada titik ini Anda dapat mencoba opsi yang berbeda untuk tanda integer dan endianness secara langsung, dan membaca dengan cara ini kita mencapai apa yang Robert lakukan dengan > 32767transformasi setelah file dibaca.

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)

Akhirnya, tetapkan proyeksi seperti yang dibaca oleh raster (dan ini akan memberikan rasio aspek yang sama dalam plot yang terlihat ketika membaca seperti itu).

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

EDIT: Ups, lupa untuk mengurangi dari atas, sekarang sudah diperbaiki - masih ada masalah setengah-sel yang belum saya dapatkan sampai ke dasar juga.

Ada beberapa masalah dengan file ini atau dengan GDAL. Saya menggunakan windows 7

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

dan

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

Perhatikan bahwa NoDataValue sama dengan nilai Bmin (-9999), yang ganjil. Yang lebih buruk adalah bahwa GDType adalah UInt16 - Integer 2 byte yang tidak ditandatangani - yang berarti Anda tidak dapat memiliki nilai lebih rendah dari nol. Ini mungkin bug yang diperbaiki di gdal 1.8.0

Masalahnya diilustrasikan ketika Anda melakukannya

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

Saya pikir cara berpuasa untuk memperbaikinya adalah:

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

Masalah ini tampaknya disebabkan oleh masalah mengenali fakta bahwa data dalam format integer 2 byte yang ditandatangani. Ini salah ditafsirkan sebagai format integer 2 byte yang tidak ditandatangani. Karenanya nilai nodata Anda dari -9999 menjadi: 2bytes = 256 * 256 -9999 = 55537

Apa yang saya temukan aneh adalah bahwa nilai min: -9999 dan nilai maks: 5483 adalah sama untuk windows dan mac. Tampaknya dalam kedua kasus tidak ada data yang tidak diidentifikasi dengan benar ketika membangun header, tetapi ketika benar-benar menggunakannya untuk nilai kesalahan terjadi.

solusi:

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

Untuk menggali lebih dalam: Tampaknya raster memanggil rgdal, yang pada gilirannya memanggil gdal itu sendiri. Kemungkinan besar Anda memiliki versi gdal yang berbeda di sistem Anda. Periksa saat memuat rgdal misalnya:

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

Saya baru saja melakukan pemeriksaan cepat di linux: gdal 1.8 memuat file dengan baik, tetapi gdal 1.6 gagal. Jadi sepertinya memang disebabkan oleh gdal.

Meskipun saya tidak yakin dengan kebutuhan Anda, Anda dapat mengonversi. File DEM menjadi file .GRID. Kemudian, arcgis geoprocessor atau R akan secara otomatis mengenali .GRIDs dengan nilai N / A selama manipulasi grid raster.