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

Jonathan Greenberg jgrn at illinois.edu
Tue Apr 29 17:50:57 CEST 2014


(cc'ing to r-sig-geo to preserve the info):

Justin:

Right now, mosaic doesn't work directly with HDFs (although I could do
a little tweaking to fix that, I suppose -- I'll look into that) --
you'll want to go ahead and first batch convert the HDFs to e.g.
GeoTIFFs via gdal_translate or (even easier) batch_gdal_translate,
then run the mosaic on the geotiffs.  Let me know how it goes!

--j

On Tue, Apr 29, 2014 at 7:56 AM, Justin Michell <jwm302 at 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 at 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 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
>



-- 
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