Problem merging rasters after conversion and reprojection of .hdf file
Also note you can use ?batch_gdal_translate to do your HDF -> TIF conversions. --j
On Mon, Apr 28, 2014 at 8:12 AM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
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
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