MOD09A1 conversion from hdf to tif- sur_relf_qc (band 8) issue?

[ropensci]

MOD09A1 conversion from hdf to tif- sur_relf_qc (band 8) issue?

[lbusett]

Hi,

Thanks for signaling this. We'll have a look ASAP. Did you try by any chance also to use MODIStsp functionality for extraction of bit-fields information, by selecting also the "MODLAND QA bits" output as shown below ?

Do you have the same results ?

Lorenzo

Hi, Yes - to explore the issue myself I selected the top 10 quality indicators to check the results - the output for the QA bits was a tiff with values that are all 0 Rob Rob Critchlow Postdoctoral Research Associate Department of Biology University of York York YO10 5DD UK

…

On 27 April 2017 at 10:03, Lorenzo Busetto ***@***.***> wrote: Hi, Thanks for signaling this. We'll have a look ASAP. Did you try by any chance also to use MODIStsp functionality for extraction of bit-fields information, by selecting also the "MODLAND QA bits" output as shown below ? [image: image] <https://cloud.githubusercontent.com/assets/4480397/25476022/0ac93d10-2b39-11e7-9991-3d6635f8dde5.png> Do you have the same results ? Lorenzo — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#72 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AYDP_EcIMLzf6VYaEgWLql0RDxYZzjFtks5r0FnrgaJpZM4NJ4mu> .

[lbusett]

Are all the QA layers always at 0, or just the one derived from conversion of bits 31-32 ?

No, values are 0 for 5 of the 10 QA outputs: adj_Cor atm_cor b1 b3 QA_bits Rob Critchlow Postdoctoral Research Associate Department of Biology University of York York YO10 5DD UK

…

On 27 April 2017 at 12:53, Lorenzo Busetto ***@***.***> wrote: Are all the QA layers always at 0, or just the one derived from conversion of bits 31-32 ? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#72 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AYDP_MXENl5F54HIOQwtQb0oVePG5vUBks5r0IG0gaJpZM4NJ4mu> .

[lbusett]

ok, thanks for the info.
I'll be able to have a look tomorrow

[lbusett]

Hi Rob,
could you please save your processing options (corresponding to the image for which you sent me the RData) from the GUI and send me the JSON file ?
Lorenzo

[lbusett]

Hi Rob,

I probably managed to replicate your issue, and it seems it could be related to a possible bug in gdalbuildvrt (not the R implementation but the gdal command itself) in "reading" 32 bit unsigned files.

To verify this, could you try to run the following simple script on your machine ? You just have to modify the path to the MOD09A1.A2005361.h21v09.006.2015112231520.hdf image and set the path to a "temporary" folder on your machine.

library(gdalUtils)
library(raster)

# set the path to the MOD09A1.A2005361.h21v09.006.2015112231520.hdf image on your machine

file_in          <- "/home/lb/Temp/buttami/MOD09A1.A2005361.h21v09.006.2015112231520.hdf"

# set this path to a "temporary" folder on your machine
temp_folder      <- "/home/lb/Temp/buttami"

MOD_proj_str     <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
outfile_vrt      <- file.path(temp_folder, "testraster_vrt.vrt")
outrep_file_1    <- file.path(temp_folder, "testraster_vrt_trans.tif")
outrep_file_2    <- file.path(temp_folder, "testraster_norm_trans.tif")
nodata_in        <- "4294967295"
nodata_out       <- "4294967295"
out_format       <- "GTiff"

# build the vrt file
gdalbuildvrt(file_in, outfile_vrt,  sd = 8,srcnodata = nodata_in, 
             vrtnodata = nodata_out, verbose = T, r = "nearest", 
             overwrite = T) 

# convert the vrt to tif
test_translate_vrt <- gdal_translate(outfile_vrt,  
                                   outrep_file_1, 
                               a_srs = MOD_proj_str, of = out_format, 
                               ot = "Uint32", a_nodata = nodata_in, 
                               overwrite = TRUE, output_Raster = T, r = "nearest", 
                               stats = TRUE, strict = TRUE, verbose = T)

# build the tif file directly from the hdf
hdf_in           <- get_subdatasets(file_in)
test_translate_norm <- gdal_translate(hdf_in[8],  
                                   outrep_file_2, r = "nearest", 
                                   a_srs = MOD_proj_str, of = out_format, 
                                   ot = "Uint32", a_nodata = nodata_out, 
                                   overwrite = TRUE, output_Raster = T, 
                                   stats = TRUE, strict = TRUE)

# compare the two files - they should be identical
diff = (getValues(test_translate_vrt) - getValues(test_translate_norm))
summary(diff)
testraster_norm_trans
 Min.   :-63.000      
 1st Qu.:  0.000      
 Median :  0.000      
 Mean   : -1.411      
 3rd Qu.:  0.000      
 Max.   : 55.000      

Since in MODIStsp we create temporary gdal vrt files to allow fast mosaicing of different tiles without creating temporary "physical" files, this problem affects also the final saved tif files .If this problem is confirmed we will have to change strategy. The strange thing is that it seems to affect only a subset of the possible values, as you can see from this histogram of the differences:

[lbusett]

This should be fixed by Commit b05e405

Feel free to re-open if you still encounter problems.

(Currently, the fix is only in the github version. will be included in CRAN on MODIStsp v1.3.3)