0

I am using mosaic_rasters from gdalUtils to combine different raster files. More specifically, I am using the 12 tiles that cover Southern Africa from the well-known Hansen et al. (2013) forest map (1.1) that can be downloaded here: http://earthenginepartners.appspot.com/science-2013-global-forest/download_v1.1.html. I am using the tree cover, forest gain and forest loss data. All is working fine for the forest gain and forest loss files but I receive an error when combing the tree cover files (which are the largest).

Below you will find my code. Note that "tree_cover_files" is a vector that points towards the 12 tiles stored in a local folder. I am receiving the following error after some time:

Checking gdal_installation...
Scanning for GDAL installations...
Checking the gdalUtils_gdalPath option...
GDAL version 2.2.4
GDAL command being used: "C:\OSGeo4W64\bin\gdalbuildvrt.exe" -input_file_list "c:\Temp\RtmpkRd30l\file4687024480a.txt" "c:\Temp\RtmpkRd30l\file4684596162f.vrt" 
Checking gdal_installation...
Scanning for GDAL installations...
Checking the gdalUtils_gdalPath option...
GDAL version 2.2.4
GDAL command being used: "C:\OSGeo4W64\bin\gdal_translate.exe" -of "GTiff" "c:\Temp\RtmpkRd30l\file4684596162f.vrt" "P:/globiom/Projects/ISWEL/data/forest/combined_tiles/tree_cover.tif"
Input file size is 120000, 1600000...10...20...30...40...
ERROR 1: TIFFFillStrip:Read error at scanline 39921; got 3204 bytes, expected 10047
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: P:/globiom/Projects/ISWEL/data/forest/tree_cover/Hansen_GFC2014_treecover2000_10S_030E.tif, band 1: IReadBlock failed at X offset 0, Y offset 39922
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 39922

I understand that this might be thread error, related to how the files are read and intermediate vrt file is constructed (https://github.com/mapnik/node-mapnik/issues/437 ) and could be solved by setting "VRT_SHARED_SOURCE" to 0. I tried to do this using setCPLConfigOption("VRT_SHARED_SOURCE", "0") in R but still receive the same error. Is this really the way to solve this or is the file perhaps corrupt (unlikely as this dataset is used by many people - I also downloaded them twice)? I hope somebody can give me advice on how to make this work.

library(purrr)
library(curl)

### FUNCTIONS
# Function to download data
# NB: need to add mode = "wb" for binary files otherwise the tif will be   corrupted!!!!
download_f <- function(url, folder){
  print(url)
  dstfile <- file.path(proj_path, paste0(folder, basename(url)))
  if(!file.exists(dstfile))  download(url, mode = "wb", destfile = dstfile)
}

### DOWNLOAD DATA
#source:  Hansen et al. (2013), http://science.sciencemag.org/content/342/6160/850
#Website: http://earthenginepartners.appspot.com/science-2013-global-forest/download_v1.1.html

# The files are at a 30m resolution and therefore very large. We only download the tiles for Southern Africa using the following
# NS and EWE combinations
NS <- c("00N", "10S", "20S", "30S")
EW <- c("010E", "020E", "030E")

NS_EW <- expand.grid(NS = NS, EW = EW, stringsAsFactors = F)
NS_EW <- paste(NS_EW$NS, NS_EW$EW, sep = "_")
NS_EW <- paste(NS_EW, collapse = "|")

## Treecover
# list of urls for treecover tiles
url <- "http://commondatastorage.googleapis.com/earthenginepartners-hansen/GFC2014/treecover2000.txt"

# Select relevant tiles
filenames = getURL(url)
filenames = strsplit(filenames, "\n")[[1]]
tiles = grep(NS_EW, filenames, value = T)

# Set location
# Create "data" subfolder manually 
proj_path <- "C:/temp"

# download
walk(tiles, download_f, "data/")


### MOSAIC    
e <- extent(20, 40, -20, -10)
template <- raster(e)
proj4string(template) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") 
writeRaster(template, file = file.path(proj_path, "data/tree_cover.tif"), format="GTiff", overwrite = T)

# Create list of tif files
tree_cover_files <- list.files(file.path(proj_path, "data"), pattern = glob2rx("*.tif"), full.names = T)

setCPLConfigOption("VRT_SHARED_SOURCE", "0")
mosaic_rasters(gdalfile = tree_cover_files,
           dst_dataset = file.path(proj_path, "data/tree_cover.tif"), of="GTiff",
           verbose = T)
  • Have you got 12 files of the form Hansen_GFC2014_treecover2000_10S_070W.tif ? And each one is about 150Mb big? And you are trying to make one TIFF over the whole area? Is the result going to be too big for GDAL? Can you make a fully reproducible example including the URLs of the files, because its going to take ages to download those for me and I don't want to waste my time. – Spacedman Mar 26 at 8:34
  • I have added a piece of code that I used to download the relevant tiles (12) to my HD fur further processing. If you create the relevant temporary folder and subfolder, it should automatically download the relevant files after which you can run the code I wrote to mosaic the files. Perhaps the resulting file is too large for GDAL. That would at least be an explanation why it does not work. – Michiel van Dijk Apr 1 at 7:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.