extractBits: Extract Bit-Encoded Information and Create Weights Raster

Description

This function applies bitAnd and bitShiftR from bitops to convert bit-encoded information. It is also possible to convert this information to a scale from 0 to 1 in order to use it as weighting information in functions like whittaker.raster or smooth.spline.raster.

Usage

extractBits(x, bitShift = 2, bitMask = 15, filename = "", ...)
makeWeights(x, bitShift = 2, bitMask = 15, threshold = NULL,
  filename = "", decodeOnly = FALSE, ...)
maskWater(X, bitShift = NULL, bitMask = NULL, keep = NULL,
  datatype = "INT1U", ...)

Arguments

matrix, vector or Raster* object.

bitShift

integer. Bit starting point, see examples and detectBitInfo.

bitMask

integer. Bit mask size, see examples and detectBitInfo.

filename

character passed to writeRaster. If not specified, output is written to a temporary file.

...

Other arguments passed to writeRaster.

threshold

integer. Threshold for valid quality.

decodeOnly

logical. If FALSE (default), convert bits to weights from 0 (not used) to 1 (best quality). If TRUE, only extract selected bits and convert to decimal system.

Raster* object.

keep

If NULL (default), bits are only encoded, else an integer vector of values you want to keep (becomes TRUE), the rest becomes NA. See examples.

datatype

character. Output datatype, see writeRaster.

Value

A Raster* object.

Functions

extractBits: Extract bit-encoded information from Raster* file
maskWater: Masks water (additional information required)

Examples

Run this code

# NOT RUN {
# example MOD13Q1 see https://lpdaac.usgs.gov/products/modis_products_table/mod13q1 
# enter in Layers
# See in TABLE 2: MOD13Q1 VI Quality
# column 1 (bit) row 2 VI usefulness 
bitShift = 2
bitMask = 15 # ('15' is the decimal of the binary '1111')  
# or try to use
detectBitInfo("MOD13Q1") # not all products are available!
viu <- detectBitInfo("MOD13Q1","VI usefulness") # not all products are available!
viu

# simulate bit info
bit <- round(runif(10*10,1,65536))

# extract from vector
makeWeights(bit,bitShift,bitMask,decodeOnly=TRUE)
# the same as
extractBits(bit,bitShift,bitMask)

# create a Raster object
VIqual <- raster(ncol=10,nrow=10)
VIqual[] <- bit

# extract from Raster object
a <- makeWeights(VIqual,bitShift,bitMask,decodeOnly=TRUE)

# linear conversion of 0 (0000) to 15 (1111) to 1 fo 0 
b <- makeWeights(VIqual,bitShift,bitMask,decodeOnly=FALSE)

threshold=6 # every thing < threshold becomes a weight = 0
c <- makeWeights(VIqual,bitShift,bitMask,threshold,decodeOnly=FALSE)

res <- round(cbind(a[],b[],c[]),2)
colnames(res) <- c("ORIG","Weight","WeightThreshold") 
res

#####
# water mask
tif = runGdal(product="MOD13A2",begin="2009001",end="2009001", extent=extent(c(-9,-3 ,54,58)),
              SDSstring="001",job="delme") # 6.4 MB
x <- raster(unlist(tif))

res1 <- maskWater(x)
plot(res1)

res2 <- maskWater(x,keep=1) # 1 = Land (nothing else)
x11()
plot(res2)

# Land (Nothing else but land) + Ocean coastlines and lake shorelines + shallow inland Water,
# the rest becomes NA
x11()
res3 <- maskWater(x,keep=c(1,2,3)) 
plot(res3)

###############

# as on Linux you can read HDF4 directly you can also do:
if(.Platform$OS.type=="unix")
{
  x <- getHdf(HdfName="MOD13A2.A2009001.h17v03.005.2009020044109.hdf", wait=0) # 6.4 MB
  
  detectBitInfo(x) # just info
  getSds(x) # just info
  
  x <- getSds(x)$SDS4gdal[3] # you need 'VI Quality'
  x  <- raster(x) 
  # plot(x)
  # ex <- drawExtent()
  ex <- extent(c(-580779,-200911,5974929,6529959))
  x  <- crop(x,ex) # just for speed-up
  
  res1 <- maskWater(x)
  plot(res1)
  
  res2 <- maskWater(x,keep=1) # 1 = Land (Nothing else but land), the rest becomes NA
  x11()
  plot(res2)
  
  # Land (Nothing else but land) + Ocean coastlines and lake shorelines + shallow inland Water,
  # the rest becomes NA
  res3 <- maskWater(x,keep=c(1,2,3))
  x11()
  plot(res3)
}
# }
# NOT RUN {
# }