Reproyectar los datos de MODIS con R (resulta en NAs o sin extensión espacial)

Question

Estoy utilizando los datos de albedo de GLASS almacenados aquí para los datos anteriores a 2000 (AVHRR) y aquí para los datos posteriores a 2000 (MODIS). Mi objetivo final es crear una pila de rásters de cada mes que contenga datos de albedo de cielo blanco de 1982 a 2015. El problema con el que me encuentro es que los datos MODIS y AVHRR están en diferentes sistemas de referencia espacial y no puedo reproyectarlos para que estén en el mismo sistema.

Convierto de hdf a tif usando R así:

fileavhrr <- ".../GLASS02B05.V04.A1990161.2018062.hdf"
filemodis<-".../GLASS02B06.V04.A2013169.2017128.hdf"
gdal_translate(get_subdatasets(filemodis)[10], dst_dataset =
        ".../modis.tif") 
gdal_translate(get_subdatasets(fileavhrr)[8], projwin = c(-180,90,180,50), dst_dataset = ".../avhrr.tif") #ideally I'd only like data north of 50 degrees

avhrr<- raster(".../avhrr.tif")

#class       : RasterLayer 
#dimensions  : 800, 7200, 5760000  (nrow, ncol, ncell)
#resolution  : 0.05, 0.05  (x, y)
#extent      : -180, 180, 50, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs 
#values      : -32768, 32767  (min, max)

modis<- raster(".../modis.tif")

#class       : RasterLayer   
#dimensions  : 3600, 7200, 25920000  (nrow, ncol, ncell) 
#resolution  : 154.4376, 308.8751  (x, y)  
#extent   : -20015109, -18903159, 8895604, 10007555  (xmin, xmax, ymin, ymax)  
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
    +b=6371007.181 +units=m +no_defs   
#values      : -32768, 32767  (min, max)

Aquí hay cosas que he probado:

1.) Utilice la herramienta de reproyección MODIS. Por alguna razón, esta herramienta parece pensar que los subconjuntos de datos de los archivos .hdf de MODIS son sólo un mosaico (el mosaico superior izquierdo, el mosaico 0,0) y no el conjunto de datos global. Tengo entendido que los datos de MODIS son globales (¿no en mosaicos?), así que no sé por qué el MRT hace esto.

2.) Utilice el paquete raster en R.

projectedMODIS <- projectRaster(modis,avhrr,method="bilinear")

Esto devuelve una trama con valores que son todos NA:

class       : RasterLayer  
dimensions  : 800, 7200, 5760000  (nrow,> ncol, ncell) 
resolution  : 0.05, 0.05  (x, y) 
extent      : -180, 180,> 50, 90  (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=longlat +ellps=clrk66 +no_defs  
values      : NA, NA  (min, max)

3.) Utilice el paquete gdalUtils en R:

gdalwarp(srcfile=get_subdatasets(filemodis)[10], dstfile= ".../gdalMODIS_avhrr.tif", s_srs = crs(modis), t_srs =crs(avhrr) )

Esto devuelve un raster con una extensión espacial esencialmente nula.

gdalMODISavhrr<-raster(".../gdalMODIS_avhrr.tif")
#class       : RasterLayer 
#dimensions  : 357, 12850, 4587450  (nrow, ncol, ncell)
#resolution  : 0.02801551, 0.02801573  (x, y)
#extent      : -180, 179.9993, 79.99838, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +ellps=clrk66 +no_defs 
#values      : -32768, 32767  (min, max)

¿Alguna idea de por qué es tan difícil reproyectar estos datos MODIS?

Answer 1

Gracias @AndreJ. Este es efectivamente el problema y Mike de ese post me estaba ayudando con este problema en un foro separado ( http://r-sig-geo.2731867.n2.nabble.com/Reproject-MODIS-data-using-R-results-in-NAs-or-no-spatial-extent-td7592500.html ).

El problema es que los metadatos son incorrectos. El producto MODIS GLASS está en la proyección cilíndrica equidistante, pero como los metadatos lo etiquetan incorrectamente como cuadrícula sinusoidal, gdalwarp está leyendo la proyección como sinusodal. Resulta que ambos conjuntos de datos AVHRR y MODIS son conjuntos de datos de rejilla de modelado climático (CMG), y por lo tanto están en equidistante cilíndrico con lat/long geográfico. Los datos MODIS, por alguna razón, utilizan unidades diferentes, pero se puede utilizar a_ullr=c(180, 90, 180, -90) en gdalwarap para que los datos MODIS estén en lat/long. Después de eso, los datos AVHRR y MODIS se apilarán.

Answer 2

Yo intentaría projectedMODIS <- projectRaster(modis, crs = as.character(avhrr@crs)) debería obtener un raster con el mismo crs que el de AVHRR.