[R-sig-Geo] Split RGB raster in multiples pieces and save file in GeoTiff with conditions
ASANTOS
@|ex@ndre@@nto@br @end|ng |rom y@hoo@com@br
Mon Jul 8 00:21:10 CEST 2019
Dear R-Sig-Geo Members,
?????? I've like to split a RGB raster in multiple pieces and save each
image in GeoTiff format, but in the file name I need the count of points
(pts) inside each piece. My code runs perfectly when I used a single
raster image, but with stack or brick object doesn't work. I've like as
final output in my directory nine GeoTiff files (I divided the original
raster in 9 parts) with files names: SplitRas1points0.tif ...
SplitRas9points0.tif
?????? My code is:
#Packages
library(raster)
library(sp)
library(rgdal)
## Create a simulated RBG rasters
r <- raster(nc=30, nr=30)
r <- setValues(r, round(runif(ncell(r))* 255))
g <- raster(nc=30, nr=30)
g <- setValues(r, round(runif(ncell(r))* 255))
b <- raster(nc=30, nr=30)
b <- setValues(r, round(runif(ncell(r))* 255))
rgb<-stack(r,g,b)
plotRGB(rgb,
?????????????? r = 1, g = 2, b = 3)
##Given interesting points coordinates
xd???????? <- c(-24.99270,45.12069,99.40321,73.64419)
yd?? <- c(-45.435267,-88.369745,-7.086949,44.174530)
pts <- data.frame(xd,yd)
pts_s<- SpatialPoints(pts)
points(pts_s, col="black", pch=16)
# This 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
# 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)
}
#Slip RGB raster in 3 parts
splitRBG<-SplitRas(raster=rgb,ppside=3,save=FALSE,plot=FALSE)
#Count number of points inside each piece of raster
res<-NULL
for(i in 1:3){
?? pointcount = function(r, pts){
?? # make a raster of zeroes like the input
?? r2 = r
?? r2[] = 0
?? # get the cell index for each point and make a table:
?? counts = table(cellFromXY(r,pts))
?? # fill in the raster with the counts from the cell index:
?? r2[as.numeric(names(counts))] = counts
?? return(r2)
??}
r2 = pointcount(splitRBG[i], pts_s)
res<-rbind(res,c(r2))
}
#
#Run a code on each piece with the number of points inside each piece in
the file name (res[i])
list2 <- list()
for(i in 1:3){ # 3 is the number of pieces
?? rx <- raster(paste("splitRBG",i,".tif",sep=""))
writeRaster(splitRBG,filename=paste("SplitRas",i,"points",res[i],sep=""),
?????????????????????????? format="GTiff",datatype="FLT4S",overwrite=TRUE)
}
#
Any ideas please?
Thanks in advanced,
Alexandre
--
======================================================================
Alexandre dos Santos
Prote????o Florestal
IFMT - Instituto Federal de Educa????o, Ci??ncia e Tecnologia de Mato Grosso
Campus C??ceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
C??ceres - MT CEP: 78.200-000
Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)
alexandre.santos using cas.ifmt.edu.br
Lattes: http://lattes.cnpq.br/1360403201088680
OrcID: orcid.org/0000-0001-8232-6722
Researchgate: www.researchgate.net/profile/Alexandre_Santos10
LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635
Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/
More information about the R-sig-Geo
mailing list