Nothing
#+ getLandsatTilesFunction
getLandsatTiles = function(x,
tiles = getLandsatCoverage(),
buffer = 1000,
crsOut = raster::projection(x)) {
xExtent =
raster(extend(extent(x), buffer),
ncol=10, nrow=10, crs=projection(x))
xExtentProj = projectExtent(xExtent, projection(tiles))
tilesCrop = raster::crop(
tiles, xExtentProj
)
tilesCrop2 = spTransform(tilesCrop, projection(xExtent))
tilesCrop3 = raster::crop(tilesCrop2, xExtent)
tiles = tiles[tiles$PR %in% tilesCrop3$PR, ]
if(identical(crsOut, 'landsat'))
crsOut = paste(
'+proj=lsat +path=',
tiles$PATH[1],
' +lsat=5 +ellps=WGS84 +datum=NAD83', sep='')
if(identical(crsOut, 'omerc'))
crsOut = mapmisc::omerc(
tiles,
angle=seq(-25,25,by=0.1)
)
tilesProjCrop = spTransform(tiles, crsOut)
tilesProjCrop$pPath = paste('p', sprintf("%03d", tilesProjCrop$PATH), sep='')
tilesProjCrop$rRow = paste('r', sprintf("%03d", tilesProjCrop$ROW), sep='')
tilesProjCrop$pathString = sprintf("%03d", tilesProjCrop$PATH)
tilesProjCrop$rowString = sprintf("%03d", tilesProjCrop$ROW)
tilesProjCrop$tileString = paste(
tilesProjCrop$pathString,
tilesProjCrop$rowString, sep='')
tilesProjCrop$tTile = paste(
tilesProjCrop$pPath,
tilesProjCrop$rRow, sep=' ')
attributes(tilesProjCrop)$landsat = xExtent
tilesProjCrop
}
getLandsatScenes = function(
timeRange,
tiles,
sceneList = getLandsatSceneList()
) {
if(any(names(tiles) == 'tileString')) {
tiles = tiles$tileString
}
sceneList$tile = sprintf("%03d%03d", sceneList$path, sceneList$row)
toKeep = intersect(
which(sceneList$tile %in% tiles),
which(sceneList$acquisitionDate >= min(timeRange) &
sceneList$acquisitionDate <= max(timeRange))
)
sceneList[toKeep,]
}
getLandsatSceneList = function(
path = gsub("[[:alnum:]]+$", "landsat", tempdir()),
url = "http://landsat-pds.s3.amazonaws.com/scene_list.gz") {
localFile = file.path(path, basename(url))
outFile = gsub("[.]gz$", "", localFile, ignore.case = TRUE)
rdataFile = file.path(path, 'landsatScenes.RData')
if(file.exists(rdataFile)) {
rdataEnv = new.env()
load(rdataFile, envir=rdataEnv)
if('landsatScenes' %in% ls(envir=rdataEnv)) {
return(get('landsatScenes', envir=rdataEnv))
}
}
if(!any(file.size(localFile)>0, na.rm=TRUE)) {
download.file(url, localFile, mode='wb')
}
landsatScenes = read.table(
gzfile(localFile), sep=',', stringsAsFactors=FALSE,
header=TRUE
)
save(landsatScenes, file=rdataFile, compress='xz')
landsatScenes
}
getLandsatCoverage = function(
path = c(
options()$mapmiscCachePath,
tempdir())[1],
url = 'https://landsat.usgs.gov/sites/default/files/documents/wrs2_descending.zip'
) {
localFile = file.path(path, basename(url))
rdataFile = file.path(path, 'landsatCover.RData')
if(file.exists(rdataFile)) {
rdataEnv = new.env()
load(rdataFile, envir=rdataEnv)
if('landsatCover' %in% ls(envir=rdataEnv)) {
return(get('landsatCover', envir=rdataEnv))
}
}
if(!file.exists(localFile)) {
download.file(url, localFile, mode='wb')
}
sfile = file.path(path,
grep("shp$", unzip(localFile, list=TRUE)$Name, value=TRUE)
)
if(!file.exists(sfile))
utils::unzip(localFile, exdir=path)
landsatCoverLL = raster::shapefile(sfile)
if(requireNamespace('rgdal')) {
landsatCover = spTransform(
raster::crop(landsatCoverLL, extent(-180,180,-88,88)),
# "+proj=eck2 +ellps=WGS84 +datum=NAD83"
"+proj=sinu +ellps=WGS84 +datum=NAD83"
)
} else {
landsatCover = landsatCoverLL
}
save(landsatCover, file=rdataFile, compress='xz')
landsatCover
}
downloadLandsatTiles = function(
scenes, band = 1:11,
path = c(
options()$mapmiscCachePath,
tempdir())[1],
mode = 'wb',
...
) {
scenesLong = scenes[
rep(1:nrow(scenes), length(band)),
]
scenesLong$band = rep(band, rep(nrow(scenes), length(band)))
scenesLong$url = paste(gsub("index[.]html$", "",
scenesLong$download_url),
scenesLong$entityId, '_B', scenesLong$band, '.TIF', sep='')
scenesLong$localFile = file.path(path, basename(scenesLong$url))
if(!dir.exists(path)) dir.create(path)
for(D in 1:nrow(scenesLong)) {
if(!any(file.size(scenesLong[D, 'localFile'])>1, na.rm=TRUE) )
download.file(scenesLong[D,'url'], scenesLong[D, 'localFile'],
mode=mode, ...)
}
result = tapply(
X=scenesLong$localFile,
INDEX=scenesLong[,'entityId', drop=FALSE],
FUN = function(X)
raster::stack(sort(X)),
simplify=FALSE
)
lapply(result, function(x) x)
}
landsatNdvi = function(x,
suffix = 1:3,
path = c(
options()$mapmiscCachePath,
tempdir())[1],
filename = 'ndvi.grd') {
if(basename(filename) == filename) {
filename = file.path(path, filename)
}
suffixString = paste("_(", paste(suffix, collapse='|'), ')$', sep='')
StilesAll = gsub(suffixString, '', names(x))
Stiles = unique(StilesAll)
tileMat = list()
for(D in suffix)
tileMat[[D]] = outer(
names(x),
paste(Stiles, '_', D, sep=''),
'==')
numerMat = tileMat[[suffix[2]]] - tileMat[[suffix[1]]]
denomMat = tileMat[[suffix[2]]] + tileMat[[suffix[2]]]
qMat = tileMat[[suffix[3]]]
colnames(numerMat) = colnames(denomMat) = colnames(qMat) =
Stiles
calc(
x,
function(xx) {
numer = xx %*% numerMat
denom = xx %*% denomMat
qual = xx %*% qMat
qual[qual > 25000 | is.na(numer) | is.na(denom) |
denom == 0] = NA
result = numer/denom
result[is.na(qual)] = NA
meanResult = mean(result, na.rm=TRUE)
meanResult[is.nan(meanResult)] = NA
result = c(result, meanResult, sum(!is.na(result)))
names(result) = c(colnames(numerMat), 'mean', 'N')
result
},
filename = filename,
overwrite = file.exists(filename)
)
}
#'
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.