extractBits | R Documentation |
This function applies bitops::bitAnd()
and bitops::bitShiftR()
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()
.
extractBits( x, bitShift = 2, bitMask = 15, filename = "", datatype = "INT1U", NAflag = 255, ... ) makeWeights( x, bitShift = 2, bitMask = 15, threshold = NULL, decodeOnly = FALSE, filename = "", datatype = "INT1U", NAflag = 255, ... ) maskWater( X, bitShift = NULL, bitMask = NULL, keep = NULL, datatype = "INT1U", NAflag = 255, ... )
x |
|
bitShift |
|
bitMask |
|
filename |
|
datatype |
|
NAflag |
|
... |
Other arguments passed to |
threshold |
|
decodeOnly |
|
X |
|
keep |
If |
A Raster*
object.
extractBits()
: Extract bit-encoded information from Raster*
file
maskWater()
: Masks water (additional information required)
makeWeights()
and extractBits()
are identical with the only difference
that makeWeights()
does additionally convert the data into weighting
information.
Matteo Mattiuzzi
detectBitInfo()
.
## Not run: # example MOD13Q1 see MODIS Vegetation Index User's Guide (MOD13 Series; # https://lpdaac.usgs.gov/sites/default/files/public/product_documentation/mod13_user_guide.pdf) # enter in Layers # See in TABLE 5: Descriptions of the VI Quality Assessment Science Data Sets (QA SDS). # 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) } ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.