- 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 therasterbrick
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)
mask(crop(x, extent(m)), m)
and you do not need thesetValues
step. Sometime ago this was necessary but not any longer. You also do not need to pasteurize your polygon mask.