Skip to content

Commit 78bdacf

Browse files
ARSimmonsARSimmons
committed
Adding initial class project work and description
1 parent 55fa913 commit 78bdacf

File tree

10 files changed

+317
-0
lines changed

10 files changed

+317
-0
lines changed

‎Students/ASimmons/class_project/code_class_project/README.md‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
2+
<html>
3+
<head>
4+
</head>
5+
<body>
6+
# Table of Contents
7+
[Team Members](#team-members)
8+
9+
[Project Summary](#project-summary)
10+
11+
# <a name="team-members"></a>Team Members
12+
* "Arielle Simmons" <[email protected]>
13+
- Data Engineer
14+
15+
# <a name="project-summary"></a>Project Summary
16+
17+
Using the following Python packages:
18+
19+
1) Fiona (1.4+) [1]
20+
2) Shapely (1.3+) [2]
21+
22+
I propose to follow the steps outlined in M. Bostock's "How To Infer Topology" [3] to design
23+
a 'PreserveTopology' module that will identify topological elements in a LineString or
24+
MultiLineString shapefile (and possibly Polygon/Multipolygons...time allowing).
25+
26+
As described in the article, the primary work that I will undertake will include:
27+
28+
1) identify junctions (intersection points) in LineString/MultiLineString shapefiles
29+
2) split (and rotate) geometric 'arcs' at their junction points.
30+
3) delete duplicate geometric 'arcs' that have no junction points but lie directly
31+
on top of an identical geometric feature.
32+
33+
34+
Key points to note:
35+
36+
- Other work that M. Bostock notes as being necessarary
37+
(...namely: "extract - decompose shapes into lines and rings"...) has already been
38+
undertaken in a previous code I wrote here: [4].
39+
- All functions will include test code.
40+
- to run from command line: python simplify_topology.py <input file>
41+
42+
References:
43+
44+
1: https://pypi.python.org/pypi/Fiona
45+
2: https://pypi.python.org/pypi/Shapely
46+
3: http://bost.ocks.org/mike/topology/
47+
4: https://github.com/ARSimmons/Shapely_Fiona_Visvalingam_Simplify
48+
49+
*NOTE: THIS PROJECT IS STILL in progress as of 11/23/2014*
50+
51+
</body>
52+
</html>

