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

[Michael Sumner]

Fwiw there shouldn't be any need to convert from hdf to tif - could you
please try this?

x <- readGDAL( ".../GLASS02B05.V04.A1990161.2018062.hdf")

If that works it at least removes some steps from your process.

Cheers, Mike.

On Fri, Nov 30, 2018, 05:03 Elizabeth Webb <webbe using ufl.edu> wrote:

> 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]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]