Discrepancies between TIFF and HDF files

[sarahmayor]

I am downloading MOD13Q1 products and extracting EVI directly from the HDF files (after conversion from h4 to h5). I was looking into the accompanying Single-band outputs (TIFF files) and realised when I extract EVI from a coordinate from the raster file directly I get a value, whereas from the HDF files I get a -3000 (NA) value at the exact same point (obviously using the same DOY and tile). I have been searching for the source of these TIFF files to know how they were derived but have not succeeded. Maybe you could help me explain the differences? I am using GDAL version: 2.2.3 and MODStsp version: 1.3.9. Thank you very much for your help with this. Sarah

[lbusett]

@sarahmayor

thanks for reporting. It would be helpful if you could save the processing options from the GUI and cut and paste the content of the resulting JSON file.

Lorenzo

[sarahmayor]

{ "sel_prod": "Vegetation Indexes_16Days_250m (M*D13Q1)", "sensor": "Terra", "prod_version": "6", "start_date": "2005-05-20", "end_date": "2005-06-10", "bandsel": [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "indexes_bandsel": [0, 0, 0, 0], "quality_bandsel": [0, 0, 0, 0, 0, 0, 0, 0, 0], "start_x": 17, "end_x": 17, "start_y": 4, "end_y": 4, "user": "jacqueline.o", "password": "xxxxxxxxx", "use_aria": false, "download_server": "http", "download_range": "full", "proj": "Native", "output_proj4": "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs", "out_res_sel": "Native", "out_res": "Native", "full_ext": "Select MODIS Tiles", "resampling": "near", "out_format": "GTiff", "ts_format": "R rasterStack", "compress": "None", "nodata_change": "No", "scale_val": "No", "delete_hdf": "No", "reprocess": "Yes", "bbox": [null, null, null, null], "out_folder": "/media/sarah/BIG1/sarah/MODIS/testing_missing_pixels", "out_folder_mod": "/media/sarah/BIG1/sarah/MODIS/testing_missing_pixels", "MODIStspVersion": "1.3.9", "custom_indexes": [] }

[lbusett]

@sarahmayor

I can not seem to be able to replicate the issue. Original HDF and "processed" tiff seem properly aligned on my output (I attached my output tiff: could you compare with yours? ).

Could you provide me the coordinates at which you are attempting to extract data ?

Lorenzo

MOD13Q1_EVI_2005_161.zip

[sarahmayor]

Maybe it's the way I'm extracting the EVI. I'm looking at coordinate lat 42.634370 long -6.004843 => vert tile 4 horiz tile 17 line 3535.00 samp 2679.00

In R 1) directly from the tiff xy <- data.frame(x=-6.004843,y=42.63437) coords<-SpatialPoints(xy, proj4string=CRS(llcrs)) extract(raster,coords) = 2591 2) from the HDF evi <- h5read(file,name = "/MODIS_Grid_16DAY_250m_500m_VI/Data Fields/250m 16 days EVI") evi [3535,2679] = -3000

[lbusett]

Could you try replicating these commands on your files?

  inr <- raster::raster("/home/lb/tmp/buttami/mstp2/VI_16Days_250m_v6/EVI/MOD13Q1_EVI_2005_161.tif")
  xy <- data.frame(x=-6.004843,y=42.63437)
  coords<-SpatialPoints(xy, proj4string=CRS("+proj=longlat +datum=WGS84"))
  extract(inr,coords)

> 2572 
>Warning message:
> In .local(x, y, ...) : Transforming SpatialPoints to the CRS of the Raster

  sds  <- gdalUtils::get_subdatasets("/home/lb/tmp/buttami/mstp2_hdf/MOD13Q1.A2005161.h17v04.006.2015157052056.hdf")
  inr2 <- rgdal::readGDAL(sds[2])
  inr2 <- raster::raster(inr2)   
  extract(inr2,coords)

[1] 25720000
Warning message:
In .local(x, y, ...) : Transforming SpatialPoints to the CRS of the Raster

As you can see, in my case I get the same result on the TIFF and HDF (though different from yours, and besides a multiplicative factor due to wrongly applied scale and offset).

Lorenzo

[sarahmayor]

Yes, using your method I get exactly the same but we are still projecting the HDF file spatially. The reason I got a different value is because I reprojected the raster to latlong before extracting the coords (which is also a bit strange we should get the same).

However, if we extract the raw EVI from the file without using an object of class SpatialGridDataFrame just using NCDF4 or h5 packages (I tried both) we get NAs for that same pixel, I really don't understand why.

e.g. library(ncdf4) nc <- nc_open("/media/sarah/BIG1/sarah/MODIS/testing_missing_pixels/h5/MOD13Q1.A2005161.h17v04.006.2015157052056.hdf") v1 <- nc$var[["Data Fields/250m 16 days EVI"]] z_all <- ncvar_get(nc, v1) z_all[3535,2679] = NA

Thanks

[lbusett]

Yes, using your method I get exactly the same but we are still projecting the HDF file spatially.

Good. Though I do not get precisely what you mean by "we are still projecting the HDF file spatially"

The reason I got a different value is because I reprojected the raster to latlong before extracting the coords (which is also a bit strange we should get the same).

Not necessarily, unless you are using NN resampling (and NoData are dealt with correctly)

However, if we extract the raw EVI from the file without using an object of class SpatialGridDataFrame just using NCDF4 or h5 packages (I tried both) we get NAs for that same pixel, I really don't understand why.

e.g. library(ncdf4) nc <- nc_open("/media/sarah/BIG1/sarah/MODIS/testing_missing_pixels/h5/MOD13Q1.A2005161.h17v04.006.2015157052056.hdf") v1 <- nc$var[["Data Fields/250m 16 days EVI"]] z_all <- ncvar_get(nc, v1) z_all[3535,2679] = NA

The most probable reason is that we are not looking at the same pixel. Wild guess: whatdo you get from:

z_all[2679,3535]

or 

z_all[2680,3536]

?

[sarahmayor]

z_all[2680,3536] = 25720000

You're absolutely right! I just assumed based on the MODIS tile converter than sample was horizontal and line was vertical. Guess that clears it up then! Thank you very much

[lbusett]

Glad we solved this. I also can never remember which dimension are rows and which are columns... ;-)

Closing now.