MODIS icon indicating copy to clipboard operation
MODIS copied to clipboard

Global extent GDAL mosaic problem

Open mstrimas opened this issue 7 years ago • 3 comments

I'm using v1.1.3 with GDAL 2.1.3, I don't have MRT. Trying to produce a global MODIS landcover raster using the following code:

runGdal(product = "MCD12Q1.006", 
        extent = raster::extent(-180, 180, -90, 90),
        begin = "2016.01.01", end = "2016.12.31", 
        SDSstring = "01", job = "global-test")

All tiles are downloaded correctly and everything appears to work fine, but when I look at the raster it's the correct height by only 1 cell wide:

class       : RasterLayer 
dimensions  : 43200, 1, 43200  (nrow, ncol, ncell)
resolution  : 463.3127, 463.3127  (x, y)
extent      : -1e-06, 463.3127, -10007555, 10007555  (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 : /Users/mes335/data/modis/processed/global-test/MCD12Q1.A2016001.LC_Type2.tif 
names       : MCD12Q1.A2016001.LC_Type2 
values      : 0, 255  (min, max)

mstrimas avatar Aug 15 '18 14:08 mstrimas

Hi I am just guessing now, but as you do global coverage, try without extent at all.

MatMatt avatar Aug 15 '18 19:08 MatMatt

With no extent given it opens an interactive map to pick tiles.

mstrimas avatar Aug 15 '18 20:08 mstrimas

This seems to be related to outProj = 'asIn' (ie. target projection is MODIS Sinusoidal). Using outProj = '4326', the output looks almost like expected. Note that I tested this with MOD13A3.006 which requires a little less memory and computation power.

tfs1 = runGdal(product = "MOD13A3.006", 
               extent = raster::extent(-180, 180, -90, 90),
               begin = "2016.01.01", end = "2016.01.01", 
               SDSstring = "01", job = "global-test-ndvi", 
               outProj = "+init=epsg:4326", overwrite = TRUE)

raster(unlist(tfs1))
# class       : RasterLayer 
# dimensions  : 21600, 43210, 933336000  (nrow, ncol, ncell)
# resolution  : 0.008331493, 0.008333333  (x, y)
# extent      : -180, 180.0038, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
# data source : /.../MODIS_ARC/PROCESSED/global-test-ndvi/MOD13A3.A2016001.1_km_monthly_EVI.tif 
# names       : MOD13A3.A2016001.1_km_monthly_EVI 
# values      : -32768, 32767  (min, max)

However, I wonder why xmax of the resulting mosaic extends minimally beyond the dateline. I'll keep investigating and come back to you once I come up with a solution.

fdetsch avatar Sep 11 '18 07:09 fdetsch