2

I need to extract the vertex list from a line inside a shapefile. I want to know which would be a better option between python and R.

4
3

1) Using Python only without a GIS with modules as Fiona (many questions/answers in GIS SE)

import fiona
lines = fiona.open("lines.shp")
# first line
first = lines.next()
print first
{'geometry': {'type': 'LineString', 'coordinates': [(270818.44773413346, 4458997.238866236), (270833.274660459, 4458983.162670357), (270830.8347865067, 4458975.280000665), (270822.3890689795, 4458967.209648361), (270823.32748203806, 4458959.702343893), (270822.7644342029, 4458958.013200387)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'ID', 1), (u'LINE_NUM_A', 11.1), (u'LINE_ST_A', u'string a')])}
print first['geometry']
{'type': 'LineString', 'coordinates': [(270818.44773413346, 4458997.238866236), (270833.274660459, 4458983.162670357), (270830.8347865067, 4458975.280000665), (270822.3890689795, 4458967.209648361), (270823.32748203806, 4458959.702343893), (270822.7644342029, 4458958.013200387)]}
print first['geometry']['coordinates']
[(270818.44773413346, 4458997.238866236), (270833.274660459, 4458983.162670357), (270830.8347865067, 4458975.280000665), (270822.3890689795, 4458967.209648361), (270823.32748203806, 4458959.702343893), (270822.7644342029, 4458958.013200387)]
vertices = first['geometry']['coordinates']

And

with fiona.open("lines.shp") as lines:
   for line in lines:
        print"vertices: ",line['geometry']['coordinates']
vertices:  [(270818.44773413346, 4458997.238866236), (270833.274660459, 4458983.162670357), (270830.8347865067, 4458975.280000665), (270822.3890689795, 4458967.209648361), (270823.32748203806, 4458959.702343893), (270822.7644342029, 4458958.013200387)]
vertices:  [(270857.2980347587, 4458989.731561768), (270866.68216534454, 4458974.904635442), (270859.92559132277, 4458964.206726574), (270852.04292163067, 4458963.080630904), (270839.8435518691, 4458948.253704579), (270832.33624740044, 4458942.99859145), (270822.7644342029, 4458938.494208769), (270822.2013863678, 4458938.681891381)]
vertices:  [(270777.90829000267, 4458973.778539772), (270778.0959726144, 4458961.954535234), (270783.7264509659, 4458943.937004508), (270789.91997715255, 4458935.866652205), (270795.3627728923, 4458931.362269524), (270802.1193469141, 4458927.233252066), (270808.50055571244, 4458923.854965055), (270815.25712973427, 4458921.790456327), (270829.5210082247, 4458918.224486704), (270837.77904314024, 4458928.359347736), (270837.96672575193, 4458933.426778253), (270834.9638039645, 4458949.755165472), (270840.0312344808, 4458962.70526568), (270840.40659970423, 4458965.895870079)]
.....

You can then use Shapely if you want (also many questions/answers in GIS SE)

from shapely.geometry import Point
for pt in vertices:
    print Point(pt)

POINT (270818.4477341335 4458997.238866236)
POINT (270833.274660459 4458983.162670357)
POINT (270830.8347865067 4458975.280000665)
POINT (270822.3890689795 4458967.209648361)
POINT (270823.3274820381 4458959.702343893)
POINT (270822.7644342029 4458958.013200387)

2) Whith R, look at Converting SpatialLinesDataFrame to SpatialPointsDataFrame in R? for example (there are others solutions)

library(rgdal)
layer = readOGR("/Users/Shared",layer="lines")
vertices  = as(layer, "SpatialPointsDataFrame")
vertices
        coordinates       ID    LINE_NUM_ A LINE_ST_A Lines.NR Lines.ID Line.NR
 0    (270818.4, 4458997)  1      11.10  string a        1        0       1
 0.1  (270833.3, 4458983)  1      11.10  string a        1        0       1
 0.2  (270830.8, 4458975)  1      11.10  string a        1        0       1
 0.3  (270822.4, 4458967)  1      11.10  string a        1        0       1
 0.4  (270823.3, 4458960)  1      11.10  string a        1        0       1
 0.5  (270822.8, 4458958)  1      11.10  string a        1        0       1
 1    (270857.3, 4458990)  2      22.20  string b        2        1       1
 1.1  (270866.7, 4458975)  2      22.20  string b        2        1       1
 ....
2

After trying OGR in python and R. I found very useful R, I was able to get the coordinates for a feature in 3 lines:

require(rgdal)
shape <- readOGR(dsn = 'input.shp', layer ='input')
coordinates(shape)[[1]][[1]]

Output:

          [,1]     [,2]
[1,] -74.03221 4.731243
[2,] -74.03235 4.725448
[3,] -74.03600 4.727154

For python: Get all vertices of a polygon using OGR and Python

0
0

Using python in Arc, you can use "Feature Vertices to Point" tool. There are multiple different options when choosing what vertices you need. Please see the following link for more information on the tool:

http://desktop.arcgis.com/en/arcmap/latest/tools/data-management-toolbox/feature-vertices-to-points.htm

Using arcpy, you could use:

import arcpy
arcpy.env.workspace = r"path/to/workspace"

arcpy.FeatureVerticesToPoints_management(polyline dataset, "name of point shapefile", "ALL")
1
0

In R you can coerce a SpatialLines object to SpatialPoints using as():

library("sp")
g1 <- rgeos::readWKT("LINESTRING(3 4,10 50,20 25)")
g1
# An object of class "SpatialLines"
# Slot "lines":
# [[1]]
# An object of class "Lines"
# Slot "Lines":
# [[1]]
# An object of class "Line"
# Slot "coords":
#       x  y
# [1,]  3  4
# [2,] 10 50
# [3,] 20 25
# 
# ...
# 

as(g1, "SpatialPoints")
# SpatialPoints:
#       x  y
# [1,]  3  4
# [2,] 10 50
# [3,] 20 25
# Coordinate Reference System (CRS) arguments: NA 

Not the answer you're looking for? Browse other questions tagged or ask your own question.