| """ |
| ********* |
| Shapefile |
| ********* |
| |
| Generates a networkx.DiGraph from point and line shapefiles. |
| |
| "The Esri Shapefile or simply a shapefile is a popular geospatial vector |
| data format for geographic information systems software. It is developed |
| and regulated by Esri as a (mostly) open specification for data |
| interoperability among Esri and other software products." |
| See http://en.wikipedia.org/wiki/Shapefile for additional information. |
| """ |
| # Copyright (C) 2004-2010 by |
| # Ben Reilly <benwreilly@gmail.com> |
| # Aric Hagberg <hagberg@lanl.gov> |
| # Dan Schult <dschult@colgate.edu> |
| # Pieter Swart <swart@lanl.gov> |
| # All rights reserved. |
| # BSD license. |
| import networkx as nx |
| __author__ = """Ben Reilly (benwreilly@gmail.com)""" |
| __all__ = ['read_shp', 'write_shp'] |
| |
| |
| def read_shp(path): |
| """Generates a networkx.DiGraph from shapefiles. Point geometries are |
| translated into nodes, lines into edges. Coordinate tuples are used as |
| keys. Attributes are preserved, line geometries are simplified into start |
| and end coordinates. Accepts a single shapefile or directory of many |
| shapefiles. |
| |
| "The Esri Shapefile or simply a shapefile is a popular geospatial vector |
| data format for geographic information systems software [1]_." |
| |
| Parameters |
| ---------- |
| path : file or string |
| File, directory, or filename to read. |
| |
| Returns |
| ------- |
| G : NetworkX graph |
| |
| Examples |
| -------- |
| >>> G=nx.read_shp('test.shp') # doctest: +SKIP |
| |
| References |
| ---------- |
| .. [1] http://en.wikipedia.org/wiki/Shapefile |
| """ |
| try: |
| from osgeo import ogr |
| except ImportError: |
| raise ImportError("read_shp requires OGR: http://www.gdal.org/") |
| |
| net = nx.DiGraph() |
| |
| def getfieldinfo(lyr, feature, flds): |
| f = feature |
| return [f.GetField(f.GetFieldIndex(x)) for x in flds] |
| |
| def addlyr(lyr, fields): |
| for findex in xrange(lyr.GetFeatureCount()): |
| f = lyr.GetFeature(findex) |
| flddata = getfieldinfo(lyr, f, fields) |
| g = f.geometry() |
| attributes = dict(zip(fields, flddata)) |
| attributes["ShpName"] = lyr.GetName() |
| if g.GetGeometryType() == 1: # point |
| net.add_node((g.GetPoint_2D(0)), attributes) |
| if g.GetGeometryType() == 2: # linestring |
| attributes["Wkb"] = g.ExportToWkb() |
| attributes["Wkt"] = g.ExportToWkt() |
| attributes["Json"] = g.ExportToJson() |
| last = g.GetPointCount() - 1 |
| net.add_edge(g.GetPoint_2D(0), g.GetPoint_2D(last), attributes) |
| |
| if isinstance(path, str): |
| shp = ogr.Open(path) |
| lyrcount = shp.GetLayerCount() # multiple layers indicate a directory |
| for lyrindex in xrange(lyrcount): |
| lyr = shp.GetLayerByIndex(lyrindex) |
| flds = [x.GetName() for x in lyr.schema] |
| addlyr(lyr, flds) |
| return net |
| |
| |
| def write_shp(G, outdir): |
| """Writes a networkx.DiGraph to two shapefiles, edges and nodes. |
| Nodes and edges are expected to have a Well Known Binary (Wkb) or |
| Well Known Text (Wkt) key in order to generate geometries. Also |
| acceptable are nodes with a numeric tuple key (x,y). |
| |
| "The Esri Shapefile or simply a shapefile is a popular geospatial vector |
| data format for geographic information systems software [1]_." |
| |
| Parameters |
| ---------- |
| outdir : directory path |
| Output directory for the two shapefiles. |
| |
| Returns |
| ------- |
| None |
| |
| Examples |
| -------- |
| nx.write_shp(digraph, '/shapefiles') # doctest +SKIP |
| |
| References |
| ---------- |
| .. [1] http://en.wikipedia.org/wiki/Shapefile |
| """ |
| try: |
| from osgeo import ogr |
| except ImportError: |
| raise ImportError("write_shp requires OGR: http://www.gdal.org/") |
| # easier to debug in python if ogr throws exceptions |
| ogr.UseExceptions() |
| |
| def netgeometry(key, data): |
| if 'Wkb' in data: |
| geom = ogr.CreateGeometryFromWkb(data['Wkb']) |
| elif 'Wkt' in data: |
| geom = ogr.CreateGeometryFromWkt(data['Wkt']) |
| elif type(key[0]).__name__ == 'tuple': # edge keys are packed tuples |
| geom = ogr.Geometry(ogr.wkbLineString) |
| _from, _to = key[0], key[1] |
| try: |
| geom.SetPoint(0, *_from) |
| geom.SetPoint(1, *_to) |
| except TypeError: |
| # assume user used tuple of int and choked ogr |
| _ffrom = [float(x) for x in _from] |
| _fto = [float(x) for x in _to] |
| geom.SetPoint(0, *_ffrom) |
| geom.SetPoint(1, *_fto) |
| else: |
| geom = ogr.Geometry(ogr.wkbPoint) |
| try: |
| geom.SetPoint(0, *key) |
| except TypeError: |
| # assume user used tuple of int and choked ogr |
| fkey = [float(x) for x in key] |
| geom.SetPoint(0, *fkey) |
| |
| return geom |
| |
| # Create_feature with new optional attributes arg (should be dict type) |
| def create_feature(geometry, lyr, attributes=None): |
| feature = ogr.Feature(lyr.GetLayerDefn()) |
| feature.SetGeometry(g) |
| if attributes != None: |
| # Loop through attributes, assigning data to each field |
| for field, data in attributes.iter(): |
| feature.SetField(field, data) |
| lyr.CreateFeature(feature) |
| feature.Destroy() |
| |
| drv = ogr.GetDriverByName("ESRI Shapefile") |
| shpdir = drv.CreateDataSource(outdir) |
| # delete pre-existing output first otherwise ogr chokes |
| try: |
| shpdir.DeleteLayer("nodes") |
| except: |
| pass |
| nodes = shpdir.CreateLayer("nodes", None, ogr.wkbPoint) |
| for n in G: |
| data = G.node[n] or {} |
| g = netgeometry(n, data) |
| create_feature(g, nodes) |
| try: |
| shpdir.DeleteLayer("edges") |
| except: |
| pass |
| edges = shpdir.CreateLayer("edges", None, ogr.wkbLineString) |
| |
| # New edge attribute write support merged into edge loop |
| fields = {} # storage for field names and their data types |
| attributes = {} # storage for attribute data (indexed by field names) |
| |
| # Conversion dict between python and ogr types |
| OGRTypes = {int: ogr.OFTInteger, str: ogr.OFTString, float: ogr.OFTReal} |
| |
| # Edge loop |
| for e in G.edges(data=True): |
| data = G.get_edge_data(*e) |
| g = netgeometry(e, data) |
| # Loop through attribute data in edges |
| for key, data in e[2].iter(): |
| # Reject spatial data not required for attribute table |
| if (key != 'Json' and key != 'Wkt' and key != 'Wkb' |
| and key != 'ShpName'): |
| # For all edges check/add field and data type to fields dict |
| if key not in fields: |
| # Field not in previous edges so add to dict |
| if type(data) in OGRTypes: |
| fields[key] = OGRTypes[type(data)] |
| else: |
| # Data type not supported, default to string (char 80) |
| fields[key] = ogr.OFTString |
| # Create the new field |
| newfield = ogr.FieldDefn(key, fields[key]) |
| edges.CreateField(newfield) |
| # Store the data from new field to dict for CreateLayer() |
| attributes[key] = data |
| else: |
| # Field already exists, add data to dict for CreateLayer() |
| attributes[key] = data |
| # Create the feature with, passing new attribute data |
| create_feature(g, edges, attributes) |
| |
| nodes, edges = None, None |
| |
| |
| # fixture for nose tests |
| def setup_module(module): |
| from nose import SkipTest |
| try: |
| import ogr |
| except: |
| raise SkipTest("OGR not available") |
