1

I am using GLASS albedo data stored here for pre-2000 (AVHRR) data and here 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. 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?

1

Thank you @AndreJ. This is indeed the problem and Mike from that post was helping me with this problem in a separate forum (http://r-sig-geo.2731867.n2.nabble.com/Reproject-MODIS-data-using-R-results-in-NAs-or-no-spatial-extent-td7592500.html).

The problem is that the metadata is incorrect. The MODIS GLASS product is in the cylindrical equidistant projection but because the metadata incorrectly label it as the sinusoidal grid, gdalwarp is reading the projection as sinusodal. It turns out that both AVHRR and MODIS data sets are climate modeling grid (CMG) datasets, and are thus in cylindrical equidistant with geographic lat/long. The MODIS data for whatever reason use different units, but one can use a_ullr=c(180, 90, 180, -90) in gdalwarap to make the MODIS data in lat/long. After that, AVHRR and MODIS data will stack.

  • The grids are in a Geographic CRS (using degrees), not cylindrical equidistant (which uses meters). The false coordinates from MODIS might be a follow-up of specifying the wrong sinu CRS. The XDIM and YDIM values are not stored as datafields, but calculated from the CRS and cell position. – AndreJ Dec 2 '18 at 11:30
0

I would try projectedMODIS <- projectRaster(modis, crs = as.character(avhrr@crs)); you should obtain a raster with the same crs as the AVHRR one.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.