Problem merging rasters after conversion and reprojection of .hdf file

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20140428/abd0a8db/attachment.pl>
Justin:

Couple of comments:

1) In general, I find it easier to mosaic FIRST and then reproject
SECOND -- when you reproject first, mosaic second, you can get weird
issues particularly along the seams.
2) Since you are already using gdalUtils, have you tried
?mosaic_rasters (which will be faster than merge() in all likelihood,
since it uses the base GDAL mosaicking capabilities)?  Make sure
separate=FALSE for this particular application.

--j
Dear Group

I am dealing with the task of preprocessing and re-projecting MODIS gridded tiles in .hdf format from Sinusoidal projection, in R, to a geographic projection. I would like to align all my other climate layers (from worldclim.org) which are at 1km resolution to my NDVI layers. I am using the gdal_utils package in R, and only one subset I?m interested in (monthly NDVI). I am working on Mac OS and grabbed gdal from here: http://www.kyngchaos.com/software/frameworks .

To do this I use the following code (previous list entries relating to gdalUtils have assisted me up till this point):

# download all available monthly images for the year 2013 - MODIS PRODUCT (MOD13A3, Terra, Vegetation Indices, Tile, 1000m, Monthly)
ModisDownload(x="MOD13A3",h=c(20,21,22),v=c(8,9),dates=c('2013.01.01','2013.12.31'),mosaic=F,proj=F)

# extract the NDVI band (sd_index=1) and converts to .tiff for 1 tile for January
gdal_translate("MOD13A3.A2013001.h20v08.005.2013042085256.hdf", "Jan_1.tiff", sd_index=1)

# converts projection to Geog WGS84
gdalwarp(?Jan_1.tiff","MODIS_NDVI/NDVIJan_1.tiff", s_srs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs",t_srs="+proj=longlat +datum=WGS84 +no_defs",srcnodata=-3000,dstnodata=-3000)

 I have 6 tiles per month (covering Kenya). For each month I?d like to merge all 6 tiles using the raster package (to cover Kenya). For January I perform the following operation using the merge function in the raster package:

NDVI1  <- (merge(raster("MODIS_NDVI/NDVIJan_1.tiff"),raster("MODIS_NDVI/NDVIJan_2.tiff"),raster("MODIS_NDVI/NDVIJan_3.tiff"), raster("MODIS_NDVI/NDVIJan_4.tiff"), raster("MODIS_NDVI/NDVIJan_5.tiff"), raster("MODIS_NDVI/NDVIJan_6.tiff"), tolerance=0.5))*0.0001
NDVI1 <- crop(NDVI1, extent(kenya))

Without including the tolerance argument I get an error:

Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE,  :
  different origin

 But after this merge a lot of NA?s are introduced and the plot doesn?t look right (as it did before for each tile).

The resolution is not exactly the same. Here are a comparison of the layers (before merging NDVI tiles) and my SpatialPolygonsDataFrame called kenya, as well as my r sessionInfo() printed below:

n1 <- raster("MODIS_NDVI/NDVIJan_1.tiff")
n1
class       : RasterLayer
dimensions  : 1219, 1275, 1554225  (nrow, ncol, ncell)
resolution  : 0.008205785, 0.008205785  (x, y)
extent      : 20, 30.46238, -0.002852273, 10  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /Users/justinmichell/MODIS_NDVI/NDVIJan_1.tiff
names       : NDVIJan_1
values      : -32768, 32767  (min, max)

meanTemp1
class       : RasterLayer
dimensions  : 1115, 960, 1070400  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 33.90833, 41.90833, -4.666667, 4.625  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
data source : in memory
names       : layer
values      : -4.4, 30  (min, max)

kenya
class       : SpatialPolygonsDataFrame
features    : 1
extent      : 33.90722, 41.90517, -4.669618, 4.622499  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
variables   : 11
names       : FIPS, ISO2, ISO3,  UN,  NAME,  AREA,  POP2005, REGION, SUBREGION,    LON,  LAT
min values  :   KE,   KE,  KEN, 404, Kenya, 56914, 35598952,      2,        14, 37.858, 0.53
max values  :   KE,   KE,  KEN, 404, Kenya, 56914, 35598952,      2,        14, 37.858, 0.53

sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] gdalUtils_0.2.0 rgdal_0.8-14    RCurl_1.95-4.1  bitops_1.0-6    raster_2.2-12   sp_1.0-14

loaded via a namespace (and not attached):
[1] codetools_0.2-8      foreach_1.4.1        grid_3.0.2           iterators_1.0.6      lattice_0.20-24      scatterplot3d_0.3-34

Thanks very much,
Justin Michell

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Jonathan A. Greenberg, PhD
Assistant Professor
Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
259 Computing Applications Building, MC-150
605 East Springfield Avenue
Champaign, IL  61820-6371
Phone: 217-300-1924
http://www.geog.illinois.edu/~jgrn/
AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
Also note you can use ?batch_gdal_translate to do your HDF -> TIF conversions.

--j
Justin: