[GIS] hdf to gtiff, reprojecting, cropping with a shapefile in batch using R

coordinate systemgdalgeotiff-tiffhdfr

I have some raster data from MODIS, with hdf format.
I want to convert them into GeoTIFF, reproject them and crop them based on my shapefile.
I was able to do it in GDAL, in batch mode, with the following code, but I don't know how to do it in R. 

gdal_translate -of GTiff input.hdf output.tif
gdalwarp -t_srs "+proj=utm +zone=16 +datum=WGS84" input.tif output.tif
gdalwarp -cutline tl_2013_18_cousub.shp -crop_to_cutline -dstalpha input.tif OUTPUT.tif

I use the following code in R, but it's just for the converting to GeoTIFF and also it does not provide the result with TIFF format. (I don't know why)

    library(gdalUtils)

files <- dir(pattern = ".hdf")
for (i in 1:length(files)){
  sds <- get_subdatasets(files[i]) # the input files have different subdatatsets, and I need only the first one
  gdal_translate(sds[1], dst_dataset = paste0("new_",basename(files[i])), of = "GTiff")
}

based on your comment, I wrote the following code though it has errors,
can anyone help me with the issues in comments?

   library(gdalUtils)
        files <- list.files(path=".../test", pattern="*.hdf", full.names=T, recursive=FALSE)
    for (i in 1:length(files)){
    sds <- get_subdatasets(files[i])
    gdal_translate(sds[1], dst_dataset = paste0(basename(files[i]),".tif"), of = "GTiff")
    } 
#first issue: the file name would be like "filename.hdf.tif" how can I remove "hdf" part to make it like "filename.tif"

    shp <- shapefile("tl_2013_18_cousub.shp")

#second issue: I had to write another loop, can it be merged with previous loop to save time in running process?
    files <- list.files(path=".....", pattern="*.tif", full.names=T, recursive=FALSE)
    for (i in 1:length(files)){
    sds <- raster(files[i])
    xp <- projectRaster(sds, crs = "+proj=utm +zone=16 +datum=WGS84") 
    writeRaster(xp, file = paste0("prj_",basename(files[i])))
# third issue: the saved file has different resolution! how can I maintain the original resolution 

    cr <- crop(xp, shp)  
    writeRaster(cr, file = paste0("crp_",basename(files[i])))
# Forth issue: it guves me error: "Error in .local(x, y, ...) : extents do not overlap"
    }

Best Answer

That command line code in R roughly is 

library(raster)
x <- raster("input.hdf")  ## 1
xp <- projectRaster(x, crs = "+proj=utm +zone=16 +datum=WGS84") # 2
shp <- shapefile("tl_2013_18_cousub.shp")
crop(xp, shp)  # 3

But (and there are endless more caveats to these caveats):

  1. requires rgdal and your system's rgdal might not have HDF4 support, and it might have multiple subdatasets so require sds string rather than file
  2. Better to have a target raster with defined extents, crs, and resolution, otherwise raster (and GDAL also) will auto-guess extents/resolution for you
  3. Really we should spTransform(sp, projection(xp)) rather than the lossy raster reprojection

I'd be trying

library(raster)
x <- raster("input.hdf")  ## 1
shp <- shapefile("tl_2013_18_cousub.shp")
crop(xp, spTransform(shp, projection(xp))
Related Question