更新时间:2021-07-03 00:31:33
这里是使用光栅"包的另一种方法.该功能在空间上聚合要切碎的栅格,将聚合的栅格像元转换为多边形,然后使用每个多边形的范围来裁剪输入栅格.
Here another approach using "raster" package. The function spatially aggregates the raster to be chopped, the aggregated raster cells are turned into polygons, then each polygon's extent is used to crop the input raster.
我敢肯定,有很多复杂的方法可以做到这一点,但是这种方法对我有用,我也觉得它很直观.我希望您也觉得它有用.注意第4部分和下面的5仅用于测试,它们不是该功能的一部分.
I am sure there are sophisticated and compact ways to do this but this approach works for me and I found it intuitive as well. I hope you find it useful too. Notice Part 4 & 5 below are only for testing and they are not part of the function.
logo <- raster(system.file("external/rlogo.grd", package="raster"))
plot(logo,axes=F,legend=F,bty="n",box=FALSE)
# The function spatially aggregates the original raster
# it turns each aggregated cell into a polygon
# then the extent of each polygon is used to crop
# the original raster.
# The function returns a list with all the pieces
# in case you want to keep them in the memory.
# it saves and plots each piece
# The arguments are:
# raster = raster to be chopped (raster object)
# ppside = pieces per side (integer)
# save = write raster (TRUE or FALSE)
# plot = do you want to plot the output? (TRUE or FALSE)
SplitRas <- function(raster,ppside,save,plot){
h <- ceiling(ncol(raster)/ppside)
v <- ceiling(nrow(raster)/ppside)
agg <- aggregate(raster,fact=c(h,v))
agg[] <- 1:ncell(agg)
agg_poly <- rasterToPolygons(agg)
names(agg_poly) <- "polis"
r_list <- list()
for(i in 1:ncell(agg)){
e1 <- extent(agg_poly[agg_poly$polis==i,])
r_list[[i]] <- crop(raster,e1)
}
if(save==T){
for(i in 1:length(r_list)){
writeRaster(r_list[[i]],filename=paste("SplitRas",i,sep=""),
format="GTiff",datatype="FLT4S",overwrite=TRUE)
}
}
if(plot==T){
par(mfrow=c(ppside,ppside))
for(i in 1:length(r_list)){
plot(r_list[[i]],axes=F,legend=F,bty="n",box=FALSE)
}
}
return(r_list)
}
SplitRas(raster=logo,ppside=3,save=TRUE,plot=TRUE)
# in this example we chopped the raster in 3 pieces per side
# so 9 pieces in total
# now the raster pieces should be ready
# to be processed in the default directory
# A feature I like about this function is that it plots
# the pieces in the original order.
# notice if you cropped a rasterbrick
# use "brick" instead of "raster" to read
# the piece back in R
list2 <- list()
for(i in 1:9){ # change this 9 depending on your number of pieces
rx <- raster(paste("SplitRas",i,".tif",sep=""))
# piece_processed <- HERE YOU RUN YOUR CODE
writeRaster(piece_processed,filename=paste("SplitRas",i,sep=""),
format="GTiff",datatype="FLT4S",overwrite=TRUE)
}
# once a code has been ran on those pieces
# we save them back in the directory
# with the same name for convenience
# read each piece back in R
list2 <- list()
for(i in 1:9){ # change this 9 depending on your number of pieces
rx <- raster(paste("SplitRas",i,".tif",sep=""))
list2[[i]] <- rx
}
# mosaic them, plot mosaic & save output
list2$fun <- max
rast.mosaic <- do.call(mosaic,list2)
plot(rast.mosaic,axes=F,legend=F,bty="n",box=FALSE)
writeRaster(rast.mosaic,filename=paste("Mosaicked_ras",sep=""),
format="GTiff",datatype="FLT4S",overwrite=TRUE)