loicdtx / bfastspatial Goto Github PK

I would like to have it as the "short description", which tag is that?

RasterStack with time written to z dimension

A bit of control flow can be added, as well as visualization in bfmSpOver example (polygons case) can be improved

Line 42: prescribed error message inside getSceneinfo() refers to a global variable (sourcefile) inside an mapply() function.

...since we have been copied these code snippets throughout various functions, why not have them as separate auxiliary functions which are then called by others?

In my version of bfmSpatial() I embedded an error flag of -9999 in both the breakpoint and magnitude rasters using tryCatch() to catch possible errors. This is useful in that it allows error catching while still being able to follow up on the error (ie. target the pixel(s) where the error is encountered for follow-up analysis or de-bugging). If only c(NA, NA) is assigned, there is no [convenient] way of discriminating pixels originally masked out and those where errors were encountered.

However, -9999 could be problematic if bfmSpatial() is run on rescaled integer values (as I currently do), since magnitude values in this case could be within that range for certain indices/bands. In that case, you could consider writing another RasterLayer in the output brick simply indicating where errors were encountered (with a value of 1 for errors and 0 or NA for successful runs, for example).

Right now a time vector has be supplied externally. 2 additonal options should be added

• check if there's anything in the Z dimention
• parses through layer names to retrieve date

Perhaps one idea for reproducibility is to crop a small extent out of a bunch of reflectance .hdf and write them to geoTiff (or tar.gz) in a data directory of the package. (That would allow to write reproducible examples for hdf2ndvi (once tiff support is added), processLandsat, and processLandsatBatch)

Sr2vi needs to retrieve sensor info using getsceneinfo and pass it as argument to .tc* function

I ran into your LandsatVIs.R calculation of the Tasseled Cap for Landsat TM while trying to figure out the correct coefficients myself. The PDF scan of the Crist 1985 paper from RSE is so horrible that it's hard to see if some coefficients are positive or negative. I found some other papers documenting the coefficients and even found an alternative source for the Crist 1985 that is digitized in better resolution, and these sources show that some of the coefficients your code uses have the wrong sign.

The Crist 1985 rescan is available here and Meddens et al (2013) also published the coefficients in their work available here.

The Wetness coefficient for Band 5 should be -0.6806 and the Greenness coefficient for Band 5 should be -0.0002. Otherwise everything looks correct! I've used the wrong coefficients many times because of the poor quality of the PDF so hope I can help you solve the same small problem.

Cheers

Currently only one function to calculate ndvi + crop and mask. That function could be extend in order to support a whole set of indices, supplied as argument
SAVI, DI, NBR, EVI, TC, NDVI

First step is to write a method and a template for that to work, later anybody can contribute his own indices

Sometimes the function does not get sourced.
Is it only when installing using install_github(), or does it also happen when building locally?

Probably sr2vi.R or LandsatVIs.R need to be adapted accordingly

Main 'bottleneck' is the hdf2ndvi function and perhaps also processLandsat

If x is a character (Landsat filename), the following line extracts the Landsat code in most cases. (lndsr..hdf, ndvi..grd, *.tar.gz, etc). That would increase the flexibility of the getSceneinfo function.

str_extract(string=basename(x), '(LT4|LT5|LE7)\d{13}')

str_extract(string='/media/jdifu/lndsr.LE72300712000262COA00.hdf', '(LT4|LT5|LE7)\d{13}')
[1] "LE72300712000262"

str_extract(basename('/media/blahblah/LE72300712000262-SC20131118103555.tar.gz'), '(LT4|LT5|LE7)\d{13}')
[1] "LE72300712000262"

From espa

NOTICE: August 19, 2014 -- Upcoming changes to ESPA output file names and metadata

By the end of August, the "fmask" band's file name suffix will permanently change to "cfmask" to accurately reflect the algorithm used in ESPA. 

Example: Cloud mask for scene LE70390372008210EDC00 will be named "LE70390372008210EDC00_cfmask" instead of "LE70390372008210EDC00_fmask".

Please update any programs or scripts you may be using to reflect this change.

Changes to the cloud mask metadata will also be made:
   1. The "name" will be "cfmask" instead of "fmask"
   2. The "short_name" will be "LXXCFMASK" instead of "LXXFMSK", where LXX is LT4, LT5 or LE7
   3. The "long_name" will be "cfmask_band" instead of "fmask_band""

Some users may not want to have the layers in the ts stack ordered chronologically, but rather alphabetically by filename so that the order of the layers reflect the order in which the individual layer files are listed in the directory (using list.files). I suggest making this step optional.

Suggested pages for the wiki:

1. Landsat Pre-processing
2. Running BFAST Monitor on a raster time series (including bfmPixel, bfmSpatial)
3. Post-processing and visualization (include code migrated from bendv/bfastSpatial)

These can eventually be copied to a tutorial (vignette) to be distributed with the completed package.

