How to clip MrSID Files in R and export to GeoTiff?

Question

I would like to use R to:

1. Read in an MrSID raster (.sid format)
2. Clip the extent using a shapefile
3. Save the new clipped raster as a geotiff

I am fine with steps #2 and 3, but I can't find any information about reading an MrSID file to R. Here is an example of what my workflow would look like if the raster was a .tif instead of a .sid:

# load packages
require(raster)
require(rgdal)

# Step 1: read in raster from .tif file
r <- raster("C:/.../raster.tif")
# Q: how do you do this with a .sid file?

# read in shapefile used for clip
shp <- readOGR(dsn="C:/.../folder", layer="shapefile")

# Step 2: clip raster to shapefile
r.clip <- crop(r, shp)

# Step 3: save clipped raster as tiff
writeRaster(r.clip, "C:/.../raster_clip.tif", format="GTiff")

I know this could be done relatively easily in something like QGIS, but I have a large number of shapefiles I would like to use as clipping extents and then do some follow-up data analysis, so I would like to use R to automate it.

I've found other questions regarding GDAL and MrSID (e.g. here), but none related to R, and I'm kind of a GDAL newbie so I haven't been able to successfully translate the knowledge.

Answer 1

To read MrSID-files, you need the GDAL MrSID-plugin. Check if you have the MrSID-plugin correctly installed by the command gdalinfo --format MrSID in your commandprompt/shell (otherwise check here).

Make sure the rgdal-package installs the GDAL with the MrSID-plugin (you might have several versions installed). Linking the rgdal package to the correct GDAL is not easy - I'm still figuring out how to do this on my system - but see here and here for some guidelines.

If this works, use the readGDAL function {rgdal package} to read your data.

library(rgdal)
r <- readGDAL("C:/.../raster.sid")

Alternatively, you can use the gdal_translate function in the {gdalUtils package}. The advantage is that you can choose the GDAL version easily with gdal_setInstallation.

library(gdalUtils)
gdal_setInstallation(search_path =...) # choose GDAL version with this function
r <- gdal_translate(src_dataset = "C:/.../raster.sid",
           dst_dataset = "raster.tif",
           output_Raster = TRUE) # create R raster object