3

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.

gisnside
  • 7,818
  • 2
  • 29
  • 73
pateto777
  • 400
  • 1
  • 4
  • 12
  • 1
    Could you tell us if you plan to use a specific software ? Python can come inside QGIS (pyqgis) or inside arcgis (arcpy) with specific algorithms that could help you out. See for example : http://gis.stackexchange.com/questions/8144/get-all-vertices-of-a-polygon-using-ogr-and-python for extraction using ogr, for example. There's lots of ways to do it. – gisnside Aug 25 '16 at 16:26
  • I'm planning to use ogr in python and R. I find very useful your comment for python, now I would like to know how can I extract the same vertices using R, I can start with: shape <- readOGR(dsn = 'input.shp', layer ='input'). How can I get the coordinates? I see them when I type: shape – pateto777 Aug 25 '16 at 16:33
  • I'm not a great R-mapper, but maybe i can help you find out :) Google tips me here : https://stat.ethz.ch/pipermail/r-sig-geo/2010-June/008520.html – gisnside Aug 25 '16 at 16:37
  • For OGR, you might be interested to read this : https://pcjericks.github.io/py-gdalogr-cookbook/geometry.html or https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html – gisnside Aug 25 '16 at 16:44

4 Answers4

3

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

Community
  • 1
pateto777
  • 400
  • 1
  • 4
  • 12
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
 ....
Community
  • 1
gene
  • 54,868
  • 3
  • 110
  • 187
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")
MacroZED
  • 2,291
  • 1
  • 12
  • 22
  • 2
    Arcpy Basic license alternative is FeatureClassToNumPyArray with explode_to_points (http://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-data-access/featureclasstonumpyarray.htm) and back again (http://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-data-access/numpyarraytofeatureclass.htm) – phloem Aug 25 '16 at 18:04
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
rcs
  • 3,894
  • 1
  • 27
  • 30