ndvi, evi, savi, nbr, etc present in (at least) some datasets delivered by USGS. Possible to make the functions a bit flexible, choosing the already pre-processed layer when it already exists and claculating it when it does not exists.

Are you done with it?
To delete remote branch:
git push origin --delete

add table of content to the tutorial
update references

To be passed to processLandsat. (If necessary)

From processLandsat example

# Get list of test data files
dir <- system.file('external', package='bfastSpatial')
list <- list.files(dir, full.names=TRUE)

# Set the location of output and intermediary directories (everything in tmpdir in that case)
srdir <- dirout <- file.path(rasterOptions()$tmpdir, 'bfmspatial')
dir.create(dirout, showWarning=FALSE)
# Generate (or extract, depending on whether the layer is already in the archive or not) NDVI for the first archive file
processLandsat(x=list[1], vi='sr_band2', outdir=dirout, srdir=srdir, delete=TRUE, mask='fmask', keep=0, overwrite=TRUE)

...just opening this issue to brainstorm a bit about the most effective way to combine raster time series from neighbouring scenes.

For example, if your AOI spans 4 different scenes, there should be a convenient way to store the pre-processed time series for that AOI in one object (RasterBrick or otherwise).

A single brick is possible, but is complicated by the fact that acquisition dates will differ (maybe only slightly) between scenes. I can think of a number of options:

1. a custom grid-based object where each cell has a pixel value and a date associated with it. This seems convenient, but would take up more storage space (to accommodate the date info per grid cell).
2. a extended RasterBrick, where NA's are assigned to cells from neighbouring scenes with no data acquired at that particular date. This seems a bit clumsy to me, but would work...
3. keep them as separate bricks and process them independently of each other. But then, how would we handle overlapping cells?

In cases where the scenes are combined (option 1 or 2), overlapping cells can easily be incorporated by simply 'densifying' the time series at those locations. In some cases, the pixel values will have been acquired on the same date, in which case either a summary value (mean) is used, or preference is given to one of the overlapping values.

Ideas?

Migrate code from bendv/bfastSpatial (and close that repo). Mostly post-processing, with a few bfm helper functions and some that could also contribute to pre-processing workflow.

See @describeIn and @rdName roclets
http://cran.r-project.org/web/packages/roxygen2/vignettes/rd.html

Also make clear what are the platform limitations (hdf driver, etc)

man/cleanMODIS.Rd:17: unknown macro '\textbf'

...rather than having a 'permanently masked' input time series brick. This could be useful in testing different forest masks (made in different years or using different methods), without having to re-create the input raster brick each time.

If a forest mask is supplied (as a raster layer with same spatial attributes as input time series brick -- ie. compareRaster(input, forest_mask) == TRUE), then an additional control flow sequence could be added which basically decides whether bfm should be applied to a pixel based on the value of the accompanying forest mask.

Renaming the function would also make sense...

see http://stackoverflow.com/questions/4124900/is-there-a-way-to-use-two-statements-in-a-function-in-r

This could be useful in (e.g.) clumpSize()

Error in xj[i, , drop = FALSE] : subscript out of bounds

In addition: Warning message:
In .doExtract(x, i, drop = drop) : some indices are invalid (NA returned)

0's are currently treated as NA's, so islands of pixels with value 0 are left alone. This should be changed.

Could this be due to the call to mosaic() in mc.calc

Is it an idea to have one vignette containing several chapters (data download, pre-processing, data inventory, data exploration, post-processing), instead of multiple thematic vignettes. Latex has a command (\include{file}) that allows to have individual files for different sections. That may make things easier for git and several people write on the same overall document at the same time (since they would most likely be working on different sections and therefore on different files)

listOut <- c(listOut, dots) prior to the do.call() call should work. To be tested.

I had a look at the USGS test data as delivered in GeoTiff; they do not order at all the same way as in the .hdf . The layer subsetting has to be done differently (for SR bands and also for cloud mask layer selection). Probably using some pattern matching technique. So the return of the LandsatVIs would become c(band3, band4, etc) instead of c(3, 4). This would also make the hdf case less prone to errors whenever USGS decides to change sds orders, etc

See e-mails from TOm

processLandsat() uses a sequence of conditions to check which vi to process. Wouldn't match.fun() be a cleaner option to do that.

Current method

                if(vi == 'ndvi') {
                    viFormula <- .ndvi()
                } else if(vi == 'evi') {
                    viFormula <- .evi()
                } else if(vi == 'nbr') {
                    viFormula <- .nbr()
                } else if(vi == 'savi') {
                    viFormula <- .savi(L=L)
                } else {
                    stop("Unsupported vi")
                }

Eg formula = 'band1' + 'band2'

USe

y <- as(y, 'SpatialPolygons')

Instead of

y <- polygonFromExtent(y)

since clump() only takes individual layers as arguments, this could be wrapped in a for or foreach for raster bricks

...instead of dropLayer() when subseting parameters are supplied. Also, check the effect this has on speed (especially for bfmSpatial()).