[R-sig-Geo] Reproject MODIS data using R (results in NAs or no spatial extent)

[Elizabeth Webb]

I am using GLASS albedo data stored here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/AVHRR/> for pre-2000 (AVHRR) data and here<ftp://ftp.glcf.umd.edu/glcf/GLASS/ABD/MODIS/0.05D/> for post-2000 data (MODIS). My end goal is to create a raster stack of each month that contains white sky albedo data from 1982-2015. The problem I have run into is that the MODIS and AVHRR data are in different spatial reference systems and I can't seem to reproject them to be in the same system.

I convert from hdf to tif using R like this:

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)

Here are things I have tried:

1.) Use the MODIS Reprojection Tool<https://lpdaac.usgs.gov/tools/modis_reprojection_tool>. For whatever reason, this tool seems to think the subdatasets of the MODIS .hdf files are only one tile (the upper left most tile, tile 0,0) and not the global dataset. My understanding is that the MODIS data are global (not in tiles?), so I do not know why the MRT is doing this.


2.) Use the raster package in R.

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

This returns a raster with values that are all 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.) Use the gdalUtils package in R:

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

This returns a raster with essentially no spatial extent.

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)

Any ideas on why reprojecting this MODIS data is so difficult?


	[[alternative HTML version deleted]]