Skip to content
Prev 7903 / 29559 Next

Converting GRID to SHAPE in R

This is the function I've made. Instead of reading the shape files
(which takes a long time) I read the dbf file.
I also read the grd file,
writing back a modified dbf including the legend.
The only problem is that both read.dbf() and
write.dbf() read and write many more cases than those actually
present in the dbf file (hence the na.omit() after read.dbf() )
Roger, do you want me to send you the actual dbf file?

This is not perfect, but works and it's fast. Don't know
if this is a general solution dor grd files.

"migrid2shp" <- 
function(dirin="/media/Transcend/MASTER_ICTA2007_2008/GEODATA_miniprojs/Kefah/ISR_cov", 
infname="ISR_cov", 
dirout="/media/Transcend/MASTER_ICTA2007_2008/GEODATA_miniprojs/Kefah", 
ofname="ISR_cov")
{
require(rgdal)
require(foreign)
comm1 <- paste("gdal_translate ",paste(dirin,"/",infname,".gri",sep=""), " 
delme.tif")
system(comm1,inter=T)
comm2 <- paste("gdal_polygonize.py -f 'ESRI Shapefile' delme.tif ", 
paste(dirout,"/",ofname,sep=""), ofname)
system(comm2,inter=T)
system(paste("rm delme.tif"))

indbf <- read.dbf(paste(dirout,"/",ofname,"/",ofname, ".dbf",sep=""))
indbf <- na.omit(indbf)
a <- matrix(scan(paste(dirin,"/",infname,".grd",sep=""),what=" 
",skip=25,sep="=",fill=T),byrow=T,ncol=2)
b <- grep("Label", a[,1])
a <- a[b,]
n <- length(unlist(strsplit(a[,1],"Label")))
b <- as.numeric(unlist(strsplit(a[,1],"Label"))[seq(2,n,by=2)])
a[1:length(b),1] <-b
a <- data.frame(class=a[1:length(b),1], description=a[1:length(b),2])
row.names(a) <- a[,1]

b <- a[match(indbf$DN,a$class),]
indbf2 <- cbind(indbf,b)
write.dbf(indbf2,file=paste(dirout,"/",ofname,"/", ofname,".dbf", sep=""))
print(paste(dirout,"/",ofname,"/", ofname,".dbf", sep=""))
}
Agustin Lobo wrote: