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

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Elizabeth Webb
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

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

Michael Sumner-2
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 <[hidden email]> 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
> [hidden email]
> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

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

Michael Sumner-2
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 <[hidden email]> 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 <[hidden email]> 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
>> [hidden email]
>> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

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

Alex Mandel-2
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 <[hidden email]> 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 <[hidden email]> 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
>>> [hidden email]
>>> 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

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

Michael Sumner-2
In reply to this post by Michael Sumner-2
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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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
>>> [hidden email]
>>> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

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

Kay Kilpatrick
FYI The link that Mike provided is for ocean color products the GLASS land products use a different grid…see this link.

https://modis-land.gsfc.nasa.gov/MODLAND_grid.html

K




On Nov 29, 2018, at 3:16 PM, Michael Sumner <[hidden email]<mailto:[hidden email]>> wrote:

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://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Foceancolor.gsfc.nasa.gov%2Fdocs%2Fformat%2Fl3bins%2F&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=r4qRfvtrxz2uNp0QQtKoz0X3eAkOcPTrUqg4fUHv5GI%3D&amp;reserved=0

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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Flpdaac.usgs.gov%2Ftools%2Fmodis_reprojection_tool&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=WnSIC3kBe%2Fc5%2FeqiAlJxVjhm2BPs9h1XC%2FiPFfp4JCg%3D&amp;reserved=0>. 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
[hidden email]
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=NLV%2FZLET%2FHxlxlsl93Laqra4z5A5bDdleHi198Uxh0I%3D&amp;reserved=0

--
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
[hidden email]
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=NLV%2FZLET%2FHxlxlsl93Laqra4z5A5bDdleHi198Uxh0I%3D&amp;reserved=0


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

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

Elizabeth Webb
Thank you Katherine and Mike.   You helped identify the major problem, which is that the metadata is incorrect.  The MODIS GLASS product is a climate modeling dataset and as such it is in cylindrical equidistant projection (the metadata incorrectly label it as the sinusoidal grid).  [I'm writing this just to clarify,  I certainly didn't understand this before your emails earlier today.]  The AVHRR data appear to be in the same projection (cylindrical equidistant).  The next step should just be to convert the MODIS units to lat/long.  However, the AVHRR data does not specify a datum per se (see below).


As a follow up, because I am new to this, when converting the MODIS units to lat/long, what should I use as a datum?  Or, given that both MODIS and AVHRR appear to be in the climate modeling grid format, can we assume their datums already match?


AVHRR Datum info:

GEOGCS[\"Unknown datum based upon the Clarke 1866 ellipsoid\","
 [6] "    DATUM[\"Not specified (based on Clarke 1866 spheroid)\","
 [7] "        SPHEROID[\"Clarke 1866\",6378206.4,294.9786982138982,"
 [8] "            AUTHORITY[\"EPSG\",\"7008\"]]],"
 [9] "    PRIMEM[\"Greenwich\",0],"
[10] "    UNIT[\"degree\",0.0174532925199433]]"
[11] "Origin = (-180.000000000000000,90.000000000000000)"
[12] "Pixel Size = (0.050000000000000,-0.050000000000000)"



*Note ESPG 7008<https://epsg.io/7008-ellipsoid> says the units are in meters.  Other options exist, such as ESPG 4008<https://epsg.io/4008>, unknown datum based on the Clarke 1866 elipsoid exist.


Thanks for the help so far and for any future help that might come my way.


Elizabeth

________________________________
From: Kilpatrick, Katherine A <[hidden email]>
Sent: Thursday, November 29, 2018 3:40 PM
To: Michael Sumner
Cc: Elizabeth Webb; [hidden email]
Subject: Re: [R-sig-Geo] Reproject MODIS data using R (results in NAs or no spatial extent)

FYI The link that Mike provided is for ocean color products the GLASS land products use a different grid�see this link.

https://modis-land.gsfc.nasa.gov/MODLAND_grid.html<https://urldefense.proofpoint.com/v2/url?u=https-3A__modis-2Dland.gsfc.nasa.gov_MODLAND-5Fgrid.html&d=DwMGaQ&c=pZJPUDQ3SB9JplYbifm4nt2lEVG5pWx2KikqINpWlZM&r=KpmI6Vdu2He-hM565ejK5Q&m=6y9mIOW6mOGqZ8k_68J89wbwswkZx5iN0QPOMcOQCIw&s=nZLE2iJ2s4P03oznAf9jhIYLo8_O4nbdqj3YRCf4irs&e=>

K




On Nov 29, 2018, at 3:16 PM, Michael Sumner <[hidden email]<mailto:[hidden email]>> wrote:

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://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Foceancolor.gsfc.nasa.gov%2Fdocs%2Fformat%2Fl3bins%2F&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=r4qRfvtrxz2uNp0QQtKoz0X3eAkOcPTrUqg4fUHv5GI%3D&amp;reserved=0<https://urldefense.proofpoint.com/v2/url?u=https-3A__na01.safelinks.protection.outlook.com_-3Furl-3Dhttps-253A-252F-252Foceancolor.gsfc.nasa.gov-252Fdocs-252Fformat-252Fl3bins-252F-26amp-3Bdata-3D02-257C01-257Ckkilpatrick-2540rsmas.miami.edu-257C08935d04ff104c1b75fd08d65637a4f5-257C2a144b72f23942d48c0e6f0f17c48e33-257C0-257C0-257C636791194328942432-26amp-3Bsdata-3Dr4qRfvtrxz2uNp0QQtKoz0X3eAkOcPTrUqg4fUHv5GI-253D-26amp-3Breserved-3D0&d=DwMGaQ&c=pZJPUDQ3SB9JplYbifm4nt2lEVG5pWx2KikqINpWlZM&r=KpmI6Vdu2He-hM565ejK5Q&m=6y9mIOW6mOGqZ8k_68J89wbwswkZx5iN0QPOMcOQCIw&s=j5HW0q-bCYjevZppbmv_FnjcMCcz3MDRIA3O7k3uGEM&e=>

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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Flpdaac.usgs.gov%2Ftools%2Fmodis_reprojection_tool&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=WnSIC3kBe%2Fc5%2FeqiAlJxVjhm2BPs9h1XC%2FiPFfp4JCg%3D&amp;reserved=0>. 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
[hidden email]
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=NLV%2FZLET%2FHxlxlsl93Laqra4z5A5bDdleHi198Uxh0I%3D&amp;reserved=0

--
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
[hidden email]
https://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&amp;data=02%7C01%7Ckkilpatrick%40rsmas.miami.edu%7C08935d04ff104c1b75fd08d65637a4f5%7C2a144b72f23942d48c0e6f0f17c48e33%7C0%7C0%7C636791194328942432&amp;sdata=NLV%2FZLET%2FHxlxlsl93Laqra4z5A5bDdleHi198Uxh0I%3D&amp;reserved=0


        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo