[Fwd: Re: raster to polygons]
-------- Original Message -------- Subject: Re: raster to polygons Date: Thu, 25 Mar 2010 18:08:44 +0100 From: Agustin Lobo <Agustin.Lobo at ija.csic.es> Reply-To: Agustin.Lobo at ija.csic.es To: Tim (LWF) <Tim.Haering at lwf.bayern.de>, sig-geo <r-sig-geo at stat.math.ethz.ch> References: <70FC67C4A585D1489E66225A4E0238BA7C8B29 at RZS-EXC-CL06.rz-sued.bayern.de> <4BAA436A.2040106 at ija.csic.es> <70FC67C4A585D1489E66225A4E0238BA7C90E3 at RZS-EXC-CL06.rz-sued.bayern.de> Tim, 1. You do not need the first conversion:
comm1 <- paste("gdal_translate ",paste(dirin,"/",infname,".tif",sep=""),
"delme.tif")
system(comm1,inter=T)
as this is transforming to tif and you already start from a tif. 2. If your input is not a grid image, then you do not need all the part starting with read.dbf() 3. Both gdal_translate and gdal_polygonize.py are part of gdal utilities (not just the gdal library that you must have if you use rgdal) and my function uses system() to execute them as external commands in the OS (linux ubuntu in my case). The good thing of gdal_polygonize.py is that it is very fast and memory efficient. But it looks like you do not have them, and I'm not sure how cumbersome installing this will be under windows now, it used to be tricky. Check in gdal.com, and/or other people from this lists might advice you, but perhaps you can carry out the same type of operation with grass or saga through R. My experience in stalling packages like gdal or grass in windows through osgeo4w is that it implies a time that is better invested on using linux. Perhaps installing saga for windows is easier, I have only tried it on linux systems. 4. Finally, check package raster, perhaps Robert Hijmans has included a function equivalent to gdal_polygonize.py That would be the very best in my opinion. Agus
H?ring wrote:
Hi Agus,
Yes, I get your message. Thanks for that! Unfortunately I was not successful applying your function to my raster. My input file is in tif-format. Hopefully I changed your code correctly... Here is what I did:
# I have a xyz-file "xyzdat":
"migrid2shp" <-
function(dirin="D:/temp",
infname=" gridfile ",
dirout="D:/temp",
ofname="grid_out")
{
require(rgdal)
require(foreign)
comm1 <- paste("gdal_translate ",paste(dirin,"/",infname,".tif",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,".tif",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=""))
}
migrid2shp()
R gives me following error message:
Error in system(comm1, inter = T) : gdal_translate not found
I`m using R 2.10.1 and rgdal based on GDAL 1.7.1 under Windows XP.
Any ideas?
Thanks a lot
TIM
Dr. Agustin Lobo Institut de Ciencies de la Terra "Jaume Almera" (CSIC) LLuis Sole Sabaris s/n 08028 Barcelona Spain Tel. 34 934095410 Fax. 34 934110012 email: Agustin.Lobo at ija.csic.es http://www.ija.csic.es/gt/obster