## ----boundary, messages=FALSE-------------------------------------------------
library(terra)
library(geodata)
ken <- gadm(country="KEN", level=1, path=".")
pol <- ken[ken$NAME_1 == "Marsabit", ]


## ----prj----------------------------------------------------------------------
datadir <- file.path(dirname(tempdir()), "_modis")
mf <- file.path(datadir, "modis_qualmasked.tif")
rmask <- rast(mf)

prj <- crs(rmask)

prj

poly <- project(pol, prj)


## ----crop---------------------------------------------------------------------
rcrop <- crop(rmask, poly)


## ----cropmap------------------------------------------------------------------
plotRGB(rcrop, r = 2, g = 1, b = 4,  main = "False color composite", stretch = "lin" )
lines(poly, col="blue")


## ----clamp--------------------------------------------------------------------
r <- clamp(rcrop, 0, 1) 


## ----ndvi---------------------------------------------------------------------
ndvi <- (r[[2]] - r[[1]]) /(r[[2]] + r[[1]])
plot(ndvi, main="NDVI")