[GIS] Clipping raster brick object in R

cliprraster

I have MODIS NDVI TIFF files for a year
I want to first clip it with a irregular polygon (not by its extent only) of my study area and save it as a multi band raster in TIFF format
The TIFF image are already projected to the needed coordinate system and datum (WGS 84, UTM 45N)
I used the writeRaster() function for the rasterbrick object but it gives me an error.

If anyone has faced similar problems or could you suggest me better way doing this task.

writeRaster(ndvit1, filename="ndvi_2015t1.tif", bandorder='BSQ', format="GTiff", overwrite=TRUE)
Error in rgdal::putRasterData(x@file@transient, vv, band = i, offset = off) :
Failure during raster IO
In addition: Warning message:
closing unused connection 4 (C:\Users\Ghimirebr\AppData\Local\Temp\RtmpwJschv\raster\r_tmp_2016-04-11_121550_3352_50703.gri)

## Clip a raster (brick) with the polygon (shape file) in its extent
## and mask value nan for the remaining area.

library(raster)
library(rgdal)

# Setting working directory:

setwd("C:/Working_Folder/R")
ndvi_path <- "ndvi_2015"

# Reading the shapefile 
myShapeInR <- readOGR("C:/working_folder/Boundary", 
                      "Chitwan_Nawal_Boundary")

# To view the summary of vector file
summary(myShapeInR)

# Getting the spatial extent of the shapefile 
e <- extent(myShapeInR) 

# Reading the raster you want to crop 
all_ndvi <- list.files(ndvi_path,
                       full.names = TRUE,
                       pattern = ".tif$")
ndvi_stack <- stack(all_ndvi)
ndvi_brick <- brick(ndvi_stack)

# Converting from -1 to 1 (NDVI)
ndvit1 <- ndvi_brick*0.0001

# Saving the RASTER BRICK
writeRaster(ndvit1, filename="ndvi_2015t1.tif", bandorder='BSQ',
            format="GTiff", overwrite=TRUE)
# myraster <- raster(myraster.filename) 

# Cropping the raster to the shapefile spatial extent 
ndvit1.crop <- crop(raster, e, snap="out") 

# Dummy raster with a spatial extension equal to the cropped raster, 
# but full of NA values 
crop <- setValues(ndvit1.crop, NA) 

# Rasterize the catchment boundaries, with NA outside the catchment boundaries 
myShapeInR.r <- rasterize(myShapeInR, crop) 

# Putting NA values in all the raster cells outside the shapefile boundaries 
ndvit1.masked <- mask(x=ndvit1.crop, mask=myShapeInR.r)

Best Answer

In order to write a raster stack or brick, the Interleave band can be used in the option argument.

    writeRaster(ndvi_stack, filename="<filename>.tif", options="INTERLEAVE=BAND", overwrite=TRUE)

Related Solutions

[GIS] Really slow extraction from raster even after using crop

1) resample results in 50% improvement

I was able to get about 50% improvement by resampling directly from the cld raster to a new raster with the same extent/resolution as r and a nearest neighbor sampling method:

system.time({
  mat<-as.data.frame(getValues(r))
  mat$landuse<- NA
  mat$landuse<-getValues(resample(cld,r,method='ngb'))
})
   user  system elapsed 
   2.60    0.00    2.61

vs.

system.time({
  mat<-as.data.frame(getValues(r)) # getting values from the stack
  xy<-xyFromCell(r,c(1:ncell(r)),spatial = TRUE)
  cells<-cellFromXY(clip1,xy)
  mat$landuse<- NA
  mat$landuse<- extract(clip1,cells) #this line takes 5 mins based on profiling
})
   user  system elapsed 
   4.98    0.00    5.02

On a smaller dataset, and with a much smaller memory footprint

2) Parallelization could get you a lot more

That will improve things significantly but you can get massive improvement if you can parallelize this. R comes with a couple parallelization backends and they all run through foreach. I assume you are going to either process mat in place or save it for later. Since it takes so much work to get that resampled data let's just assume we'll save it for later. The most convenient form is probably a raster alongside the data_robin files.

Unfortunately, Windows and Unix parallelization options differ. On linux, use doMC, on Windows use doSNOW. Assuming we employ 4 workers:

linux initialization:

library(doMC)
registerDoMC(4) # number of workers should be less than number of CPU cores

windows initialization:

library(doSNOW)
cluster<-makeCluster(4, type = "SOCK") # num workers should be < num CPU cores
registerDoSNOW(cluster)

library(foreach)
library(tools)

# Assume you have an array of filenames called 'files'
foreach (i=1:length(filenames), .packages=c('raster')) %dopar% {
    r <- stack(paste0(path, "/data_robin/", files[i]))
    outFilename=paste0(path, "/data_robin/", file_path_sans_ext(files[i]), "_cld.tif")
    cldResampled <- resample(cld,r,method='ngb')
    writeRaster(cldResampled, filename=outFilename, format="GTiff")
}

One of the drawbacks of the parallel foreach is that it's hard to tell when something goes wrong. It would be good to do this serially first by replacing the %dopar% with %do% until you know it is working, and then let it run through the whole thing.

Caveats In my simple example above (each raster had 1/100 the pixels in cld and r, respectively) I only improved an additional 30% by engaging 5 workers over doing it serially with just the single process. I was unable to parallelize the example with large rasters without getting Error in { : task 1 failed - "cannot allocate vector of size xx.x Mb". I think conceptually this should work but I wasn't able to get it working at the scale you are working at.

[GIS] Crop Raster from using Polygon shapefile in R

It looks like you have a single point in the shapefile you are trying to crop to and the error message isn't exactly right. Lets' try and simulate that with a reproducible example:

b <- brick(system.file("external/rlogo.grd", package="raster"))
p=data.frame(x=50,y=35) # a point in the middle of b
coordinates(p)=~x+y

Then:

> extent(b)
class       : Extent 
xmin        : 0 
xmax        : 101 
ymin        : 0 
ymax        : 77 
> extent(p)
class       : Extent 
xmin        : 50 
xmax        : 50 
ymin        : 35 
ymax        : 35 
> crop(b, extent(p), snap="out")
Error in .local(x, y, ...) : extents do not overlap

Perhaps the crop algorithm is computing the area of overlap, getting zero, and assuming the areas don't overlap. Anyway, its a minor bug in crop, I'd say. To check, lets just make sure it works with a larger extent:

> p2=data.frame(x=c(50,51),y=c(35,38))
> coordinates(p2)=~x+y
> crop(b, extent(p2), snap="out")
class       : RasterBrick 
dimensions  : 3, 1, 3, 3  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 50, 51, 35, 38  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : red, green, blue 
min values  : 155,   165,  222 
max values  : 157,   168,  224

Of course we can't be sure this is your situation without your data, but it looks like it. Track back your SE.geo object and see what the original is like:

se <- readOGR(dsn = "C:/ArquivosR/EVI/Dados_anuais/2000", layer = "Estado")

Best Answer

Related Solutions

[GIS] Really slow extraction from raster even after using crop

[GIS] Crop Raster from using Polygon shapefile in R

Related Question