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

[Michael Sumner]

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