Converting GRID to SHAPE in R (was: Re: [Fwd: Re: rgdal: problem with grid data])

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

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

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

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

Hi Roger and the others,

I've noticed that readRAST6() uses closeAllConnection(), as it is trying
to close a raw "socket" type connection I had opened for parallel
processing (snow package).

Isn't this a side effect and perhaps closeAllConnection() should be
replaced here by any other alternative?

Thanks (also for the package) and best regards,

Javier

---
Javier Garc?a-Pintado
75 Trotsworth Court
Christchurh Road
Virginia Water
GU25 4AH
Surrey
United Kingdom

On Sun, 28 Mar 2010, jgarcia at ija.csic.es wrote:

Hi Roger and the others,

I've noticed that readRAST6() uses closeAllConnection(), as it is trying
to close a raw "socket" type connection I had opened for parallel
processing (snow package).

Isn't this a side effect and perhaps closeAllConnection() should be
replaced here by any other alternative?

Thanks (also for the package) and best regards,

Javier,

Could you try a change in the package source from sourceforge, check out
the CVS source from:

cvs -z3 \
-d:pserver:anonymous at r-spatial.cvs.sourceforge.net:/cvsroot/r-spatial co \
-P spgrass6

I've added a flag to permit the call to closeAllConnections() to be
avoided; you would use close_OK=FALSE. For most users it is still needed,
as the number of connections opened in trying to guess what the GDAL GRASS
plugin needs is not well-specified.

Hope this helps,

Roger

Javier

---
Javier Garc?a-Pintado
75 Trotsworth Court
Christchurh Road
Virginia Water
GU25 4AH
Surrey
United Kingdom

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no

On Sun, 28 Mar 2010, jgarcia at ija.csic.es wrote:

Hi Roger and the others,

I've noticed that readRAST6() uses closeAllConnection(), as it is trying
to close a raw "socket" type connection I had opened for parallel
processing (snow package).

Isn't this a side effect and perhaps closeAllConnection() should be
replaced here by any other alternative?

Thanks (also for the package) and best regards,

Javier,

Could you try a change in the package source from sourceforge, check out
the CVS source from:

cvs -z3 \
-d:pserver:anonymous at r-spatial.cvs.sourceforge.net:/cvsroot/r-spatial co \
-P spgrass6

I've added a flag to permit the call to closeAllConnections() to be
avoided; you would use close_OK=FALSE. For most users it is still needed,
as the number of connections opened in trying to guess what the GDAL GRASS
plugin needs is not well-specified.

Hope this helps,

Roger

Javier

---
Javier Garc?a-Pintado
75 Trotsworth Court
Christchurh Road
Virginia Water
GU25 4AH
Surrey
United Kingdom

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no

Hi Roger,
Just to confirm that everything goes OK.

Thanks again,
Javier
---

On Sun, 28 Mar 2010, jgarcia at ija.csic.es wrote:

Hi Roger and the others,

I've noticed that readRAST6() uses closeAllConnection(), as it is trying
to close a raw "socket" type connection I had opened for parallel
processing (snow package).

Isn't this a side effect and perhaps closeAllConnection() should be
replaced here by any other alternative?

Thanks (also for the package) and best regards,

Javier,

Could you try a change in the package source from sourceforge, check out
the CVS source from:

cvs -z3 \
-d:pserver:anonymous at r-spatial.cvs.sourceforge.net:/cvsroot/r-spatial co \
-P spgrass6

I've added a flag to permit the call to closeAllConnections() to be
avoided; you would use close_OK=FALSE. For most users it is still needed,
as the number of connections opened in trying to guess what the GDAL GRASS
plugin needs is not well-specified.

Hope this helps,

Roger

Javier

---
Javier Garc?a-Pintado
75 Trotsworth Court
Christchurh Road
Virginia Water
GU25 4AH
Surrey
United Kingdom

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no