Skip to content

6 Recipes

Mateusz Żółtak edited this page Oct 16, 2019 · 1 revision

Here you can find configuration examples for running a whole data processing workflows for few sample output indicators.

Please read the init step documentation first.

Monthly NDVI composites

To compute monthly NDVI composites we must run the composite step on the NDVI band using the n type indicators computed by the which step.

The which step also depends on the NDVI band.

To compute the NDVI band we must run the indicator step which depends on the mask step.

And of course as for all workloads the init, download and tile steps must be run.

Therefore we need to configure and run init, download, mask, indicator, which, composite and tile steps.

Configuration

To compute NDVI we need bands B04 and B08 and to compute the mask an SCL band is required. This gives us bands parameter value.

bands = c('B04', 'B08', 'SCL')

For the mask step we assume a simple configuration leaving only vegetation, non-vegetation and water pixels.

maskParam = list(
  list(bandName = 'MASK', minArea = 0L, bufferSize = 0L, invalidValues = c(0L:3L, 7L:11L), bufferedValues = integer())
)

The NDVI is computed by (B04 - B08) / (B04 + B08). We want it in the native Sentinel 2 resolution of 10m. Taking into account corner cases (division by 0) the indicator step config is:

indicatorIndicators = list(
  list(bandName = 'NDVI',  resolution = 10, mask = 'MASK', factor = 10000, bands = c('A' = 'B04', 'B' = 'B08'), equation = '(A.astype(float) - B) / (0.0000001 + A + B)')
)

Other steps are straightforward to configure - we just need to pick up band names and put them into right places.

A common config file for all scripts involved should look as follows (let's assume it's saved in myConfig.R):

cubeRpath = 'ADJUST'
gridFile = 'ADJUST'
tmpDir = 'ADJUST'
rawDir = 'ADJUST'
periodsDir = 'ADJUST'
tilesDir = 'ADJUST'
cacheTmpl = 'ADJUST/{region}_{dateFrom}_{dateTo}_{cloudCovMax}_{bands}'

bands = c('B04', 'B08', 'SCL')
cloudCov = ADJUST
nCores = ADJUST
chunksPerCore = 10

dwnldMethod = 'download'
dwnldMaxRemovals = 1000
dwnldNCores = 4
dwnldTimeout = 120
dwnldSkipExisting = 'samesize'
dwnldTries = 3

maskSkipExisting = TRUE
maskParam = list(
  list(bandName = 'MASK', minArea = 0L, bufferSize = 0L, invalidValues = c(0L:3L, 7L:11L), bufferedValues = integer())
)

indicatorSkipExisting = TRUE
indicatorIndicators = list(
  list(bandName = 'NDVI',  resolution = 10, mask = 'MASK', factor = 10000, bands = c('A' = 'B04', 'B' = 'B08'), equation = '(A.astype(float) - B) / (0.0000001 + A + B)')
)

whichSkipExisting = TRUE
whichBands = c('NDVI')
whichPrefix = 'NMAX'
whichDoyPrefix = 'DOYMAX'
whichBlockSize = 512

compositeSkipExisting = TRUE
compositeBands = list(
  band      = c('NDVI'),
  whichBand = c('NMAXNDVI'),
  outBand   = c('NDVI')
)
compositeBlockSize = 512

tileSkipExisting = TRUE
tileRawBands = character()
tilePeriodBands = list(
  '1 month' = c('NDVI'),
)
tileResamplingMethod = 'near'
tileGdalOpts = '--config GDAL_CACHEMAX 1024 -co "COMPRESS=DEFLATE" -co "TILED=YES" -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512"'

Workflow

Let's assume the processing is done for the region myRegion for the whole year 2018.

Rscript      init.R myConfig.R myRegion 2018-01-01 2018-13-31 {yourApiLogin} {yourApiPswd}
Rscript     dwnld.R myConfig.R myRegion 2018-01-01 2018-13-31
Rscript      mask.R myConfig.R myRegion 2018-01-01 2018-13-31
Rscript indicator.R myConfig.R myRegion 2018-01-01 2018-13-31
Rscript     which.R myConfig.R myRegion 2018-01-01 2018-13-31 '1 month'
Rscript composite.R myConfig.R myRegion 2018-01-01 2018-13-31 '1 month'
Rscript      tile.R myConfig.R myRegion 2018-01-01 2018-13-31

Remarks:

  • The dwnld.R, mask.R and indicator.R scripts can be run regularly, e.g. daily or weekly to spread the workload over a longer period of time.
  • The which.R, composite.R and tile.R scripts can be run monthly.

