[R-sig-Geo] [Fwd: Re: raster to polygons]

[Agustin Lobo]

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