# rank change log
# 20180524 - add TIGER "S1500" (vehicle trails) to Roads class 2
# 20180524 - add argument for TIGER data year, even state lines move
# so e sure sure years match when merging/stacking rasters
buildFeatureStack <- function(baseLayer,mapshape,
spList,
filterList=NULL,
maxRasterize=10000,
sliceFeatureBuffer=0,
polySimplify=0,polyMethod="vis",
polyWeighting=0.85,polySnapInt=0.0001,
silent=FALSE,noisy=FALSE) {
# baseLayer a rasterLayer
# returns a rasterStack which has 4 layers which can be stored as ints
# collapse this down to spList and filterList
bLrect <- as(raster::extent(baseLayer), "SpatialPolygons")
sp::proj4string(bLrect) <- raster::crs(baseLayer)
spList <- filterFeatures(spList,filterList)
featureNames <- c("spTown","spRoads","spWaterA","spWaterL")
for (x in featureNames) {
if( x %in% names(spList) &&
length(spList[[x]])>0 ) {
assign(x,spList[[x]])
} else {
assign(x,NULL)
}
}
if (length(spRoads)>0) {
if (!silent) print("roads")
spRoads <- sxdfMask(spRoads,bLrect)
}
if (length(spRoads)>0) {
if (!silent) print(paste0(nrow(spRoads)," roads to process"))
rlayer <- shapeToRasterLayer(sxdf=spRoads,
templateRaster=baseLayer,
maxRasterize=maxRasterize,
silent=silent,noisy=noisy)
if (!is.finite(rlayer@data@min)) {
warning("rlayer mess-up")
print(rlayer@data@min)
}
} else {
rlayer <- raster::raster(nrows=nrow(baseLayer),
ncols=ncol(baseLayer),
ext=extent(baseLayer),
crs=crs(baseLayer),
vals=0)
if (!silent) print("no roads to add")
}
if (length(spWaterL)>0) {
if (!silent) print("water lines")
spWaterL <- sxdfMask(spWaterL,bLrect)
}
if (length(spWaterL)>0) {
if (!silent) print(paste0(nrow(spWaterL)," water lines to process"))
wLlayer <- shapeToRasterLayer(sxdf=spWaterL,
templateRaster=baseLayer,
maxRasterize=maxRasterize,
silent=silent,noisy=noisy)
if (!is.finite(wLlayer@data@min)) {
warning("wLlayer mess-up")
print(wLlayer@data@min)
}
} else {
wLlayer <- raster::raster(nrows=nrow(baseLayer),
ncols=ncol(baseLayer),
ext=extent(baseLayer),
crs=crs(baseLayer),
vals=0)
if (!silent) print("no water lines to add")
}
if (length(spWaterA)>0) {
if (!silent) print("water polygons")
if (sliceFeatureBuffer>0) {
bLrectA <- bufferUnion(bLrect,
mapbuffer=sliceFeatureBuffer,
mapunion=NULL)
} else {
bLrectA <- bLrect
}
spWaterA <- sxdfMask(spWaterA,bLrectA)
}
if (length(spWaterA)>0) {
if (!silent) print(paste0(nrow(spWaterA)," water areas to process"))
wAlayer <- shapeToRasterLayer(sxdf=spWaterA,
templateRaster=baseLayer,
maxRasterize=maxRasterize,
polySimplify=polySimplify,
polyMethod=polyMethod,
polyWeighting=polyWeighting,
polySnapInt=polySnapInt,
silent=silent,noisy=noisy)
if (!is.finite(wAlayer@data@min)) {
warning("wAlayer mess-up")
print(wAlayer@data@min)
}
} else {
wAlayer <- raster::raster(nrows=nrow(baseLayer),
ncols=ncol(baseLayer),
ext=extent(baseLayer),
crs=crs(baseLayer),
vals=0)
if (!silent) print("no water polygons to add")
}
if (length(spTown)>0) {
if (!silent) print("towns")
if (sliceFeatureBuffer>0) {
bLrectA <- bufferUnion(bLrect,
mapbuffer=sliceFeatureBuffer,
mapunion=NULL)
} else {
bLrectA <- bLrect
}
spTown <- sxdfMask(spTown, bLrectA)
}
if (length(spTown)>0) {
if (!silent) print(paste0(nrow(spTown)," towns to process"))
tlayer <- shapeToRasterLayer(sxdf=spTown,
templateRaster=baseLayer,
maxRasterize=maxRasterize,
polySimplify=polySimplify,
polyMethod=polyMethod,
polyWeighting=polyWeighting,
polySnapInt=polySnapInt,
silent=silent,noisy=noisy)
if (!is.finite(tlayer@data@min)) {
warning("tlayer mess-up")
print(tlayer@data@min)
}
} else {
tlayer <- raster::raster(nrows=nrow(baseLayer),
ncols=ncol(baseLayer),
ext=extent(baseLayer),
crs=crs(baseLayer),
vals=0)
if (!silent) print("no towns to add")
}
s <- raster::stack(tlayer,wAlayer,wLlayer,rlayer)
names(s) <- c("town","waterA","waterL","roads")
return(s)
}
shapeToRasterLayer <- function(sxdf,templateRaster,
maxRasterize=10000,
polySimplify=0,polyMethod="vis",
polyWeighting=0.85,polySnapInt=0.0001,
keepTouch=FALSE,
silent=FALSE,noisy=FALSE) {
# return a rasterLayer based on templateRaster, rasterizing sxdf lines/polygons
# using values in sxdf@data[,"rank"]
zeroRaster <- templateRaster
raster::values(zeroRaster) <- 0
retRaster <- zeroRaster
CP <- as(extent(templateRaster), "SpatialPolygons")
sp::proj4string(CP) <- raster::crs(templateRaster)
sxdf <- sxdfMask(sxdf,CP,keepTouch=keepTouch)
if (!is.null(sxdf)) {
sxdf <- sxdf[order(sxdf$value),] # sort for rasterize
if (class(sxdf)=="SpatialPolygonsDataFrame") {
if (polySimplify>0) {
sxdf <- rmapshaper::ms_simplify(sxdf,keep=polySimplify,
method=polyMethod,
weighting=polyWeighting,
snap=TRUE,snap_interval=polySnapInt)
sxdf <- sp::spTransform(sxdf,raster::crs(templateRaster))
}
}
nloops <- ceiling(nrow(sxdf)/maxRasterize)
for (i in 1:nloops) {
if (nloops > 1) print(paste0(i," / ",nloops))
firstrow <- (i-1)*maxRasterize + 1
lastrow <- min(i*maxRasterize,nrow(sxdf))
gc() # cleanup, this takes a lot of memory
#rlayer <- velox::velox(zeroRaster)
rlayer <- zeroRaster
temptime <- round(system.time(
#rlayer$rasterize(spdf=sxdf[firstrow:lastrow,],field="value",background=0)
rlayer <- rasterize(sxdf[firstrow:lastrow,],rlayer,field="value",background=0)
)[[3]],digits=2)
if (!silent) print(paste0(" ",temptime," rasterizing"))
if (nloops>1) {
temptime <- round(system.time(
#retRaster <- maxLayerValue(retRaster,rlayer$as.RasterLayer(band=1))
retRaster <- maxLayerValue(retRaster,rlayer)
)[[3]],digits=2)
if (!silent) print(paste0(" ",temptime," combining"))
} else {
#retRaster <- rlayer$as.RasterLayer(band=1)
retRaster <- rlayer
}
}
}
return(retRaster)
}
maxLayerValue <- function(rasterLayer1,rasterLayer2) {
return(max(raster::stack(list(rasterLayer1,rasterLayer2))))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.