Reprojecting big raster files in R
Ivan Lizarazo
27.11.2020
1. Introduction
Reprojecting big raster files is a common task needed for GIS analysis. It may be challenging to conduct it by using the projectRaster functionality provided by the raster library: it can take several hours.
In this notebook, I explore another option to overcome the challenge: the gdalUtils library which provides wrapper functions for GDAL (Geospatial Data Abstraction Library).
Let’s start by cleaning the memory and loading the required libraries:
rm(list=ls())library(raster)## Loading required package: splibrary(ggplot2)
library(gdalUtils)2. Data reading and visualization
Let’s read a Landsat 5 image:
(l5_w1 <- stack("landsat5A4-0000000000-0000000000.tif")) ## class : RasterStack
## dimensions : 9782, 10496, 102671872, 5 (nrow, ncol, ncell, nlayers)
## resolution : 4.491576e-06, 4.491576e-06 (x, y)
## extent : -74.1405, -74.09335, 4.329395, 4.373331 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : B1, B2, B3, B4, B5Let’s make a color composite:
plotRGB(l5_w1, r=4, g=5, b=3, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
Now, let’s read a second big raster file:
(l5_w2 <- stack("landsat5A4-0000000000-0000010496.tif")) ## class : RasterStack
## dimensions : 9782, 2851, 27888482, 5 (nrow, ncol, ncell, nlayers)
## resolution : 4.491576e-06, 4.491576e-06 (x, y)
## extent : -74.09335, -74.08055, 4.329395, 4.373331 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : B1, B2, B3, B4, B5cellStats(l5_w2, 'min')## B1 B2 B3 B4 B5
## 135 187 157 471 323Another plot:
plotRGB(l5_w2, r=4, g=5, b=3, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
We will read another big raster which uses a different CRS:
(photo <- stack("img_area4.tif")) ## class : RasterStack
## dimensions : 9713, 13287, 129056631, 3 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : 993012.7, 999656.2, 970493.7, 975350.2 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## names : img_area4.1, img_area4.2, img_area4.3
## min values : 0, 0, 0
## max values : 255, 255, 255What is the photo CRS?
crs(photo)## CRS arguments:
## +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1
## +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs3. Raster reprojection
Now, let’s use the gdalwarp function to reproject one file. In case of doubt about the parameters for using this function, please visit https://gdal.org/programs/gdalwarp.html.
system.time(
gdalwarp("landsat5A4-0000000000-0000010496.tif",
dstfile="l5w1rep.tif",
t_srs=crs(photo),
tr=c(0.5,0.5),
ot = 'Byte',
dstnodata = 0,
output_Raster=TRUE,
overwrite=TRUE,verbose=TRUE))## Checking gdal_installation...## Scanning for GDAL installations...## Checking Sys.which...## Checking common locations...## GDAL version 2.4.2## GDAL command being used: "/Library/Frameworks/GDAL.framework/Versions/2.4/Programs/gdalwarp" -overwrite -tr 0.5 0.5 -t_srs "+proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" -ot "Byte" -dstnodata "0" -of "GTiff" "landsat5A4-0000000000-0000010496.tif" "l5w1rep.tif"## user system elapsed
## 48.462 5.084 54.338Now, let’s reproject the second file:
system.time(
gdalwarp("landsat5A4-0000000000-0000000000.tif",
dstfile="l5w1rep.tif",
t_srs=crs(photo),
tr=c(0.5,0.5),
ot = 'Byte',
dstnodata = 0,
output_Raster=TRUE,
overwrite=TRUE,verbose=TRUE))## Checking gdal_installation...## GDAL version 2.4.2## GDAL command being used: "/Library/Frameworks/GDAL.framework/Versions/2.4/Programs/gdalwarp" -overwrite -tr 0.5 0.5 -t_srs "+proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" -ot "Byte" -dstnodata "0" -of "GTiff" "landsat5A4-0000000000-0000000000.tif" "l5w1rep.tif"## user system elapsed
## 17.785 3.110 21.8694. Check time
(l5_w2 <- stack("l5w1rep.tif")) ## class : RasterStack
## dimensions : 9718, 10466, 101708588, 5 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : 993007.8, 998240.8, 970496.2, 975355.2 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## names : B1, B2, B3, B4, B5
## min values : 0, 0, 0, 0, 0
## max values : 255, 255, 255, 255, 2555. Reproducibility
sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## 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_2.0.3.2 ggplot2_3.3.2 raster_3.3-13 sp_1.4-4
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 compiler_3.6.3 pillar_1.4.6 R.methodsS3_1.8.1
## [5] R.utils_2.10.1 iterators_1.0.13 tools_3.6.3 digest_0.6.27
## [9] evaluate_0.14 lifecycle_0.2.0 tibble_3.0.4 gtable_0.3.0
## [13] lattice_0.20-41 pkgconfig_2.0.3 rlang_0.4.8 foreach_1.5.1
## [17] rgdal_1.5-18 yaml_2.2.1 xfun_0.18 withr_2.3.0
## [21] stringr_1.4.0 dplyr_1.0.2 knitr_1.30 generics_0.1.0
## [25] vctrs_0.3.4 grid_3.6.3 tidyselect_1.1.0 glue_1.4.2
## [29] R6_2.5.0 rmarkdown_2.4.6 purrr_0.3.4 magrittr_1.5
## [33] scales_1.1.1 codetools_0.2-16 ellipsis_0.3.1 htmltools_0.5.0
## [37] colorspace_1.4-1 stringi_1.5.3 munsell_0.5.0 crayon_1.3.4
## [41] R.oo_1.24.0