[R-sig-Geo] Get coordinates from huge shapefiles.

[Edzer Pebesma]

Mike, I am not sure, but the following toy example, which builds on the
toy data in vignette("sp"), might help:

library(sp)
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)
Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3, Sr4), "s3/4")
SpP = SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)
spdf = SpatialPolygonsDataFrame(SpP, data.frame(a=1:3), match.ID=F)
as(spdf, "SpatialLinesDataFrame")
as(as(spdf, "SpatialLinesDataFrame"), "SpatialPointsDataFrame")
as.data.frame(as(as(spdf, "SpatialLinesDataFrame"),
"SpatialPointsDataFrame"))


On 03/01/2012 01:40 AM, Christian Jansson wrote:
> Hi Mike.
> 
> Thanks a lot! It helped one step further.
> 
> I understand that this: "lapply(p at polygons, function(x)
> dim(x at Polygons[[1]]@coords))"
> gives me the number of points per polygon, and that xy is a table of all
> coordinates. So if for example the first value in the lapply above is 14,
> then the first 14 rows of xy belong to polygon 1. How can I efficiently add
> a third column to add the polygon number per coordinate pair?
> 
> Best Regards,
> /Chris
> 
> On Wed, Feb 29, 2012 at 5:44 PM, Michael Sumner <mdsumner at gmail.com> wrote:
> 
>> Sorry, that should have been "x at y" for the raw access to S4 elements.
>>
>> On Thu, Mar 1, 2012 at 7:44 AM, Michael Sumner <mdsumner at gmail.com> wrote:
>>> You can do this:
>>>
>>> https://stat.ethz.ch/pipermail/r-sig-geo/2010-February/007619.html
>>>
>>> (The email archive, or my browser has auto-converted "x @ y" symbols
>>> to "x at y" rather helpfully so you will need to fix that if that's
>>> what you see).
>>>
>>> The fortify tools (in ggplot2) and the geometry/topology-sense in
>>> rgeos are important to consider though, fortify is trying to build a
>>> table that still stores all the topology so that (for e.g.) ggplot2
>>> can work with it in a standard way.
>>>
>>> Raw extraction does not take any care with whether the polygons make
>>> sense, but that takes work since shapefiles just don't store that
>>> information, they are just sets of ring boundaries and it's up to the
>>> software using them to build the actual topology. My code at that post
>>> just drops the final duplicated coordinate . . .
>>>
>>> Cheers, Mike.
>>>
>>> On Thu, Mar 1, 2012 at 5:15 AM, Christian Jansson <chrjan70 at gmail.com>
>> wrote:
>>>> Hi,
>>>>
>>>> Thanks to everyone on this list. It is a true source of help!
>>>>
>>>> I'm reading huge shapefiles like this:
>>>>
>>>> shp_0 <- rgdal::readOGR(dsn="dir", layer="TheLayer")
>>>>
>>>> There are two columns in the shp_0 at data: bin1 and bin2.
>>>>
>>>> I need to get a table with three columns: id, X, Y.
>>>>
>>>> First I do this:
>>>> shp_0 at data$id = rownames(shp_0 at data)
>>>>
>>>> When I then run this:
>>>> tbl_coords <-fortify.SpatialPolygonsDataFrame(shp_0)
>>>> I get: "Using bin1 to define regions.
>>>> Error in createPolygonsComment(p) :
>>>>  rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon
>>>> for hole at index 3"
>>>>
>>>> After reading in previous messages on this email-list, I tried:
>>>>
>>>> slot(shp_0, "polygons") <- lapply(slot(shp_0, "polygons"),
>>>> checkPolygonsHoles), which proceeds to a new R-prompt with no issue
>>>> mentioned.
>>>>
>>>> But everytime I then run
>>>> shp_1 <- unionSpatialPolygons(shp_0, as.character(shp_0$id))
>>>> I get another crasch:
>>>> "This application has requested the Runtime to terminate it in an
>> unusual
>>>> way.
>>>> Please contact the application's support team for more information."
>>>>
>>>> How can I get the coordinates per polygon in one large table? I know
>> how to
>>>> get them per polygon, but to rbind more than 30 000 polygons take
>> forever.
>>>> To get the areas is easy "sapply(shp_0 at polygons, function(x) x at area)",
>> but
>>>> the coords are deeper into the sp-object and I don't know how to
>> "apply-it
>>>> out".
>>>>
>>>> Thank you very much for all help.
>>>>
>>>> /Chris
>>>>
>>>>        [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>
>>>
>>> --
>>> Michael Sumner
>>> Institute for Marine and Antarctic Studies, University of Tasmania
>>> Hobart, Australia
>>> e-mail: mdsumner at gmail.com
>>
>>
>>
>> --
>> Michael Sumner
>> Institute for Marine and Antarctic Studies, University of Tasmania
>> Hobart, Australia
>> e-mail: mdsumner at gmail.com
>>
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de