Hi Jonathan

Thanks, very helpful!

I have followed your suggestions (except I already had a loop going so I didn’t use batch_gdal_translate). But after this process when I I try to match these NDVI values (for Jan as shown below) to the coordinates of my dataset, I get mostly NA’s. I am not sure if this makes sense, although the reprojected map when I plot it looks fine I think.

My coordinates in my spatial dataset represent 464 unique sample sites scattered across Kenya with a count of the number people examined and positive at each site, and downloaded from the malaria atlas project (http://www.map.ox.ac.uk/browse-resources/). Any thoughts would be appreciated!   

# Jan tiles (4 tiles cover Kenya)

> for (i in 1 : 4){
+   # extract the NDVI band (sd_index=1) and converts to .tiff 
+   gdal_translate(out.files[i], paste0("Jan", "_", i, ".tiff"), sd_index=1) 
+ }

> mosaic_rasters(gdalfile=c("Jan_1.tiff", "Jan_2.tiff", "Jan_3.tiff", "Jan_4.tiff"),dst_dataset="NDVI1_mosaic.tiff", separate=F, verbose=TRUE)

> gdalwarp("NDVI1_mosaic.tiff","MODIS_NDVI/NDVIJan.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)   

> NDVI1 <- raster("MODIS_NDVI/NDVIJan.tiff")*0.0001
> NDVI1 <- crop(NDVI1, extent(kenya))

> NDVI1
class       : RasterLayer 
dimensions  : 1107, 953, 1054971  (nrow, ncol, ncell)
resolution  : 0.008397857, 0.008397857  (x, y)
extent      : 33.905, 41.90816, -4.671056, 4.625372  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : NDVIJan 
values      : -0.2, 0.9031  (min, max)

> summary(NDVI1)
            NDVIJan
Min.        -0.2000
1st Qu.      0.1813
Median       0.2892
3rd Qu.      0.4550
Max.         0.9031
NA's    652208.0000

> extract(NDVI1, sitesPolyPoints@coords)

 extract(NDVI1, sitesPolyPoints@coords)
  [1] 0.5094     NA 0.6557     NA     NA 0.6575 0.6493 0.6469     NA 0.6607     NA     NA     NA 0.6814     NA     NA 0.6620     NA     NA     NA     NA 0.6749 0.6887     NA 0.6880 0.5959     NA
 [28]     NA     NA 0.6570 0.6212     NA     NA     NA     NA     NA 0.6647     NA     NA 0.6296     NA     NA 0.5831     NA     NA     NA 0.6147     NA 0.6167     NA 0.6258 0.6615     NA 0.6613
 [55]     NA 0.6626     NA     NA 0.6406     NA 0.6451 0.6718     NA 0.6641     NA     NA     NA     NA 0.6479     NA 0.6506     NA 0.6478     NA 0.6629     NA     NA     NA 0.6519     NA 0.6699
 [82]     NA     NA     NA     NA 0.5996     NA 0.6694     NA 0.6855     NA 0.6625     NA     NA     NA     NA     NA     NA     NA 0.6389     NA     NA     NA     NA     NA     NA     NA     NA
[109]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA 0.6203     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
[136] 0.6159     NA     NA 0.6538     NA 0.6613 0.6731     NA     NA     NA 0.6781     NA 0.6090     NA     NA 0.6054     NA     NA     NA 0.6810     NA     NA     NA     NA     NA     NA     NA
[163] 0.6486     NA     NA 0.6328     NA     NA 0.6681     NA 0.6681 0.6766     NA     NA 0.6910     NA     NA     NA     NA     NA     NA     NA     NA     NA 0.6706     NA     NA     NA     NA
[190]     NA 0.6839     NA 0.6632     NA     NA     NA     NA     NA     NA 0.6716     NA     NA     NA     NA     NA     NA 0.6282     NA     NA     NA     NA 0.6607     NA     NA     NA     NA
[217]     NA     NA     NA     NA     NA 0.2105     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA 0.7046     NA 0.7367     NA 0.6968     NA     NA
[244]     NA 0.5180 0.3300 0.7082 0.5610 0.3725     NA     NA     NA     NA     NA     NA 0.2636 0.7032 0.6201 0.6996     NA     NA 0.3071 0.3292 0.4628     NA     NA     NA     NA     NA     NA
[271]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
[298]     NA     NA     NA     NA 0.7367     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
[325]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA 0.7073     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
[352]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
[379]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
[406]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
[433]     NA     NA     NA     NA     NA     NA     NA     NA     NA 0.3652     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA

> head(sitesPolyPoints@coords)
     coords.x1 coords.x2
[1,]   33.9940    0.1580
[2,]   34.1164   -0.0115
[3,]   34.1175    0.3770
[4,]   34.1277   -0.4222
[5,]   34.1318   -0.0035
[6,]   34.1350    0.4016

Cheers
Justin


On Apr 29, 2014, at 2:56 PM, Justin Michell <jwm302@gmail.com> wrote:

> Hi Jonathan
> 
> Thanks for you help here. I get the following error when I tried the following as per you advice (and the example for mosaic_rasters). Do I not have the right GDALDriver?:
> 
>> gdal_setInstallation()
>> valid_install <- !is.null(getOption("gdalUtils_gdalPath"))
>> valid_install
> [1] TRUE
>> if(require(raster) && require(rgdal) && valid_install)
> + {
> + NDVI_1 <- mosaic_rasters(gdalfile=c(out.files[1],out.files[2], out.files[3],out.files[4], out.files[5], out.files[6]),dst_dataset="NDVI1_mosaic.hdf", separate=F, of=c("HDF4","HDF5"), verbose=TRUE)
> + gdalinfo("NDVI1_mosaic.hdf")
> + }
> Checking gdal_installation...
> GDAL command being used: "/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdalbuildvrt" -separate  "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "MOD13A3.A2013001.h20v08.005.2013042085256.hdf" "MOD13A3.A2013001.h20v09.005.2013042085401.hdf" "MOD13A3.A2013001.h21v08.005.2013042085514.hdf" "MOD13A3.A2013001.h21v09.005.2013042085607.hdf" "MOD13A3.A2013001.h22v08.005.2013042085522.hdf" "MOD13A3.A2013001.h22v09.005.2013042085359.hdf"
> Checking gdal_installation...
> GDAL command being used: "/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_translate" -of "HDF4" -of "HDF5" "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "NDVI1_mosaic.hdf"
> Input file size is 3600, 24000
> character(0)
> attr(,"status")
> [1] 1
> Warning messages:
> 1: running command '"/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_translate" -of "HDF4" -of "HDF5" "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "NDVI1_mosaic.hdf"' had status 1 
> 2: running command '"/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdalinfo" "NDVI1_mosaic.hdf"' had status 1 
> ERROR 6: GDALDriver::Create() ... no create method implemented for this format.
> 
> ERROR 4: `NDVI1_mosaic.hdf' does not exist in the file system,
> and is not recognised as a supported dataset name.
> 
> gdalinfo failed - unable to open 'NDVI1_mosaic.hdf’.
> 
> 
> Cheers
> Justin
> 
> 
> 
> 
> 
> 
> On Apr 28, 2014, at 3:17 PM, Jonathan Greenberg <jgrn@illinois.edu> wrote:
> 
>> 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@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@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@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@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@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
> 


	[[alternative HTML version deleted]]

