Possible error during unpacking of MOD13Q1 hdf files?

[calvus17]

I have just updated to the new version 1.1.2 of the package and downloaded with runGdal MOD13Q1 hdf files (e.g. MOD13Q1.A2018081.h19v04.006.2018097234738.hdf). The hdf seems to be OK, but the unpacked NDVI layer is strange, the dimension and resolution differs from the images created by the previous version of the package. What could be the problem?
Thanks for helping!

The modis package extracted version of NDVI looks like this:

class : RasterLayer
dimensions : 4800, 7249, 34795200 (nrow, ncol, ncell)
resolution : 231.6418, 231.6564 (x, y)
extent : 933039.2, 2612211, 4447802, 5559753 (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
data source : /home/norbi/Documents/MODIS/Monitoring/Kicsomagolt_raszterek/MODIS/MOD13Q1.A2018081.250m_16_days_NDVI.tif
names : MOD13Q1.A2018081.250m_16_days_NDVI
values : -32768, 32767 (min, max)

The manually extracted version of NDVI band in R:

class : RasterLayer
dimensions : 4800, 4800, 23040000 (nrow, ncol, ncell)
resolution : 231.6564, 231.6564 (x, y)
extent : 1111951, 2223901, 4447802, 5559753 (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
data source : /home/norbi/Documents/MODIS/Monitoring/Letoltott_HDF/MOD13Q1.A2018081.h19v04.006.2018097234738 NDVI .tif
names : MOD13Q1.A2018081.h19v04.006.2018097234738_NDVI_
values : -32768, 32767 (min, max)

[fdetsch]

Awesome, thanks for reporting – and sorry for the long delay.

Dimensions and resolution issues should be fixed as of 72fac2d when working with 'tileH,tileV':

# devtools::install_github("MatMatt/MODIS", ref = "develop")
library(MODIS)

tfs = runGdal("MOD13Q1", begin = "2018081", end = "2018081"
              , tileH = 19, tileV = 4, SDSstring = "100000000000", job = "tmp"
              , overwrite = TRUE, outProj = "asIn")

rst = raster(unlist(tfs))
rst
# class       : RasterLayer 
# dimensions  : 4800, 4800, 23040000  (nrow, ncol, ncell)
# resolution  : 231.6564, 231.6564  (x, y)
# extent      : 1111951, 2223901, 4447802, 5559753  (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 
# data source : C:\...\MODIS_ARC\PROCESSED\tmp\MOD13Q1.A2018081.250m_16_days_NDVI.tif 
# names       : MOD13Q1.A2018081.250m_16_days_NDVI 
# values      : -32768, 32767  (min, max)

Moreover as of 0ab9a66, this also works for spatial subsets defined as Extent or bbox objects from raster and sf, respectively.

library(sf)
library(mapview)

tfs_ext = runGdal("MOD13Q1", begin = "2018081", end = "2018081"
              , extent = st_bbox(franconia), SDSstring = "100000000000", job = "tmp2"
              , overwrite = TRUE, outProj = "asIn")

rst_ext = raster(unlist(tfs_ext))
rst_ext
# class       : RasterLayer 
# dimensions  : 817, 1141, 932197  (nrow, ncol, ncell)
# resolution  : 231.6564, 231.6564  (x, y)
# extent      : 633811.8, 898131.7, 5433268, 5622531  (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 
# data source : C:\...\MODIS_ARC\PROCESSED\tmp2\MOD13Q1.A2018081.250m_16_days_NDVI.tif 
# names       : MOD13Q1.A2018081.250m_16_days_NDVI 
# values      : -32768, 32767  (min, max)

@calvus17 you're working on a Linux OS, right? Could you please test the above code?

[calvus17]

Thank you for your reply!
I tried the above code but unfortunately I got the same results as before. The x coordinate deviate from the y, latter seems to be good. I will try it later on an other computer. I use R studio Version 1.1.442 with R 3.4.4.
Kind regards,
Norbert

tfs = runGdal("MOD13Q1", begin = "2018081", end = "2018081"

•           , tileH = 19, tileV = 4, SDSstring = "100000000000", job = "tmp"
  

•           , overwrite = TRUE, outProj = "asIn")
  

No collection specified, getting the newest for MOD13Q1
########################
outProj = asIn
pixelSize = asIn
resamplingType = near
Output Directory = /home/norbi/Documents/MODIS/Monitoring/Kicsomagolt_raszterek/tmp
########################
Local structure is up-to-date. Using offline information!
Creating output file that is 7249P x 4800L.
Processing input file HDF4_EOS:EOS_GRID:/home/norbi/Documents/MODIS/Monitoring/Letoltott_HDF/MODIS/MOD13Q1.006/2018.03.22/MOD13Q1.A2018081.h19v04.006.2018097234738.hdf:MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI.
0...10...20...30...40...50...60...70...80...90...100 - done.

rst = raster(unlist(tfs))
rst
class : RasterLayer
dimensions : 4800, 7249, 34795200 (nrow, ncol, ncell)
resolution : 231.6418, 231.6564 (x, y)
extent : 933039.2, 2612211, 4447802, 5559753 (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
data source : /home/norbi/Documents/MODIS/Monitoring/Kicsomagolt_raszterek/tmp/MOD13Q1.A2018081.250m_16_days_NDVI.tif
names : MOD13Q1.A2018081.250m_16_days_NDVI
values : -32768, 32767 (min, max)

[fdetsch]

Hm, I cannot reproduce that...

But you did install the MODIS 'develop' version before running the code (see first commented line in my first code chunk), right? In addition, what's your output of sessionInfo().

[calvus17]

OK, I missed that line, it works now perfectly!
Thanks a lot!

[fdetsch]

Great, glad it finally works. Again, thanks for notifying - it's really important that we fixed this!

[calvus17]

Dear Florian, I have just updated the MODIS package to version 1.2.2 and experience the same problem during unpacking the hdf file to geotiff that was already solved two years ago. I use MOD13Q1 and the resulting tiff files have a wrong resolution, not 4800*4800. I am able to get around the problem by extracting the tiff files with my own code, but if you have the chance then please take a look at on this. Thanks! Kind regards, Norber Móricz Florian Detsch <notifications@github.com> ezt írta (időpont: 2018. júl. 11., Sze, 10:02):

…

Awesome, thanks for reporting – and sorry for the long delay. Dimensions and resolution issues should be fixed as of 72fac2d <72fac2d> when working with 'tileH,tileV': # devtools::install_github("MatMatt/MODIS", ref = "develop") library(MODIS) tfs = runGdal("MOD13Q1", begin = "2018081", end = "2018081" , tileH = 19, tileV = 4, SDSstring = "100000000000", job = "tmp" , overwrite = TRUE, outProj = "asIn") rst = raster(unlist(tfs))rst# class : RasterLayer # dimensions : 4800, 4800, 23040000 (nrow, ncol, ncell)# resolution : 231.6564, 231.6564 (x, y)# extent : 1111951, 2223901, 4447802, 5559753 (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 # data source : C:\...\MODIS_ARC\PROCESSED\tmp\MOD13Q1.A2018081.250m_16_days_NDVI.tif # names : MOD13Q1.A2018081.250m_16_days_NDVI # values : -32768, 32767 (min, max) Moreover as of 0ab9a66 <0ab9a66>, this also works for spatial subsets defined as Extent or bbox objects from *raster* and *sf*, respectively. library(sf) library(mapview) tfs_ext = runGdal("MOD13Q1", begin = "2018081", end = "2018081" , extent = st_bbox(franconia), SDSstring = "100000000000", job = "tmp2" , overwrite = TRUE, outProj = "asIn") rst_ext = raster(unlist(tfs_ext))rst_ext# class : RasterLayer # dimensions : 817, 1141, 932197 (nrow, ncol, ncell)# resolution : 231.6564, 231.6564 (x, y)# extent : 633811.8, 898131.7, 5433268, 5622531 (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 # data source : C:\...\MODIS_ARC\PROCESSED\tmp2\MOD13Q1.A2018081.250m_16_days_NDVI.tif # names : MOD13Q1.A2018081.250m_16_days_NDVI # values : -32768, 32767 (min, max) @calvus17 <https://github.com/calvus17> you're working on a Linux OS, right? Could you please test the above code? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#46 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AkiUB5G0Eovx4LuDiJT68SDCottTtWy4ks5uFbEtgaJpZM4TRQJJ> .