Estoy tratando de leer un archivo ráster en formato .DEM en Windows usando el paquete 'ráster' en R.
Tengo problemas con los valores de NA al cargar los datos en R en Windows 7, pero no tengo el problema en una Mac con OSX Lion. En Windows, los valores de NA no parecen leerse correctamente. La pregunta es ¿por qué sucede esto?
El archivo ráster utilizado se descargó de USGS con el siguiente código R:
download.file('http://edcftp.cr.usgs.gov/pub/data/gtopo30/global/e020n90.tar.gz', 'e020n90.tar.gz')
untar('e020n90.tar.gz')
Luego leí el ráster en R usando el paquete 'ráster'. En OSX Lion y R64 versión 2.13.1, se reconocen los valores de 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
Pero en Windows 7 (64 bits, la misma versión R) convierte los valores de celda que deberían ser NA en números:
> 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
¿Por qué no hay valores de NA en el ráster cuando lo leo en Windows? ¿Cómo podría solucionarlo? Supongo que tiene que ver con la forma en que se almacenan los números, muchos de los valores de NA se convierten a 55540.
Información de Windows (después de cargar el ráster):
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
Información de OSX (después de cargar el ráster):
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
sessionInfo()
en tu publicación?Respuestas:
Una solución alternativa es elegir los datos sin procesar, ya que este es un formato de archivo muy simple.
No es para todos, pero puede ser esclarecedor ver lo que está sucediendo.
En este punto, puede probar las diferentes opciones para el signo entero y la endianidad directamente, y al leer de esta manera, logramos lo que Robert hace con la
> 32767
transformación después de leer el archivo.Finalmente, configure la proyección tal como la lee el ráster (y esto daría la misma relación de aspecto en el gráfico que se ve cuando se lee de esa manera).
EDITAR: Whoops, se había olvidado de restar desde la parte superior, ahora solucionado: todavía hay un problema de media celda que no he llegado al fondo también.
fuente
r <- raster('E020N90.DEM')
y luego ejecutarvalues(r)<-readBin("E020N90.DEM", "integer", size = 2, signed = TRUE, n = nrows(r) * ncols(r), endian = "big")
y luegovalues(r)[values(r)==-9999]<-NA
Hay algunos problemas con este archivo o con GDAL. Estoy usando windows 7
y
Tenga en cuenta que NoDataValue es el mismo que el valor de Bmin (-9999), que es impar. Lo que es peor es que GDType es UInt16: enteros de 2 bytes sin signo, lo que significa que no puede tener valores inferiores a cero. Este es probablemente un error que se corrigió en gdal 1.8.0
El problema se ilustra cuando lo haces
Creo que la forma más rápida de solucionar esto es:
fuente
El problema parece ser causado por un problema que reconoce el hecho de que los datos están en formato entero de 2 bytes con signo. Se interpreta erróneamente como un formato entero de 2 bytes sin signo. Por lo tanto, su valor de nodata de -9999 se convierte en: 2bytes = 256 * 256 -9999 = 55537
Lo que me parece extraño es que el valor mínimo: -9999 y el valor máximo: 5483 son los mismos para Windows y Mac. Parece que en ambos casos no se identificaron correctamente los datos al compilar los encabezados, pero al usarlos realmente para los valores se produjo un error.
solución alterna:
Para profundizar: parece que la trama llama a rgdal, que a su vez llama a gdal. Lo más probable es que tenga una versión diferente de gdal en su sistema. Compruebe al cargar rgdal, por ejemplo:
Loaded GDAL runtime: GDAL 1.8.0, released 2011/01/12
Acabo de hacer una comprobación rápida en Linux: gdal 1.8 carga bien el archivo, pero falla gdal 1.6. Por lo tanto, parece ser causado por gdal.
fuente
rgdal
encontrar unagdal
instalación actualizada en Win7? Descargué e instalé losgdal
binarios más recientes (32 y 64). Estos se instalaron en la ubicación predeterminada perorgdal
aún usan 1.7.2, incluso después de la actualización.Aunque no estoy seguro acerca de su requisito, puede convertir. Archivos DEM en archivos .GRID. Luego, el geoprocesador arcgis o R reconocerá automáticamente .GRIDs con valores N / A durante la manipulación de la trama de cuadrícula.
fuente
system2
.