Yearly NDVI centiles and day of year with a maximum NDVI

To compute NDVI centiles we must run the aggregate step on the NDVI band.

For the day of year with maximum NDVI we must run the which step which also depends on the NDVI band.

To compute the NDVI indicator we must run the indicator step which depends on the mask step.

And of course as for all workloads the init, download and tile steps must be run.

Therefore we need to configure and run init, download, mask, indicator, which, aggregate and tile steps.

Configuration

To compute NDVI we need bands B04 and B08 and to compute the mask an SCL band is required. This gives us bands parameter value.

bands = c('B04', 'B08', 'SCL')

For the mask step we assume a simple configuration leaving only vegetation, non-vegetation and water pixels.

maskParam = list(
  list(bandName = 'MASK', minArea = 0L, bufferSize = 0L, invalidValues = c(0L:3L, 7L:11L), bufferedValues = integer())
)

The NDVI is computed by (B04 - B08) / (B04 + B08). We want it in the native Sentinel 2 resolution of 10m. Taking into account corner cases (division by 0) the indicator step config is:

indicatorIndicators = list(
  list(bandName = 'NDVI',  resolution = 10, mask = 'MASK', factor = 10000, bands = c('A' = 'B04', 'B' = 'B08'), equation = '(A.astype(float) - B) / (0.0000001 + A + B)')
)

Other steps are straightforward to configure - we just need to pick up band names and put them into right places (and choose the centiles to be computed, e.g. aggregateQuantiles = c(0.05, 0.5, 0.98)).

A common config file for all scripts involved should look as follows (let's assume it's saved in myConfig.R):

cubeRpath = 'ADJUST'
gridFile = 'ADJUST'
tmpDir = 'ADJUST'
rawDir = 'ADJUST'
periodsDir = 'ADJUST'
tilesDir = 'ADJUST'
cacheTmpl = 'ADJUST/{region}_{dateFrom}_{dateTo}_{cloudCovMax}_{bands}'

bands = c('B04', 'B08', 'SCL')
cloudCov = ADJUST
nCores = ADJUST
chunksPerCore = 10

dwnldMethod = 'download'
dwnldMaxRemovals = 1000
dwnldNCores = 4
dwnldTimeout = 120
dwnldSkipExisting = 'samesize'
dwnldTries = 3

maskSkipExisting = TRUE
maskParam = list(
  list(bandName = 'MASK', minArea = 0L, bufferSize = 0L, invalidValues = c(0L:3L, 7L:11L), bufferedValues = integer())
)

indicatorSkipExisting = TRUE
indicatorIndicators = list(
  list(bandName = 'NDVI',  resolution = 10, mask = 'MASK', factor = 10000, bands = c('A' = 'B04', 'B' = 'B08'), equation = '(A.astype(float) - B) / (0.0000001 + A + B)')
)

whichSkipExisting = TRUE
whichBands = c('NDVI')
whichPrefix = 'NMAX'
whichDoyPrefix = 'DOYMAX'
whichBlockSize = 512

aggregateSkipExisting = TRUE
aggregateBands = c('NDVI2')
aggregateBlockSize = 512
aggregateQuantiles = c(0.05, 0.5, 0.98)
aggregateCounts = FALSE
aggregateCountsBand = 'NDVI2'
aggregateCountsOutBand = 'N2'

tileSkipExisting = TRUE
tileRawBands = character()
tilePeriodBands = list(
  '1 month' = c('NDVI'),
)
tileResamplingMethod = 'near'
tileGdalOpts = '--config GDAL_CACHEMAX 1024 -co "COMPRESS=DEFLATE" -co "TILED=YES" -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512"'

Workflow

Let's assume the processing is done for the region myRegion for the whole year 2018.

Rscript      init.R myConfig.R myRegion 2018-01-01 2018-13-31 {yourApiLogin} {yourApiPswd}
Rscript     dwnld.R myConfig.R myRegion 2018-01-01 2018-13-31
Rscript      mask.R myConfig.R myRegion 2018-01-01 2018-13-31
Rscript indicator.R myConfig.R myRegion 2018-01-01 2018-13-31
Rscript     which.R myConfig.R myRegion 2018-01-01 2018-13-31 '1 year'
Rscript aggregate.R myConfig.R myRegion 2018-01-01 2018-13-31 '1 year'
Rscript      tile.R myConfig.R myRegion 2018-01-01 2018-13-31

Remarks:

  • The dwnld.R, mask.R and indicator.R scripts can be run regularly, e.g. daily or weekly to spread the workload over a longer period of time.
  • The which.R, aggregate.R and tile.R scripts can be run yearly.