‎Students/ASimmons/class_project/code_class_project/simplify_topology.py‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
__author__ = 'asimmons'
2+
3+
import shapely
4+
import fiona
5+
from shapely.geometry import shape, LineString, Polygon, MultiLineString, MultiPolygon
6+
import math
7+
8+
9+
# g = PreserveTopology();g.find_junctions(r'I:\It_23\simplify\test_lines_simplify_topology\sample_lines.shp')
10+
11+
class PreserveTopology:
12+
13+
quantitizationFactor = (10000,10000)
14+
15+
def set_quantitization_factor(self, quantValue):
16+
self.quantitizationFactor = (quantValue, quantValue)
17+
18+
def find_junctions(self, inFile):
19+
"""
20+
Takes an 'inFile' of an ESRI shapefile, converts it into a Shapely geometry,
21+
returns a dictionary of junction points found in all shapes.
22+
23+
- A point is considered a junction if it shares the same point
24+
with another shape AND has different neighbors.
25+
26+
- Identical features on top of one another DO NOT have different
27+
neighbors, and therefor DO NOT have any junctions.
28+
"""
29+
30+
# key = quantitized junction points, value = 1
31+
dictJunctions = {}
32+
# key = checked quantitized points, value = array of quantitized neighbors
33+
dictNeighbors = {}
34+
35+
#Open input file
36+
37+
# loop over each
38+
with fiona.open(inFile, 'r') as input:
39+
40+
# Read shapely geometries from file
41+
# Loop through all shapely objects
42+
for myGeom in input:
43+
44+
myShape = shape(myGeom['geometry'])
45+
46+
if isinstance(myShape, Polygon):
47+
myShape = self.find_junctions_polygon(myShape, dictJunctions, dictNeighbors)
48+
elif isinstance(myShape, MultiPolygon):
49+
myShape = self.find_junctions_multipolygon(myShape, dictJunctions, dictNeighbors)
50+
elif isinstance(myShape, LineString):
51+
myShape = self.find_junctions_line(myShape, dictJunctions, dictNeighbors)
52+
elif isinstance(myShape, MultiLineString):
53+
myShape = self.find_junctions_multiline(myShape, dictJunctions, dictNeighbors)
54+
else:
55+
raise ValueError('Unhandled geometry type: ' + repr(myShape.type))
56+
57+
return dictJunctions
58+
59+
def quantitize(self, point):
60+
# the default quantization factor is 10,000 (-q 1e4)
61+
# Divide by quantitiztion factor, round(int), multiply by quantitization factor
62+
x_quantitized = int(round(point[0]/self.quantitizationFactor[0]) * self.quantitizationFactor[0])
63+
y_quantitized = int(round(point[1]/self.quantitizationFactor[1]) * self.quantitizationFactor[1])
64+
65+
return (x_quantitized,y_quantitized)
66+
67+
def append_junctions(self, dictJunctions, dictNeighbors, pointsList):
68+
"""
69+
70+
Builds a global dictionary of all the junctions and neighbors found in all
71+
the shapes. It determines if a point is a junction based on if it shares the same
72+
point AND has different neighbors.
73+
74+
"""
75+
# updates dictJunctions & dictNeighbors
76+
for index, point in enumerate(pointsList):
77+
quant_point = self.quantitize(point)
78+
quant_neighbors = []
79+
# append the previous neighbor
80+
if index - 1 > 0:
81+
quant_neighbors.append(self.quantitize(pointsList[index - 1]))
82+
# append the next neighbor
83+
if index + 1 < len(pointsList):
84+
quant_neighbors.append(self.quantitize(pointsList[index + 1]))
85+
86+
# check if point is in dictNeighbors, if it is
87+
# check if the neighbors are equivalent to what
88+
# is already in there, if not equiv. append to
89+
# dictJunctions
90+
91+
if quant_point in dictNeighbors:
92+
# check if neighbors are equivalent
93+
if set(dictNeighbors[quant_point]) != set(quant_neighbors):
94+
dictJunctions[quant_point] = 1
95+
else:
96+
# Otherwise, add to neighbors
97+
dictNeighbors[quant_point] = quant_neighbors
98+
99+
100+
101+
def find_junctions_line(self, myShape, dictJunctions, dictNeighbors):
102+
103+
pointsLineList = list(myShape.coords)
104+
105+
self.append_junctions(dictJunctions, dictNeighbors, pointsLineList)
106+
107+
108+
def find_junctions_multiline(self, dictJunctions, dictNeighbors):
109+
110+
raise ValueError('simplify_topology not implemented yet for: ' + repr(myShape.type))
111+
112+
def find_junctions_polygon(self, dictJunctions, dictNeighbors):
113+
114+
raise ValueError('simplify_topology not implemented yet for: ' + repr(myShape.type))
115+
116+
def find_junctions_multipolygon(self, dictJunctions, dictNeighbors):
117+
118+
raise ValueError('simplify_topology not implemented yet for: ' + repr(myShape.type))
119+
120+
# def ring_coords(self, inFile):
121+
122+
123+
# return ringPts
124+
125+
# def lines_2_arcs(self):
126+
# break lines into arcs by junctions
127+
# return
128+
129+
# def polygons_2_arcs(self):
130+
# return
131+
# break polygon rings into arcs by junctions
132+
133+
134+
# def delete_dup_ring_arcs(self):
135+
136+
# delete ring arcs that are duplicate boundaries
137+
# return
138+
139+
# def simplify_arcs(self):
140+
# simplify lines and rings by arc segments
141+
# return

‎Students/ASimmons/class_project/code_class_project/test_lines_simplify_topology/sample_lines.dbf‎

1.17 KB
Binary file not shown.

