[R-sig-Geo] [R] polygon shapefile from line edge coordinate list

Murray Richardson murray.richardson at utoronto.ca
Wed Mar 12 21:51:30 CET 2008


Thakns Roger

Since I am more familiar with postGIS I took that approach.  It worked 
like a charm.

I can post the script if you think it's worthwhile.

Murra

Roger Bivand wrote:
> On Sun, 9 Mar 2008, Murray Richardson wrote:
>
>> Hi Roger
>>
>> Thanks so much for your help.  I can definitely work with this.  On your
>> advice I'm responding via r-sig-geo.
>>
>> Unfortunately, there is no guarantee that the line segments will come in
>> any particular order from the alpha shapes routine.  In fact, there will
>> also likely be line segments from  different polys mixed together, so I
>> need to write the code so that the polygons are assembled only for
>> connected line segments.  I was naive when I started thinking about the
>> problem.  As you suggest, this could get messy.
>
> I understand. I might be tempted to use AWK externally to number the 
> unique (x,y) pairs, and reorder the text file externally. There is 
> some AWK code and shell scripts in the infrastructure in the maps 
> package and the two Bell Labs reports referenced there. They treat the 
> short line segment spaghetti scenario directly. Maybe that might 
> provide a way forward?
>
>>
>> As an alternative I may load the segments into a postGIS linestring
>> table and polygonize that, all via RODBC/postGIS.  An added benefit is
>> that it would clean up dangling segments, which may also exist.
>>
>
> This would be more robust, and via ogr2ogr would let you make the 
> shapefiles too. I'd prototype in AWK because I know it adequately and 
> don't know PostGIS, but if you already know PostGIS, the underlying 
> code is certainly stronger.
>
> Please let us know how you get on, others may face the same problem at 
> some time.
>
> Roger
>
>> Thanks again
>>
>> Murray
>>
>>
>>
>> Roger Bivand wrote:
>>> Murray Richardson <murray.richardson <at> utoronto.ca> writes:
>>>
>>>
>>>> Hello,
>>>>
>>>> I am looking for advice on a task I am trying to complete.
>>>>
>>>> I have a 4 column dataframe defining the start and end coordinates of
>>>> line edges (from a CGAL alpha shapes function to define concave hulls
>>>> from point clusters).  I would like to create polygon shapefiles from
>>>> these line edges, presumably creating lines first and then polygons.
>>>>
>>>
>>> How many polygons in each data.frame? If more than one, how are they
>>> separated? You note below that there are multiples, are they flagged by
>>> for example an NA row, or is there a jump from (endX endY) to the next
>>> (startX startY)?
>>>
>>>
>>>> e.g. columns are:
>>>>
>>>> startX  startY endX endY
>>>>
>>>> where each row represents start and end coordinates of a line segment.
>>>>
>>>
>>> Are the line segments ordered in sequence. If they are, something like:
>>>
>>> xy0 <- data.frame(sX=c(1,2,2,1,3,4,4,3), sY=c(1,1,2,2,3,3,4,4),
>>>   eX=c(2,2,1,1,4,4,3,3), eY=c(1,2,2,1,3,4,4,3))
>>>
>>> brks <- logical(nrow(xy0))
>>> for (i in 2:nrow(xy0)) brks[i] <- (xy0$sX[i] != xy0$eX[i-1]) &
>>>   (xy0$sY[i] != xy0$eY[i-1])
>>> cbrks <- cumsum(brks)+1
>>> # find the separate polygons
>>>
>>> xy1 <- split(xy0, cbrks)
>>>
>>> library(sp)
>>> Plist <- vector(mode="list", length=length(xy1))
>>> for (i in seq(along=Plist)) {
>>>   crds <- rbind(cbind(xy1[[i]]$sX, xy1[[i]]$sY),
>>>     c(xy1[[i]]$eX[nrow(xy1[[i]])], xy1[[i]]$eY[nrow(xy1[[i]])]))
>>> #   close the polygon
>>>   Plist[[i]] <- Polygons(list(Polygon(crds)), ID=names(xy1[i]))
>>> }
>>> # make a list of Polygons objects with IDs
>>>
>>> Plist1 <- SpatialPolygons(Plist)
>>> plot(Plist1, axes=TRUE)
>>> SPDF <- SpatialPolygonsDataFrame(Plist1, data=data.frame(i=names(xy1),
>>>   row.names=names(xy1)))
>>> # convert this to a SpatialPolygonsDataFrame object for export
>>>
>>> library(maptools)
>>> writePolyShape(SPDF, "chulls")
>>>
>>> If the segments are not in "join up the dots" sequence, it will be
>>> (even) more messy. Please also consider whether the R-sig-geo list
>>> might not be more relevant.
>>>
>>> Roger
>>>
>>>
>>>> I am new to R spatial packages and I am not sure which one is best
>>>> suited to this task. There seems to be a lot of options, which is 
>>>> great
>>>> but hard to know which one to start with.  Any suggestions on the best
>>>> way to proceed with this?  One further challenge is that the list
>>>> includes line segments defining multiple polygons.
>>>>
>>>> Thanks for any advice.
>>>>
>>>> Murray
>>>>
>>>> ______________________________________________
>>>> R-help <at> r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide 
>>>> http://www.R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>>>
>>>>
>>>>
>>>
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide 
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>
>> _______________________________________________
>> 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