Reproject MODIS data using R (results in NAs or no spatial extent)

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.

Reproject MODIS data using R (results in NAs or no spatial extent)

Thread (7 messages)