[R-sig-Geo] Converting GRID to SHAPE in R
Agustin Lobo
alobolistas at gmail.com
Wed Mar 24 14:39:48 CET 2010
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:
> One further problem: it seems that readGDAL()
> does not read the info from the grd file, just the gri.
> Thus the legend is ignored.
>
> Agus
>
> Agustin Lobo wrote:
>> Thanks for all the inputs I've got. I've
>> tried the following solution, but the steps
>> with as.SpatialPolygons.GridTopology() and
>> unionSpatialPolygons() take too long a time, had to abort the
>> second one after 1h, so
>> the procedure does not seem to be a really working
>> solution (note than I'm using the landcover map of a tiny country as
>> example)
>>
>> #Data from http://www.diva-gis.org/data/cov/ISR_cov.zip
>>
>> isrlc <-
>> readGDAL("/media/Transcend/MASTER_ICTA2007_2008/GEODATA_miniprojs/Kefah/ISR_cov/isr_cov.gri")
>>
>>
>> > class(isrlc)
>> [1] "SpatialGridDataFrame"
>> attr(,"package")
>> [1] "sp"
>>
>> isrlcSP <- as.SpatialPolygons.GridTopology(getGridTopology(isrlc))
>>
>> > class(isrlcSP)
>> [1] "SpatialPolygons"
>> attr(,"package")
>> [1] "sp"
>>
>> isrlcSPU <- unionSpatialPolygons(isrlcSP, isrlc at data[,1])
>> Aborted after 1h
>>
>> A second solution is saving as gtif, then gdal_polygonize:
>>
>> gdal_translate isr_cov.gri isr_cov.tif
>> gdal_polygonize.py -f 'ESRI Shapefile' isr_cov.tif isr_cov
>>
>> which takes longer for writing the commands than the actual execution
>> time!
>> Then we still have to attach the legend, which can be easily done in
>> R. So I think the best is writing a function using system() for the
>> above 2 commands, then
>> reading the grid and the shape, attach the legend to the SPDF and save
>> again as
>> shape.
>>
>> So, problem solved. I'll post the function to the list as soon as I
>> make it.
>>
>> Agus
>>
>>
>> Michael Sumner wrote:
>>> Hi, I'm not sure on the license status of gpclib in maptools - not for
>>> commercial use now - but I'm sure Roger will comment on that on the
>>> list. This seems to work:
>>>
>>> library(maptools)
>>>
>>> gdata <- image2Grid(list(x = 1:nrow(volcano), y = 1:ncol(volcano), z
>>> = volcano))
>>>
>>> SpP_grd <- as.SpatialPolygons.GridTopology(getGridTopology(gdata))
>>>
>>>
>>> unionpoly <- unionSpatialPolygons(SpP_grd, gdata$z)
>>>
>>>
>>> trdata <- data.frame(z = sort(unique(gdata$z)), row.names =
>>> sapply(slot(unionpoly, "polygons"), function(i) slot(i, "ID")))
>>>
>>>
>>> SpPDF <- SpatialPolygonsDataFrame(unionpoly, trdata)
>>>
>>>
>>> I checked the result in GIS and it seems to work perfectly.
>>>
>>> Cheers, Mike.
>>>
>>> On Tue, Mar 23, 2010 at 6:01 PM, Agustin Lobo <alobolistas at gmail.com>
>>> wrote:
>>>>
>>>> -------- Original Message --------
>>>> Subject: Re: [R-sig-Geo] rgdal: problem with grid data
>>>> Date: Tue, 23 Mar 2010 07:53:38 +0100
>>>> From: Agustin Lobo <aloboaleu at gmail.com>
>>>> Reply-To: Agustin.Lobo at ija.csic.es
>>>> To: Michael Sumner <mdsumner at gmail.com>
>>>> CC: sig-geo <r-sig-geo at stat.math.ethz.ch>
>>>> References: <4BA86074.8010709 at gmail.com>
>>>> <522664f81003222349g39cc0c4ya756547d80c1d2f6 at mail.gmail.com>
>>>>
>>>> ok, I see, I have to convert to Spatial Polygons data frame first,
>>>> right?
>>>> Ideally, I should be able of clumping those cells with the same value,
>>>> so this is a bit more complicated than I expected. But should be
>>>> possible.
>>>> I'll try and let the list know.
>>>> Thanks
>>>>
>>>> Agus
>>>>
>>>> Michael Sumner wrote:
>>>>> What are you expecting to get - a point for every raster cell? If so,
>>>>> use gridded(isrlc) <- FALSE to convert to a SPDF first - then the
>>>>> write will work.
>>>>>
>>>>> More complicated conversions are possible, such as with
>>>>> as.SpatialPolygons.GridTopology.
>>>>>
>>>>> Cheers, Mike.
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Tue, Mar 23, 2010 at 5:32 PM, Agustin Lobo <alobolistas at gmail.com>
>>>>> wrote:
>>>>>> Hi!
>>>>>>
>>>>>> I read a grid file into R with
>>>>>> isrlc <-
>>>>>>
>>>>>> readGDAL("/media/Transcend/MASTER_ICTA2007_2008/GEODATA_miniprojs/Kefah/ISR_cov/isr_cov.gri")
>>>>>>
>>>>>>
>>>>>> While isrlc seems correct, I'm getting an empty shapefile when I
>>>>>> attempt
>>>>>> to
>>>>>> write:
>>>>>>
>>>>>>
>>>>>> writeOGR(isrlc,dsn="/media/Transcend/MASTER_ICTA2007_2008/GEODATA_miniprojs/Kefah",
>>>>>>
>>>>>> layer="ISR_cov2",driver="ESRI Shapefile")
>>>>>>
>>>>>> Data from
>>>>>> http://www.diva-gis.org/data/cov/ISR_cov.zip
>>>>>>
>>>>>> Thanks
>>>>>>
>>>>>> Agus
>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> R-sig-Geo at stat.math.ethz.ch
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>
>>>> --
>>>> 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
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at stat.math.ethz.ch
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>
>>
>
>
More information about the R-sig-Geo
mailing list