‎Students/ASimmons/class_project/code_class_project/test_lines_simplify_topology/sample_lines.prj‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
PROJCS["WGS_1984_World_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]

‎Students/ASimmons/class_project/code_class_project/test_lines_simplify_topology/sample_lines.sbn‎

140 Bytes
Binary file not shown.

‎Students/ASimmons/class_project/code_class_project/test_lines_simplify_topology/sample_lines.sbx‎

116 Bytes
Binary file not shown.

‎Students/ASimmons/class_project/code_class_project/test_lines_simplify_topology/sample_lines.shp‎

372 Bytes
Binary file not shown.

‎Students/ASimmons/class_project/code_class_project/test_lines_simplify_topology/sample_lines.shp.xml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
<?xml version="1.0" encoding="UTF-8"?>
2+
<metadata xml:lang="en"><Esri><CreaDate>20140727</CreaDate><CreaTime>13292800</CreaTime><ArcGISFormat>1.0</ArcGISFormat><SyncOnce>TRUE</SyncOnce><DataProperties><lineage><Process ToolSource="c:\program files (x86)\arcgis\desktop10.1\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131027" Time="175108">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="c:\program files (x86)\arcgis\desktop10.1\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131027" Time="184901">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="c:\program files (x86)\arcgis\desktop10.1\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131027" Time="213937">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="c:\program files (x86)\arcgis\desktop10.1\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131027" Time="214757">CalculateField oceana_us_admin1_borders ADMIN1_EDI "Edited Admin1 Class 0 to create adjacency to Class 1, ASIMMONS 10/27/2013" VB #</Process><Process ToolSource="c:\program files (x86)\arcgis\desktop10.1\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131027" Time="221853">CalculateField oceana_us_admin1_borders ADMIN1_EDI "Moved Admin1 Class 0 to create adjacency with Class 1, ASIMMONS 10/26/2013" VB #</Process><Process ToolSource="c:\program files (x86)\arcgis\desktop10.1\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131027" Time="222245">CalculateField oceana_us_admin1_borders ADMIN1_EDI "Moved Admin1 Class 0 to create adjacency with Class 1, ASIMMONS 10/26/2013" VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131028" Time="143344">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131028" Time="150735">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131028" Time="154103">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131028" Time="161939">CalculateField oceana_us_admin1_borders ADMIN1_EDI "Moved Class 0 line to be more coincident with processed p, ASimmons 10/26/2013" VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131028" Time="162228">CalculateField oceana_us_admin1_borders ADMIN1_EDI "Edited Class 0 to be coincident with processed p, ASimmons 10/28/2013" VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131028" Time="184058">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131029" Time="135611">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131029" Time="150422">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131029" Time="151113">CalculateField oceana_us_admin1_borders ADMIN1_EDI "Edited Admin1 Class 0 to Processed P, ASimmons 10/29/2013" VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131029" Time="151122">CalculateField oceana_us_admin1_borders ADMIN1_EDI "Edited Admin1 Class 0 to Processed P, ASimmons 10/29/2013" VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131029" Time="153940">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131029" Time="154926">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131029" Time="160011">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131029" Time="161035">CalculateField oceana_us_admin1_borders class 1 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131030" Time="234139">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131030" Time="234325">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131030" Time="234700">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131030" Time="234824">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="093244">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="093319">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="093900">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="095830">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="101136">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="101613">CalculateField oceana_us_admin1_borders ADMIN1_EDI "Edited Class 0 to match Processed P, ASimmons, 10/31/2013" VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="102349">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="103007">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="103615">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="104056">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="104347">CalculateField oceana_us_admin1_borders class 0 VB #</Process><Process ToolSource="d:\esri\desktop10.2\ArcToolbox\Toolboxes\Data Management Tools.tbx\CalculateField" Date="20131031" Time="104836">CalculateField oceana_us_admin1_borders class 0 VB #</Process></lineage><itemProps><itemLocation><linkage Sync="TRUE">file://\\ASIMMONS\C$\Users\asimmons\Documents\It_13\oceana_admin_1\Oceana.gdb</linkage><protocol Sync="TRUE">Local Area Network</protocol></itemLocation></itemProps></DataProperties></Esri></metadata>

‎Students/ASimmons/class_project/code_class_project/test_lines_simplify_topology/sample_lines.shx‎

116 Bytes
Binary file not shown.

‎Students/ASimmons/class_project/code_class_project/test_simplify_topology.py‎

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
__author__ = 'asimmons'
2+
3+
from simplify_topology import *
4+
from nose.tools import *
5+
import unittest
6+
7+
8+
## Test to Identify Junctions - Shared points where two or more lines intersect
9+
class test_PreserveTopology(unittest.TestCase):
10+
11+
"""
12+
cases to cover:
13+
1) one line - no junctions
14+
2) two lines - no junctions
15+
3) two lines - with junctions
16+
4) two lines that cross, but do not share a point - no junctions.
17+
5) two lines that are identical, and on top of each other, but
18+
because the neighbors are the same - no junctions.
19+
"""
20+
21+
## A
22+
## \
23+
## \
24+
## X
25+
## \
26+
## \
27+
## D
28+
29+
def test_one_line_no_junction_append_junctions(self):
30+
g = PreserveTopology()
31+
# set the quantitization factor
32+
g.set_quantitization_factor(1)
33+
# Line AXD
34+
array = [(1,3),(1.4,2),(2,0)]
35+
# create global dictionary
36+
test_dictJunctions = {}
37+
test_dictNeighbors = {}
38+
g.append_junctions(test_dictJunctions,test_dictNeighbors,array)
39+
assert(not test_dictJunctions ==True) # test that test_dictJunctions is empty
40+
assert(not test_dictNeighbors ==True) # test that test_dictNeighbors is empty
41+
42+
## A
43+
## \
44+
## \
45+
## X
46+
##
47+
##
48+
## B-----C
49+
50+
def test_two_lines_no_junctions_append_junctions(self):
51+
52+
g = PreserveTopology()
53+
# set the quantitization factor
54+
g.set_quantitization_factor(1)
55+
# LINE BC & AX
56+
array = [(0,0),(1,0)]
57+
array2 = [(1,3),(1.4,2)]
58+
# create global dictionary
59+
test_dictJunctions = {}
60+
test_dictNeighbors = {}
61+
g.append_junctions(test_dictJunctions,test_dictNeighbors,array)
62+
g.append_junctions(test_dictJunctions,test_dictNeighbors,array2)
63+
assert(not test_dictJunctions ==True) # test that test_dictJunctions is empty
64+
assert(not test_dictNeighbors ==True) # test that test_dictNeighbors is empty
65+
66+
## A
67+
## \
68+
## \
69+
## X
70+
## \
71+
## \
72+
## B-----C-----D-----F
73+
74+
def test_two_lines_one_junction(self):
75+
g = PreserveTopology()
76+
# set the quantitization factor
77+
g.set_quantitization_factor(1)
78+
# LINE BCDF & AXD
79+
array = [(0,0),(1,0),(2,0),(3,0)]
80+
array2 = [(1,3),(1.4,2),(2,0)]
81+
# create global dictionary
82+
test_dictJunctions = {}
83+
test_dictNeighbors = {}
84+
g.append_junctions(test_dictJunctions,test_dictNeighbors,array)
85+
g.append_junctions(test_dictJunctions,test_dictNeighbors,array2)
86+
assert test_dictJunctions == {(2,0): 1}
87+
88+
##
89+
## X
90+
## /
91+
## /
92+
## B-----C-----D-----F
93+
## /
94+
## Y
95+
96+
def test_two_lines_cross_but_do_not_intersect_no_junctions(self):
97+
g = PreserveTopology()
98+
# set the quantitization factor
99+
g.set_quantitization_factor(1)
100+
# these lines do NOT have a junction because there is no shared point
101+
# LINE BCDF
102+
array = [(0,0),(1,0),(2,0),(3,0)]
103+
# LINE XY
104+
array2 = [(3,3),(-1,.5)]
105+
# create global dictionary
106+
test_dictJunctions = {}
107+
test_dictNeighbors = {}
108+
g.append_junctions(test_dictJunctions,test_dictNeighbors,array)
109+
g.append_junctions(test_dictJunctions,test_dictNeighbors,array2)
110+
assert(not test_dictJunctions ==True) # test that test_dictJunctions is empty
111+
assert(not test_dictNeighbors ==True) # test that test_dictNeighbors is empty
112+
113+
114+
def test_quantitize(self):
115+
g = PreserveTopology()
116+
result = g.quantitize((12345,12345))
117+
assert_equal(result,(10000, 10000))
118+
119+
120+
if __name__ == "__main__":
121+
unittest.main()

0 commit comments

Comments
 (0)