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

Jonathan Greenberg jgrn at illinois.edu
Thu May 1 16:30:54 CEST 2014


That's a good trick -- I'll test to see how it performs vs. the
two-step method I was using and maybe tweak gdalUtils.  Thanks!

--j

On Wed, Apr 30, 2014 at 7:49 AM, Loïc Dutrieux <loic.dutrieux at wur.nl> wrote:
> Dear Justin,
>
> gdalwarp is able to reproject/resample and mosaic all at once, and that is
> probably the preferred option in that case (since there are no overlapping
> areas).
> See the simplified example below:
>
>
> ###
> library(raster)
> library(gdalUtils)
>
> dir <- file.path(rasterOptions()$tmpdir, 'MODIS')
> dir.create(dir, recursive=TRUE)
>
> # Download two neighboring modis tiles
> download.file('http://e4ftl01.cr.usgs.gov/MOLT/MOD13A2.005/2000.04.22/MOD13A2.A2000113.h20v08.005.2006275210548.hdf',
> destfile=file.path(dir, 'NDVI.h20v08.hdf'))
> download.file('http://e4ftl01.cr.usgs.gov/MOLT/MOD13A2.005/2000.04.22/MOD13A2.A2000113.h21v08.005.2006275210549.hdf',
> destfile=file.path(dir, 'NDVI.h20v09.hdf'))
>
> # Get list of files downloaded
> list <- list.files(dir, full.names=TRUE)
> # Get list of sds1
> sdsList <- sapply(X=list, FUN=function(x){get_subdatasets(x)[1]})
> # Warp
> gdalwarp(srcfile=sdsList, t_srs="+proj=longlat +datum=WGS84 +no_defs",
> dstfile=file.path(dir, 'mosaic_latlong.tif'))
>
> r <- raster(file.path(dir, 'mosaic_latlong.tif'))
> plot(r)
>
> Hope this helps
> Best regards,
>
> --
> Loïc Dutrieux
> Laboratory of Geo-Information Science and Remote Sensing
> Wageningen University
> The Netherlands
>
>
>
>
> On 14-04-30 01:10 AM, Justin Michell wrote:
>>
>> 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 at coords)
>>
>>
>>   extract(NDVI1, sitesPolyPoints at 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 at 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 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
>>>
>>>
>>
>>
>>          [[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
>
>
> _______________________________________________
> 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