[R-sig-Geo] Converting GRID to SHAPE in R

Agustin Lobo alobolistas at gmail.com
Wed Mar 24 11:09:58 CET 2010


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