1

Im trying to calculate Land surface temperature estimating with Landsat 8 images, but i got stuck in one place and i cant move on. I m following split window algorithm.
I calculated NDVI and I have to apply criteria to calculate land surface emmisivity. here my script in R studio :

#Ndvi
NDVI=(Palfa.B5-Palfa.B4)/(Palfa.B4+Palfa.B5)
NDVI <- NDVI@data@values
#vegetation fraction
NDVImin=0.2 #pre set value
NDVImax=0.5 #pre set value
Pv <- ((NDVI-NDVImin)/(NDVImax-NDVImin))^2
#roughness of the surface
F=0.55 #pre set value of geometric factor
C10 <- (1-es10) * ev10 * (1-Pv) * F
C11 <- (1-es11) * ev11 * (1-Pv) * F
#LSE for band 10
e10 <- if (NDVI<=0.2) {-0.047 * Palfa.B4+0.973} 
       else if (NDVI>0.2 && NDVI<=0.5) {ev10 * Pv+es10 * (1-Pv)+C10} 
       else (NDVI>0.5) ev10+C10

#LSE for band 11
e11 <- if (NDVI<=0.2) { -0.026 * Palfa.B5+0.984 } 
       else if (NDVI>0.2 && NDVI<=0.5) { ev11 * Pv+es11 * (1-Pv)+C11 }
       else (NDVI>0.5) ev11+C11 

something is wrong bcs i got this :

Warning message: In if (NDVI <= 0.2) { : the condition has length > 1 and only the first element will be used ev11+C11

I also tried with ifelse but then my result is vector but i need to export raster
have u got any ideas how can i solve it?
Attached you will find criteria : image enter image description here

  • Are ev10, es10, ev11 and es11 constants or rasters? – Mikkel Lydholm Rasmussen Aug 20 '18 at 14:01
  • @MikkelLydholmRasmussen ev10,es10,ev11 es111 are constants but Pv and C10 and C11 are rasters – 1milenac Aug 20 '18 at 14:15
0

I did not notice in my initial answer that NDVI was a vector created from the values of a raster. That is why you are receiving that warning message. It will not evaluate the condition for more than one element. Since you need to evaluate an entire raster and use other rasters in the calculations you can use the method shown here: https://gis.stackexchange.com/a/167514/70980

library(raster)
#Create sample rasters
NDVI <- raster(nrow=5, ncol=5)
NDVI <- setValues(NDVI, runif(ncell(NDVI)))
Pv <- setValues(NDVI, runif(25 ,0.6,0.9))
Palfa.B4 <- setValues(NDVI, runif(25 ,0.2,0.4))
C10 <- setValues(NDVI, runif(25 ,1,2))

#Stack rasters 
s <- stack(NDVI, Pv, Palfa.B4, C10)
#Sample values for constants
es10 <- 0.85
ev10 <- 0.54

#if-else reclassify statement
rc <- function(x1,x2,x3,x4) {
  #x1 - NDVI
  #x2 - Pv
  #x3 - Palfa.B4
  #x4 - C10
  ifelse(x1 <= 0.2, (-0.047 * x3 + 0.973), 
         ifelse((x1 > 0.2) & (x1 <= 0.5), (ev10 * x2 + es10 * (1 - x2) + x4), 
                ev10 + x4))
}

r.class <- overlay(s, fun=rc)

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.