Mercurial > repos > shellac > sam_consensus_v3
comparison env/lib/python3.9/site-packages/networkx/readwrite/nx_shp.py @ 0:4f3585e2f14b draft default tip
"planemo upload commit 60cee0fc7c0cda8592644e1aad72851dec82c959"
author | shellac |
---|---|
date | Mon, 22 Mar 2021 18:12:50 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:4f3585e2f14b |
---|---|
1 """ | |
2 ********* | |
3 Shapefile | |
4 ********* | |
5 | |
6 Generates a networkx.DiGraph from point and line shapefiles. | |
7 | |
8 "The Esri Shapefile or simply a shapefile is a popular geospatial vector | |
9 data format for geographic information systems software. It is developed | |
10 and regulated by Esri as a (mostly) open specification for data | |
11 interoperability among Esri and other software products." | |
12 See https://en.wikipedia.org/wiki/Shapefile for additional information. | |
13 """ | |
14 import networkx as nx | |
15 | |
16 __all__ = ["read_shp", "write_shp"] | |
17 | |
18 | |
19 def read_shp(path, simplify=True, geom_attrs=True, strict=True): | |
20 """Generates a networkx.DiGraph from shapefiles. Point geometries are | |
21 translated into nodes, lines into edges. Coordinate tuples are used as | |
22 keys. Attributes are preserved, line geometries are simplified into start | |
23 and end coordinates. Accepts a single shapefile or directory of many | |
24 shapefiles. | |
25 | |
26 "The Esri Shapefile or simply a shapefile is a popular geospatial vector | |
27 data format for geographic information systems software [1]_." | |
28 | |
29 Parameters | |
30 ---------- | |
31 path : file or string | |
32 File, directory, or filename to read. | |
33 | |
34 simplify: bool | |
35 If True, simplify line geometries to start and end coordinates. | |
36 If False, and line feature geometry has multiple segments, the | |
37 non-geometric attributes for that feature will be repeated for each | |
38 edge comprising that feature. | |
39 | |
40 geom_attrs: bool | |
41 If True, include the Wkb, Wkt and Json geometry attributes with | |
42 each edge. | |
43 | |
44 NOTE: if these attributes are available, write_shp will use them | |
45 to write the geometry. If nodes store the underlying coordinates for | |
46 the edge geometry as well (as they do when they are read via | |
47 this method) and they change, your geomety will be out of sync. | |
48 | |
49 strict: bool | |
50 If True, raise NetworkXError when feature geometry is missing or | |
51 GeometryType is not supported. | |
52 If False, silently ignore missing or unsupported geometry in features. | |
53 | |
54 Returns | |
55 ------- | |
56 G : NetworkX graph | |
57 | |
58 Raises | |
59 ------ | |
60 ImportError | |
61 If ogr module is not available. | |
62 | |
63 RuntimeError | |
64 If file cannot be open or read. | |
65 | |
66 NetworkXError | |
67 If strict=True and feature is missing geometry or GeometryType is | |
68 not supported. | |
69 | |
70 Examples | |
71 -------- | |
72 >>> G = nx.read_shp("test.shp") # doctest: +SKIP | |
73 | |
74 References | |
75 ---------- | |
76 .. [1] https://en.wikipedia.org/wiki/Shapefile | |
77 """ | |
78 try: | |
79 from osgeo import ogr | |
80 except ImportError as e: | |
81 raise ImportError("read_shp requires OGR: http://www.gdal.org/") from e | |
82 | |
83 if not isinstance(path, str): | |
84 return | |
85 | |
86 net = nx.DiGraph() | |
87 shp = ogr.Open(path) | |
88 if shp is None: | |
89 raise RuntimeError(f"Unable to open {path}") | |
90 for lyr in shp: | |
91 fields = [x.GetName() for x in lyr.schema] | |
92 for f in lyr: | |
93 g = f.geometry() | |
94 if g is None: | |
95 if strict: | |
96 raise nx.NetworkXError("Bad data: feature missing geometry") | |
97 else: | |
98 continue | |
99 flddata = [f.GetField(f.GetFieldIndex(x)) for x in fields] | |
100 attributes = dict(zip(fields, flddata)) | |
101 attributes["ShpName"] = lyr.GetName() | |
102 # Note: Using layer level geometry type | |
103 if g.GetGeometryType() == ogr.wkbPoint: | |
104 net.add_node((g.GetPoint_2D(0)), **attributes) | |
105 elif g.GetGeometryType() in (ogr.wkbLineString, ogr.wkbMultiLineString): | |
106 for edge in edges_from_line(g, attributes, simplify, geom_attrs): | |
107 e1, e2, attr = edge | |
108 net.add_edge(e1, e2) | |
109 net[e1][e2].update(attr) | |
110 else: | |
111 if strict: | |
112 raise nx.NetworkXError( | |
113 "GeometryType {} not supported".format(g.GetGeometryType()) | |
114 ) | |
115 | |
116 return net | |
117 | |
118 | |
119 def edges_from_line(geom, attrs, simplify=True, geom_attrs=True): | |
120 """ | |
121 Generate edges for each line in geom | |
122 Written as a helper for read_shp | |
123 | |
124 Parameters | |
125 ---------- | |
126 | |
127 geom: ogr line geometry | |
128 To be converted into an edge or edges | |
129 | |
130 attrs: dict | |
131 Attributes to be associated with all geoms | |
132 | |
133 simplify: bool | |
134 If True, simplify the line as in read_shp | |
135 | |
136 geom_attrs: bool | |
137 If True, add geom attributes to edge as in read_shp | |
138 | |
139 | |
140 Returns | |
141 ------- | |
142 edges: generator of edges | |
143 each edge is a tuple of form | |
144 (node1_coord, node2_coord, attribute_dict) | |
145 suitable for expanding into a networkx Graph add_edge call | |
146 """ | |
147 try: | |
148 from osgeo import ogr | |
149 except ImportError as e: | |
150 raise ImportError( | |
151 "edges_from_line requires OGR: " "http://www.gdal.org/" | |
152 ) from e | |
153 | |
154 if geom.GetGeometryType() == ogr.wkbLineString: | |
155 if simplify: | |
156 edge_attrs = attrs.copy() | |
157 last = geom.GetPointCount() - 1 | |
158 if geom_attrs: | |
159 edge_attrs["Wkb"] = geom.ExportToWkb() | |
160 edge_attrs["Wkt"] = geom.ExportToWkt() | |
161 edge_attrs["Json"] = geom.ExportToJson() | |
162 yield (geom.GetPoint_2D(0), geom.GetPoint_2D(last), edge_attrs) | |
163 else: | |
164 for i in range(0, geom.GetPointCount() - 1): | |
165 pt1 = geom.GetPoint_2D(i) | |
166 pt2 = geom.GetPoint_2D(i + 1) | |
167 edge_attrs = attrs.copy() | |
168 if geom_attrs: | |
169 segment = ogr.Geometry(ogr.wkbLineString) | |
170 segment.AddPoint_2D(pt1[0], pt1[1]) | |
171 segment.AddPoint_2D(pt2[0], pt2[1]) | |
172 edge_attrs["Wkb"] = segment.ExportToWkb() | |
173 edge_attrs["Wkt"] = segment.ExportToWkt() | |
174 edge_attrs["Json"] = segment.ExportToJson() | |
175 del segment | |
176 yield (pt1, pt2, edge_attrs) | |
177 | |
178 elif geom.GetGeometryType() == ogr.wkbMultiLineString: | |
179 for i in range(geom.GetGeometryCount()): | |
180 geom_i = geom.GetGeometryRef(i) | |
181 yield from edges_from_line(geom_i, attrs, simplify, geom_attrs) | |
182 | |
183 | |
184 def write_shp(G, outdir): | |
185 """Writes a networkx.DiGraph to two shapefiles, edges and nodes. | |
186 Nodes and edges are expected to have a Well Known Binary (Wkb) or | |
187 Well Known Text (Wkt) key in order to generate geometries. Also | |
188 acceptable are nodes with a numeric tuple key (x,y). | |
189 | |
190 "The Esri Shapefile or simply a shapefile is a popular geospatial vector | |
191 data format for geographic information systems software [1]_." | |
192 | |
193 Parameters | |
194 ---------- | |
195 outdir : directory path | |
196 Output directory for the two shapefiles. | |
197 | |
198 Returns | |
199 ------- | |
200 None | |
201 | |
202 Examples | |
203 -------- | |
204 nx.write_shp(digraph, '/shapefiles') # doctest +SKIP | |
205 | |
206 References | |
207 ---------- | |
208 .. [1] https://en.wikipedia.org/wiki/Shapefile | |
209 """ | |
210 try: | |
211 from osgeo import ogr | |
212 except ImportError as e: | |
213 raise ImportError("write_shp requires OGR: http://www.gdal.org/") from e | |
214 # easier to debug in python if ogr throws exceptions | |
215 ogr.UseExceptions() | |
216 | |
217 def netgeometry(key, data): | |
218 if "Wkb" in data: | |
219 geom = ogr.CreateGeometryFromWkb(data["Wkb"]) | |
220 elif "Wkt" in data: | |
221 geom = ogr.CreateGeometryFromWkt(data["Wkt"]) | |
222 elif type(key[0]).__name__ == "tuple": # edge keys are packed tuples | |
223 geom = ogr.Geometry(ogr.wkbLineString) | |
224 _from, _to = key[0], key[1] | |
225 try: | |
226 geom.SetPoint(0, *_from) | |
227 geom.SetPoint(1, *_to) | |
228 except TypeError: | |
229 # assume user used tuple of int and choked ogr | |
230 _ffrom = [float(x) for x in _from] | |
231 _fto = [float(x) for x in _to] | |
232 geom.SetPoint(0, *_ffrom) | |
233 geom.SetPoint(1, *_fto) | |
234 else: | |
235 geom = ogr.Geometry(ogr.wkbPoint) | |
236 try: | |
237 geom.SetPoint(0, *key) | |
238 except TypeError: | |
239 # assume user used tuple of int and choked ogr | |
240 fkey = [float(x) for x in key] | |
241 geom.SetPoint(0, *fkey) | |
242 | |
243 return geom | |
244 | |
245 # Create_feature with new optional attributes arg (should be dict type) | |
246 def create_feature(geometry, lyr, attributes=None): | |
247 feature = ogr.Feature(lyr.GetLayerDefn()) | |
248 feature.SetGeometry(g) | |
249 if attributes is not None: | |
250 # Loop through attributes, assigning data to each field | |
251 for field, data in attributes.items(): | |
252 feature.SetField(field, data) | |
253 lyr.CreateFeature(feature) | |
254 feature.Destroy() | |
255 | |
256 # Conversion dict between python and ogr types | |
257 OGRTypes = {int: ogr.OFTInteger, str: ogr.OFTString, float: ogr.OFTReal} | |
258 | |
259 # Check/add fields from attribute data to Shapefile layers | |
260 def add_fields_to_layer(key, value, fields, layer): | |
261 # Field not in previous edges so add to dict | |
262 if type(value) in OGRTypes: | |
263 fields[key] = OGRTypes[type(value)] | |
264 else: | |
265 # Data type not supported, default to string (char 80) | |
266 fields[key] = ogr.OFTString | |
267 # Create the new field | |
268 newfield = ogr.FieldDefn(key, fields[key]) | |
269 layer.CreateField(newfield) | |
270 | |
271 drv = ogr.GetDriverByName("ESRI Shapefile") | |
272 shpdir = drv.CreateDataSource(outdir) | |
273 # delete pre-existing output first otherwise ogr chokes | |
274 try: | |
275 shpdir.DeleteLayer("nodes") | |
276 except: | |
277 pass | |
278 nodes = shpdir.CreateLayer("nodes", None, ogr.wkbPoint) | |
279 | |
280 # Storage for node field names and their data types | |
281 node_fields = {} | |
282 | |
283 def create_attributes(data, fields, layer): | |
284 attributes = {} # storage for attribute data (indexed by field names) | |
285 for key, value in data.items(): | |
286 # Reject spatial data not required for attribute table | |
287 if key != "Json" and key != "Wkt" and key != "Wkb" and key != "ShpName": | |
288 # Check/add field and data type to fields dict | |
289 if key not in fields: | |
290 add_fields_to_layer(key, value, fields, layer) | |
291 # Store the data from new field to dict for CreateLayer() | |
292 attributes[key] = value | |
293 return attributes, layer | |
294 | |
295 for n in G: | |
296 data = G.nodes[n] | |
297 g = netgeometry(n, data) | |
298 attributes, nodes = create_attributes(data, node_fields, nodes) | |
299 create_feature(g, nodes, attributes) | |
300 | |
301 try: | |
302 shpdir.DeleteLayer("edges") | |
303 except: | |
304 pass | |
305 edges = shpdir.CreateLayer("edges", None, ogr.wkbLineString) | |
306 | |
307 # New edge attribute write support merged into edge loop | |
308 edge_fields = {} # storage for field names and their data types | |
309 | |
310 for e in G.edges(data=True): | |
311 data = G.get_edge_data(*e) | |
312 g = netgeometry(e, data) | |
313 attributes, edges = create_attributes(e[2], edge_fields, edges) | |
314 create_feature(g, edges, attributes) | |
315 | |
316 nodes, edges = None, None |