Output raster files shifted by a couple of pixels?

[dk-forestry]

I downloaded and processed some scenes from MOD13Q1. While visually exploring the results in ArcGIS, I noticed that MODIStsp outputs seemed to be slightly shifted to the to east and to the south compared to ArcGIS base maps.

To further explore this issue, I imported the original MODIS hdf in ArcGIS. Visually, the original data seems to be aligned with the features of the ArcGIS base maps.

Here's a direct comparison between the original hdf and MODIStsp outputs:

Original MODIS hdf:

The same scene with the MODIStsp output:

It seems that all pixels are shifted two positions to the East and one position to the South.

Is this a problem with ArcGIS, MODIStsp or my workflow? In MODIStsp, I used native projection and resolution to produce GTiff outputs with no compression.

sessionInfo() R version 3.6.1 (2019-07-05) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 17134) MODIStsp_1.3.9

[lbusett]

Hi. This sould definitely not happen. I'll have a look. Could you please save your options file from the GUI and cut and paste here the contents of the JSON file?

Lorenzo

[lbusett]

Another thing: Did you process a full tile, or a subset area?

[dk-forestry]

Thanks for the fast reply!

I processed a custom area, with the extent loaded from a shapefile. The custom area spans over 4 tiles (but none of the tiles is covered entirely by the area).

Here is the content of the JSON file:

{ "sel_prod": "Vegetation Indexes_16Days_250m (M*D13Q1)", "sensor": "Terra", "prod_version": "6", "start_date": "2019-01-01", "end_date": "2019-02-01", "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": 21, "end_x": 22, "start_y": 10, "end_y": 11, "user": "[...]", "password": "[...]", "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": "Define Custom Area", "resampling": "near", "out_format": "GTiff", "ts_format": [], "compress": "None", "nodata_change": "No", "scale_val": "No", "delete_hdf": "No", "reprocess": "No", "bbox": [4330656.1, -2847288.3, 5492217.7, -1328641.9], "out_folder": "[...]", "out_folder_mod": "[...]", "MODIStspVersion": "1.3.9", "custom_indexes": [] }

[dk-forestry]

I just did a test run with the option 'Select MODIS tiles' (and not 'Define custom area'). Now, the output of MODIStsp is properly aligned to the original hdf files. This is the case for both 1 tile and multiple tiles. I guess this means that the issue is limited to custom areas?

[lbusett]

Ok. I'll have a look ASAP. I already checked on Full Tiles, and I do not see misalignment. Let's see.

[lbusett]

@dk-forestry

I cannot seem to be able to replicate the issue.

• Result from HDF of the 001 Image in QGIS ( i had to convert it to GTiff beforehanf, though, because QGIS3 does not seem to be able to open HDF4)

• Result from MODIStsp:

Could you also try to convert an HDF to Geotiff using gdal_translate BEFORE looking in ARCGIS? You can use something like:

gdl_translate HDF4_EOS:EOS_GRID:"/home/lb/tmpmodis/MOD13Q1.A2019001.h21v10.006.2019024152717.hdf":MODIS_Grid_16DAY_250m_500m_VI:"250m 16 days EVI" /home/lb/tmpmodis/pippo.tif

(you'll just have to update the paths).

Lorenzo

• Result from HDF of the 001 Image

[dk-forestry]

Thanks for looking into this! As you suggested, I converted from HDF to Geotiff using gdal_translate, but it made no difference.

I tried a couple of things. First, I removed ArcGIS from the equation by visualizing the issue in R like this:

library(raster)

raster_from_hdf <- raster('pippo.tif')  # pippo.tif was produced using gdal_translate
raster_from_modistsp <- raster('raster_from_modistsp.tif')

crop_extent <- extent(4425250, 4427000, -2702000, -2700500)

raster_from_hdf <- crop(raster_from_hdf, crop_extent)
raster_from_modistsp <- crop(raster_from_modistsp, crop_extent)

modis_stack <- stack(raster_from_hdf, raster_from_modistsp)

png(filename = 'image.png', width=500,height=300)
plot(modis_stack)
dev.off()

The crop_extent shows an area of interest that illustrates the issue. It lays in the tile MOD13Q1.A2019001.h21v11.006.2019024152558.

Then I used MODIStsp() to process the data. Here's the result using a custom area that corresponds to crop_extent (i.e. "bbox": [4425250, -2702000, 4427000, -2700500]). This results in a correct output:

Next, I used MODIStsp() with a custom area that spans over two tiles ("bbox": [4425250, -2702000, 4460000, -2700500]). This produced an output where the pixels are shifted two positions to the east:

Finally, I used a custom area that spans over 4 tiles ("bbox": [4330656.1, -2847288.3, 5492217.7, -1328641.9]). In the result, pixels are shifted two to the east and one to the south:

[lbusett]

This is very strange. My test yesterday was done on "your" subset (therefore, mosaicing 4 tiles and making a subset.

If I send you my MODIStsp result could you compare it with yours?

[dk-forestry]

Maybe there's some issue with my local setup. I'll try to look into that.

If I send you my MODIStsp result could you compare it with yours?

Sure! :)

[lbusett]

Ok. Here it is.

MOD13Q1_EVI_2019_001.zip

These results were obtained yesterday based on your json (all Madagascar, therefore - cutting on 4330656.1,5492217.7,-2847288.3,-1328641.9 )

I already checked using the code you provided and starting from the file I sent you. This is what I got:

raster_from_hdf <- raster("/home/lb/tmp/buttami/mstp1_hdf/pippo2.tif")  # pippo.2tif was produced using gdal_translate
raster_from_modistsp <- raster("/home/lb/tmp/buttami/mstp1/VI_16Days_250m_v6/EVI/MOD13Q1_EVI_2019_001.tif")

crop_extent <- extent(4425250, 4427000, -2702000, -2700500)

raster_from_hdf      <- crop(raster_from_hdf, crop_extent)
raster_from_modistsp <- crop(raster_from_modistsp, crop_extent)

modis_stack <- stack(raster_from_hdf, raster_from_modistsp)

plot(modis_stack)

, so, it appears to be ok.

Could you also share:

1) your session info (devtools::session_info()) 2) the results of

gdal_setInstallation() 
getOption("gdalUtils_gdalPath")[[1]]

3) the results of gdalwarp --version from a terminal

[dk-forestry]

The problem was indeed an old gdal version.

On Windows, rgdal comes currently with GDAL 2.2.3. I somehow (mistakenly) assumed that all R packages would use this gdal version, but this is not the case. Instead, gdal_setInstallation() used GDAL 2.1.0 from an old OSGeo4W in my PATH. After installing an up-to-date version of OSGeo4W, gdal_setInstallation() now uses GDAL 2.4.1 and MODIStsp produces beautiful and correctly aligned outputs.

Thanks again for helping to figure out this weird issue!

[lbusett]

Glad we solved this, though it is weird and problemetic that old gdal versions give that problem. I'll probably need to introduce a check on GDAL version.