The MODIS grid seems to be shifted leading to wrong results returned by the getTile()

[zozlak]

I was searching for MODIS tiles intersecting with my regions of interest

bbox = matrix(c(39.995535714, 40.995535714, -0.995535714, 0.004464286), nrow = 2, ncol = 2, byrow = TRUE)
getTile(bbox)

The result should be h21v08, h21v09, h22v08 & h22v09 but only last two are returned.

While investigating the getTile() function I exported the MODIS grid used by the package (the sr variable in https://github.com/MatMatt/MODIS/blob/master/R/getTile.R#L277) and compared it in the QGIS with a few MODIS HDFs (MOD09A1.006 from 2018-11-17 tiles h21v08, h21v09, h22v08 & h22v09). The result suggests the grid used by the package is a little shifted comparing to actual tiles location (white box - my region of interest, red lines - package grid, colorful tiles - actual MODIS HDFs coverage):

Is there a reason for that? How can I properly get the MODIS tiles my region intersects with?

[fdetsch]

I agree that this needs further investigation. The boundaries in MODIS:::sr and actual G-Polygon coordinates from the granules' metadata differ quite a bit. Until this is fixed, one possible workaround would be to extend your target spatial extent, for example by 0.5 degree. This makes sure that all adjacent tiles are downloaded as well. Finally, just crop() the downloaded and extracted images using your original extent as input.

MODISoptions(outProj = "+init=epsg:4326")

## extend target extent by .5 deg
bbox = extent(c(39.995535714, 40.995535714, -0.995535714, 0.004464286))
bbox_xtd = extend(bbox, .5) 

## download and extract images
tfs = runGdal(product = "MOD09A1", collection = NULL
              , begin = as.Date("2018-11-17"), end = as.Date("2018-11-17")
              , extent = bbox_xtd
              , SDSstring = "1100000000000", job = "issue#59"
              , overwrite = TRUE)
tfs = unlist(tfs, use.names = FALSE)

## import images and clip to target extent
rst = stack(tfs)
crp = crop(rst, bbox, snap = "out")

[zozlak]

@fdetsch thanks for the advice.

[zozlak]

I looked at different grids:

• (blue) the MODIS:::sr grid mentioned above (one used by the MODIS::getTile())
• (red) the official MODIS bounding boxes (https://modis-land.gsfc.nasa.gov/pdf/sn_bound_10deg.txt)
• (green) the official MODIS g-rings (https://modis-land.gsfc.nasa.gov/pdf/sn_gring_10deg.txt)
• (violet) a grid constructed by combining neighbour bounding boxes coordinates

The comparison for a sample tile (h24v07) looks as follows (the projection is +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs):

The bbox grid, for obvious reasons, has the worse overall match but it matches lower left and upper right corner really well. The MODIS:::sr grid is shifted in both directions. The g-ring grid matches upper and right bounds well but has problems with left and lower bounds. The best fit is provided by the grid constructed by combining neighbour bounding boxes (but it doesn't work for tiles intersecting the end of the world).

An R code used to obtain the grids:

library(dplyr)

toGeoJson = function(df, file) {
  df = df %>%
    dplyr::group_by(iv, ih) %>%
    dplyr::mutate(
      polygon = list(sp::Polygons(polygon, paste(iv, ih, sep = '_')))
    )
  spdf = sp::SpatialPolygonsDataFrame(
    sp::SpatialPolygons(df$polygon),
    data.frame(iv = df$iv, ih = df$ih, row.names = paste(df$iv, df$ih, sep = '_'))
  )
  rgdal::writeOGR(spdf, file, '', driver = 'GeoJSON')
}

# grid used by MODIS::getTile()
rgdal::writeOGR(MODIS:::sr, 'modis_sr.geojson', '', driver = 'GeoJSON')

# https://modis-land.gsfc.nasa.gov/pdf/sn_bound_10deg.txt
system("curl -k https://modis-land.gsfc.nasa.gov/pdf/sn_bound_10deg.txt | tr '\15' '\12' | head -n -2 > modis_bound.txt")
x1 = read.fwf('modis_bound.txt', c(3, 4, 11, 11, 10, 10), skip = 7, col.names = c('iv', 'ih', 'lon_min', 'lon_max', 'lat_min', 'lat_max')) %>%
  dplyr::as.tbl() %>%
  dplyr::filter(lon_min > -200) %>%
  dplyr::group_by(iv, ih) %>%
  dplyr::mutate(
    polygon = list(sp::Polygon(cbind(c(lon_min, lon_min, lon_max, lon_max, lon_min), c(lat_min, lat_max, lat_max, lat_min, lat_min))))
  )
toGeoJson(x1, 'modis_bound.geojson')

# https://modis-land.gsfc.nasa.gov/pdf/sn_gring_10deg.txt
system("curl -k https://modis-land.gsfc.nasa.gov/pdf/sn_gring_10deg.txt | tr '\15' '\12' | head -n -2 > modis_gring.txt")
x2 = read.fwf('modis_gring.txt', c(3, 4, 11, 10, 11, 10, 11, 10, 11, 10), skip = 7, col.names = c('iv', 'ih', 'll_x', 'll_y', 'ul_x', 'ul_y', 'ur_x', 'ur_y', 'lr_x', 'lr_y'), stringsAsFactors = FALSE) %>%
  dplyr::as.tbl() %>%
  dplyr::filter(ll_x > -200) %>%
  dplyr::group_by(iv, ih) %>%
  dplyr::mutate(
    polygon = list(sp::Polygon(cbind(c(ll_x, lr_x, ur_x, ul_x, ll_x), c(ll_y, lr_y, ur_y, ul_y, ll_y))))
  )
toGeoJson(x2, 'modis_gring.geojson')

# variation of modis bound
x3 = x1 %>%
  dplyr::group_by(iv) %>%
  dplyr::arrange(iv, ih) %>%
  dplyr::mutate(
    ll_x = lon_min,
    lr_x = lead(lon_min),
    ul_x = lag(lon_max),
    ur_x = lon_max
  ) %>%
  dplyr::filter(!is.na(lr_x) & !is.na(ul_x)) %>%
  dplyr::group_by(iv, ih) %>%
  dplyr::mutate(
    polygon = list(sp::Polygon(cbind(c(ll_x, lr_x, ur_x, ul_x, ll_x), c(lat_min, lat_min, lat_max, lat_max, lat_min))))
  )
toGeoJson(x3, 'modis_bound2.geojson')