[R-sig-Geo] Reproject MODIS data using R (results in NAs or no spatial extent)
Michael Sumner
md@umner @ending from gm@il@com
Thu Nov 29 21:16:41 CET 2018
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 <mdsumner using 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 <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
>
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list