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

10

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
żółta czapka
źródło
wersja rastrowa 1.9-12 na oba systemy
yellowcap
Czy możesz zamieścić sessionInfo()swój post?
Roman Luštrik,
Mam różne wartości na raster_1.8-12 (ale identyczne z twoimi na 1.9-12) w winXP.
Roman Luštrik,
Czy działało dobrze z rastrem_1.8-12, czy było inaczej?
yellowcap

Odpowiedzi:

11

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)

wprowadź opis zdjęcia tutaj

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.

mdsumner
źródło
W rzeczywistości możesz połączyć obie metody (ta odpowiedź i odpowiedzi my / roberts): r <- raster('E020N90.DEM')a następnie uruchomić, values(r)<-readBin("E020N90.DEM", "integer", size = 2, signed = TRUE, n = nrows(r) * ncols(r), endian = "big")a następnie values(r)[values(r)==-9999]<-NA
johanvdw 30.09.11
Ha tak, ale to herezja
mdsumner,
6

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')
RobertH
źródło
1
Ta poprawka jest lepsza niż moja, ponieważ punkty danych w Morzu Kaspijskim również są konwertowane (te punkty są również ujemne). Miły!
johanvdw
6

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.

johanvdw
źródło
Załadowano środowisko wykonawcze GDAL: GDAL 1.7.2, wydany 2010/04/23
Roman Luštrik
W systemie Windows moja wersja GDAL jest również tą cytowaną powyżej (1.7.2.), W OSX mam 1.8.0. Ale dlaczego nie mogę odczytać pliku DEM za pomocą 1.7.2.? Czy jest jakieś obejście?
yellowcap
Otrzymałem różne wyniki w różnych wersjach rastra (patrz moje komentarze powyżej), więc nie jestem całkowicie przekonany, że to GDAL per se .
Roman Luštrik,
Czy możesz opisać, jak rgdalznaleźć zaktualizowaną gdalinstalację na Win7? Pobrałem i zainstalowałem najnowsze gdalpliki binarne (zarówno 32, jak i 64). Zostały one zainstalowane w domyślnej lokalizacji, ale rgdalnadal korzystają z wersji 1.7.2, nawet po aktualizacji.
yellowcap
Aktualizacja rgdal nie jest oczywista i będzie wymagać ponownej kompilacji rgdal. Więcej informacji tutaj .
johanvdw
0

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.

Zły w środku
źródło
Używanie innego oprogramowania do pierwszej konwersji pliku jest możliwe, ale nie tak, jak zamierzałem. Chodziło o to, aby używać R tylko do pobierania, czytania i analizowania pliku.
yellowcap
w zasadzie możesz uruchomić gdaltranslate przez R. używając system2.
johanvdw