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

[Alex M]

gdalUtils has a get_subdatasets() which is really helpful here, it's
just a gdalinfo with grep.

Here's what my code looks like for reprojecting, what I found is that
even if the hdf has the proj definition, gdal tools don't seem to use it
right, so I specify the proj string as the s_srs or a_srs.

.modis_warp <- function(output, t_srs="EPSG:4326", ...){
  # Multiple sources say this is the proj definition of the MODIS Sinusoidal
  MODIS_SRS <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs"

  CO <- c("COMPRESS=DEFLATE","PREDICTOR=2","ZLEVEL=9")

  catcher <- tryCatch(
    gdalwarp(output[1]
             , output[2]
             , s_srs=MODIS_SRS
             , t_srs=t_srs
             , of="GTiff"
             #, srcnodata
             #, dstnodata
             , "-multi"
             #, paste0("-wo NUM_THREADS=",ncpu)#Multithread computatation
             , CO
             , output)
  )

  return(catcher)
}

Sure you could transform the raster in "memory" but that will probably
right a temp grd file to the system anyways, might as well just convert
it to a multi-band tiff for ease of use.

Thanks,
Alex

On 11/29/18 11:39, Michael Sumner wrote:
> Ah, never mind - it's the subdataset discovery that's probably not easy
> with rgdal.
> Sorry for the noise.
> 
> Mike.
> 
> On Fri, 30 Nov 2018 at 06:38 Michael Sumner <mdsumner using gmail.com> wrote:
> 
>> 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
>>
>> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
> 
> 	[[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
>