[R-sig-Geo] Problem merging rasters after conversion and reprojection of .hdf file

Jonathan Greenberg jgrn at illinois.edu
Mon Apr 28 15:12:59 CEST 2014 

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

On Mon, Apr 28, 2014 at 7:18 AM, Justin Michell <jwm302 at gmail.com> wrote:
> 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

More information about the R-sig-Geo mailing list