Jônatan, excelente!

no seu script apenas substitui rasterToPoints para rasterToPolygons no final.

O objetivo do raster é apenas para extrair as informações do índice (i, j). Os valores do pixel funcionam como ID's.

Agora, o resultado final já é georreferenciado e as quadrículas estão prontas para uma posterior análise geoestatística.

Att,

Alex Santos

### ENTRADA DE DADOS ###

r <- raster(nc=20,nr=10)

r[] <- rep(1:5,each=40)

plot(r)

res <- data.frame(rowColFromCell(r,1:ncell(r)), 1:ncell(r),getValues(r))

names(res) <- c("row","col","cell","v")

## criando stack com 4 layers

s <- stack(r,r,r,r)

## atribuindo valores para cada layer baseado no res

s[] <- as.matrix(res)

plot(s)

#p<-rasterToPoints(s, spatial=T)

pp<- rasterToPolygons(s,dissolve=F)

#spplot(pp)

# writeOGR(p,getwd(),"p2",driver="ESRI Shapefile")

writeOGR(pp,getwd(),"matriz_pl",driver="ESRI Shapefile")

pescrito <- readOGR(getwd(), "matriz_pl")

pescrito@data