Cloud masking ECOSTRESS LST data

Question

I have downloaded a land surface temperature (LST) raster from ECOSTRESS as well as it's associate cloud mask (you can sign-up for free here and download your own data).

NASA's tutorial (minute 31:15) on YouTube uses Python to remove cloudy pixels from the LST image. I'd like to work with R for the same task. So, I did:

library(terra)

# working dir
setwd <- "path/"

# read the raw LST
lst_raw <- rast(paste0(wd, "ECO2LSTE.001_SDS_LST_doy2022152172444_aid0001.tif"))

# read the cloud mask
cloud_mask <- rast(paste0(wd, "ECO2CLD.001_SDS_CloudMask_doy2022152172444_aid0001.tif"))

# convert the cloud mask to binary rast
cloud_mask_binary <- ?????????

The part of converting the cloud mask to binary raster in Python is:

cloud.data = (cloud.data >> 2) & 1

I quote from the YouTube video:

To process this quality flag image into a cloud mask we are going to access the data array of this dataset with .data and we are going to bit shift this array two bits to the right and then and with one and this will produce a binary cloud mask.

What is the equivalent Python's code for converting the cloud mask to binary raster in R?

R 4.4.0, RStudio 2024.04.2 Build 764, Windows 11.

Answer 1

@Spacedman I found another solution, which seems to work okay and there is no need to install other libraries.

library(terra)

# working dir
wd <- "path/"

# read the raw LST
lst_raw <- rast(paste0(wd, "ECO2LSTE.001_SDS_LST_doy2022152172444_aid0001.tif"))
plot(lst_raw)

# read the cloud mask
cloud_mask <- rast(paste0(wd, "ECO2CLD.001_SDS_CloudMask_doy2022152172444_aid0001.tif"))
plot(cloud_mask)

 # convert the cloud mask to binary raster
    cloud_mask_binary <- (cloud_mask %/% 4) %% 2
    plot(cloud_mask_binary)

lst_masked <- ifel(cloud_mask_binary == 0, lst_raw, NA)
plot(lst_masked)

I haven't converted the LST image to Celsius nor I applied the scale factor yet, that's why the "strange" pixel values.

Answer 2

Using the bitops package you can work on values at the bit level:

> library(bitops)

First test a matrix with all 8-bit values from 0 to 255 in a 16x16 raster:

> r = rast(matrix(0:255,16,16))

We can't operate directly on this so let's copy it and work on the values of a copy:

> rmask = r
> rmask[] = rmask[] %&% 2^4
> plot(rmask)

Here's I've used the bitwise "AND" operator with binary pattern 00010000 to see if the 5th bit from the right (4th counting from zero) is a one.

Then you can compare rmask with 0 to get a 0/1 mask raster.

That python code looks like its testing against bit 3 (0-indexed) so you want to do cloudmask[] = cloudmask[] %&% 2^2 but check this (I've not downloaded the data to test). Ooh, and you can use terra::app to save the initial making a copy step:

> rmask = app(r, \(x){x %&% 4})
> plot(rmask)