To answer the "Any ideas on why reprojecting this MODIS data is so difficult?" - it's that the sinusoidal projection is actually pretty exotic, it matches the daily aggregation of L2 geophysical variables from individual swaths (multiple per day) into L3 bins (aggregated daily, then into longer time periods with bin identity maintained). I personally think it's a bad choice for a regular gridded product, but it's a way of "regularizing" the L3 bins, which are quite abstract and hard to use. The L3 bins are a ragged array of bins at regular latitude intervals, a diminishing number of bins as absolute latitude increases towards each pole. Described here:
https://oceancolor.gsfc.nasa.gov/docs/format/l3bins/ My answer to the "how to reproject the sinusoidal grid to regular longlat" is - you should not do that. The sinusoidal grid is faithful to the statistical properties of the L3 bins, though not to the bin geometries themselves. One is meant to query the sinusoidal grid for the data needed, it's really inappropriate to stretch those data out by resampling/ reprojection. (I don't know if the sinusoidal-regular-grid creators have written down this rationale). Generally one can reproject vector points, polygons, lines into the sinusoidal projection and use extractions in that coordinate system, and that's what I'd recommend. (Using the L3 bins is also possible in R but it's a long story about how to do that). Just a question: do you download the auxialiary .hdf.xml files that accompany the .hdf files? They perhaps contain extra information that GDAL can interpret. Using some inside knowledge, delve into subdatasets with stars/sf: (sds <- sf::gdal_subdatasets("GLASS02B06.V04.A2013169.2017128.hdf")) x <- stars::read_stars(sds[[1]]) x stars object with 2 dimensions and 1 attribute attribute(s), summary of first 1e+05 cells: HDF4_EOS:EOS_GRID:"GLASS02B06.V04.A2013169.2017128.hdf":GLASS02B06:ABD_BSA_VIS Min. : NA 1st Qu.: NA Median : NA Mean :NaN 3rd Qu.: NA Max. : NA NA's :1e+05 dimension(s): from to offset delta refsys point values x 1 7200 -20015109 154.438 +proj=sinu +lon_0=0 +x_0=... NA NULL [x] y 1 3600 10007555 -308.875 +proj=sinu +lon_0=0 +x_0=... NA NULL [y] plot(x) sf::st_bbox(x) xmin ymin xmax ymax -20015109 8895604 -18903159 10007555 And OMG! To my disgust we see these data are not in sinusoidal projection at all, but in Equidistant Cylindrical (the -20015109 is approximately the distance from the prime meridian to the anti meridian along the equator (i.e about pi * 6378137m). So, I don't trust these in-situ metadata at all. I see this stuff gets messed up fairly frequently, perhaps it's a GDAL problem interpreting the metadata from the file, or perhaps it's actually mis-applied in the hdf - finding out requires careful examination of the metadata, as does choice of datum (as follows). Using stars to read and raster to "fix it" (raster just because I'm more sure of my steps that way): r <- setExtent(stars:::st_as_raster(x), extent(-180, 180, -90, 90)) projection(r) <- "+proj=longlat +datum=WGS84" plot(r) maps::map(add = T) plot(extent(r), add = TRUE) ## convince yourself the extent is correct with things like: abline(h = c(-90, 90), col = "red") That use of "datum=WGS84" may well be incorrect, but compared to trying to project data from a completely mis applied projection it's an improvement. HTH Cheers, Mike. On Fri, 30 Nov 2018 at 06:39 Michael Sumner <mdsum...@gmail.com> 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 <mdsum...@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 <we...@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@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 > > -- 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@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo