Split RGB raster in multiples pieces and save file in GeoTiff with conditions

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Split RGB raster in multiples pieces and save file in GeoTiff with conditions

ASANTOS
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)

         [hidden email]
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/

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo