I have a set of data which represents a set of nodes, each node is associated with a value (represented by a color in the image). What I want to achieve is selectively changing those values.
The mesh represents a porous system (say a rock for example) model. The pressure in my system is specified at the nodes. I want to be able to assign initial conditions for the pressure at specific nodes only (the ones located inside the polygon). So the weights of my nodes are the pressure at that node.
I actually want to define a polygon, and attribute a value to each vertex (think of it as a weight), and using the weight of the vertex and the distance from the vertices to each node INSIDE the polygon to correct the value for that node.
This is what my output look like:
I am working on an algorithm which takes in a set of values of the form [x,y,z] and another with [value,value,value,value]. Both have the same number of rows. e.i, the rows in the first input are the location of the node, and the rows in the second the values associated with that node.
I made an algorithm which takes in a set of points forming a polygon and a set of weights corresponding to each vertex of that polygon.
I then scan my merged inputs and replace the value of any node found to be inside the polygon. The value is defined by the algorithm in this post. If the node is not inside the polygon it's value is preserved.
I then write all the new_values to a file.
My concern is that I will not be able to handle large inputs of a few million nodes without a MemoryError. At the moment handling a 9239 rows of inputs, takes me 9 seconds.
This is my code:
# for PIP problem
import shapely.geometry as shapely
# for plot
import matplotlib.pyplot as plt
# for handling data
import csv
import itertools
# for timing
import time
#=================================================================================
# POINT IN POLYGONE PROBLEM
#=================================================================================
class MyPoly(shapely.Polygon):
def __init__(self,points):
closed_path = list(points)+[points[0]]
super(MyPoly,self).__init__(closed_path)
self.points = closed_path
self.points_shapely = [shapely.Point(p[0],p[1]) for p in closed_path]
def convert_to_shapely_points_and_poly(poly,points):
poly_shapely = MyPoly(poly)
points_shapely = (shapely.Point(p[0],p[1]) for p in points)
return poly_shapely,points_shapely
def isBetween(a, b, c): #is c between a and b ?
crossproduct = (c.y - a.y) * (b.x - a.x) - (c.x - a.x) * (b.y - a.y)
if abs(crossproduct) > 0.01 : return False # (or != 0 if using integers)
dotproduct = (c.x - a.x) * (b.x - a.x) + (c.y - a.y)*(b.y - a.y)
if dotproduct < 0 : return False
squaredlengthba = (b.x - a.x)*(b.x - a.x) + (b.y - a.y)*(b.y - a.y)
if dotproduct > squaredlengthba: return False
return True
def get_edges(poly):
# get edges
edges = []
for i in range(len(poly.points)-1):
t = [poly.points_shapely[i],poly.points_shapely[i+1]]
edges.append(t)
return edges
def inPoly(poly,point, inclusive):
if poly.contains(point) == True:
return 1
elif inclusive:
for e in get_edges(poly):
if isBetween(e[0],e[1],point):
return 1
return 0
def plot(poly_init,points_init, inclusive = True):
#convert to shapely poly and points
poly,points = convert_to_shapely_points_and_poly(poly_init,points_init)
#plot polygon
plt.plot(*zip(*poly.points))
#plot points
xs,ys,cs = [],[],[]
for point in points:
xs.append(point.x)
ys.append(point.y)
color = inPoly(poly,point, inclusive)
cs.append(color)
print point,":", color
plt.scatter(xs,ys, c = cs , s = 20*4*2)
#setting limits
axes = plt.gca()
axes.set_xlim([min(xs)-5,max(xs)+50])
axes.set_ylim([min(ys)-5,max(ys)+10])
plt.show()
# TESTS ========================================================================
#set up poly
polys = {
1 : [[10,10],[10,50],[50,50],[50,80],[100,80],[100,10]], # test rectangulary shape
2 : [[20,10],[10,20],[30,20]], # test triangle
3 : [[0,0],[0,10],[20,0],[20,10]], # test bow-tie
4 : [[0,0],[0,10],[20,10],[20,0]], # test rect clockwise
5 : [[0,0],[20,0],[20,10],[0,10]] # test rect counter-clockwise
}
#points to check
points = {
1 : [(10,25),(50,75),(60,10),(20,20),(20,60),(40,50)], # rectangulary shape test pts
2 : [[20,10],[10,20],[30,20],[-5,0],[20,15]] , # triangle test pts
3 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,5]], # bow-tie shape test pts
4 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,2],[30,8]], # rect shape test pts
5 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,2],[30,8]] # rect shape test pts
}
for data in zip(polys.itervalues(),points.itervalues()):
plot(data[0],data[1], True)
#================================================================================
# WEIGHTING FUNCTION
#================================================================================
def add_weights(poly, weights):
poly.weights = [float(w) for w in weights]+[weights[0]] #need to add the first weight
# at the end to account for
# the first point being added to close the loop
def distance(a,b):
dist = ( (b.x - a.x)**2 + (b.y - a.y)**2 )**0.5
if dist == 0: dist = 0.000000001
return dist
def get_weighted_sum(poly, point):
return sum([poly.weights[n]/distance(point,p) for n,p in enumerate(poly.points_shapely) if poly.weights[n] != 'nan'])
def get_weighted_dist(poly, point):
return sum([1/distance(point,p) for n,p in enumerate(poly.points_shapely) if poly.weights[n] != 'nan'])
def get_point_weighted_value(poly, point):
return get_weighted_sum(poly,point)/get_weighted_dist(poly,point)
#==============================================================================
# GETTING THE DATA inside the Polygone
#==============================================================================
'''Function Definitions'''
def data_extraction(filename,start_line,node_num,span_start,span_end):
with open(filename, "r") as myfile:
file_= csv.reader(myfile, delimiter=' ') #extracts data from .txt as lines
return (x for x in [filter(lambda a: a != '', row[span_start:span_end]) \
for row in itertools.islice(file_, start_line, node_num)])
def merge_data(msh_data,reload_data):
return (zip(msh_data,reload_data))
def edit_value(data, poly_init, weights):
#make x,y coordinates of the data into points
points_init = ([float(pair[0][0]),float(pair[0][2])]for pair in data)
#convert to shapely poly and points
poly,points = convert_to_shapely_points_and_poly(poly_init,points_init)
add_weights(poly,weights)
#fliter out points in polygon
new_pair = []
for n,point in enumerate(points):
if inPoly(poly, point, True):
value = str(get_point_weighted_value(poly, point))
new_pair.append([value,value,value,value])
else:
new_pair.append(data[n][1])
return (x for x in new_pair)
def make_file(filename,data):
with open(filename, "w") as f:
f.writelines('\t'.join(i) + '\n' for i in data)
f.close()
def run(directory_path,poly_list,weight_list):
''' Directory and File Names'''
msh_file = '\\nodes.txt'
reload_file = '\\values.txt'
new_reload_file = '\\new_values.txt'
dir_path = directory_path
msh_path = dir_path+msh_file
reload_path = dir_path+reload_file
'''Running the Code with your data'''
mesh_data = data_extraction(msh_path,0,54,1,4)
reload_data = data_extraction(reload_path,0,54,0,7)
data = merge_data(mesh_data,reload_data)
new_values = edit_value(data,poly_list,weight_list)
make_file(dir_path+new_reload_file,new_values)
t0 = time.time()
run("M:\\MyDocuments"
,([75,-800],[50,-900],[50,-1350],[90,-1000],[100,-900])
,(5.0e5,1e6,1e8,5.0e5,1.0e7))
t1= time.time()
print t1-t0
And this is sample data for you to play with (you will need to save it as a values.txt in the folder of interest):
1.067896706746556e+006 8.368595971460000e+006 1.068728658407457e+006 8.368595971460000e+006
2.844581459224940e+005 8.613334125294963e+006 2.846631849296530e+005 8.613337865004616e+006
1.068266636556349e+006 8.368595971460000e+006 1.069097067800019e+006 8.368595971460000e+006
2.844306728256134e+005 8.613334269264592e+006 2.846366503960088e+005 8.613338015263893e+006
2.646871122251647e+003 9.280390372578276e+006 2.647124079593603e+003 9.279361848151272e+006
2.645513962411728e+003 9.280388336827532e+006 2.645732877622660e+003 9.279359747270351e+006
1.067996132697676e+006 8.368595971460000e+006 1.068827019510901e+006 8.368595971460000e+006
1.068040363056876e+006 8.368595971460000e+006 1.068870797759632e+006 8.368595971460000e+006
1.068068562573701e+006 8.368595971460000e+006 1.068898735336173e+006 8.368595971460000e+006
1.068088894288788e+006 8.368595971460000e+006 1.068918905897983e+006 8.368595971460000e+006
1.068104407561180e+006 8.368595971460000e+006 1.068934323713974e+006 8.368595971460000e+006
1.068116587634527e+006 8.368595971460000e+006 1.068946455001944e+006 8.368595971460000e+006
1.068126287610437e+006 8.368595971460000e+006 1.068956140100951e+006 8.368595971460000e+006
1.068134059350058e+006 8.368595971460000e+006 1.068963921144020e+006 8.368595971460000e+006
1.068140293705664e+006 8.368595971460000e+006 1.068970181121935e+006 8.368595971460000e+006
1.068145286907994e+006 8.368595971460000e+006 1.068975209852979e+006 8.368595971460000e+006
1.068149274285654e+006 8.368595971460000e+006 1.068979237556288e+006 8.368595971460000e+006
1.068152448234754e+006 8.368595971460000e+006 1.068982452745734e+006 8.368595971460000e+006
1.068154968237062e+006 8.368595971460000e+006 1.068985012157432e+006 8.368595971460000e+006
1.068156966832556e+006 8.368595971460000e+006 1.068987046589365e+006 8.368595971460000e+006
1.068158553584154e+006 8.368595971460000e+006 1.068988664694503e+006 8.368595971460000e+006
1.068159818109515e+006 8.368595971460000e+006 1.068989955820237e+006 8.368595971460000e+006
1.068160832713088e+006 8.368595971460000e+006 1.068990992449035e+006 8.368595971460000e+006
1.068161654841110e+006 8.368595971460000e+006 1.068991832481593e+006 8.368595971460000e+006
1.068162329417389e+006 8.368595971460000e+006 1.068992521432464e+006 8.368595971460000e+006
1.068162891039026e+006 8.368595971460000e+006 1.068993094522149e+006 8.368595971460000e+006
1.068163365980015e+006 8.368595971460000e+006 1.068993578612505e+006 8.368595971460000e+006
1.068163773959426e+006 8.368595971460000e+006 1.068993993936813e+006 8.368595971460000e+006
1.068164129647837e+006 8.368595971460000e+006 1.068994355591033e+006 8.368595971460000e+006
1.068164443906155e+006 8.368595971460000e+006 1.068994674772899e+006 8.368595971460000e+006
1.068164724771358e+006 8.368595971460000e+006 1.068994959776965e+006 8.368595971460000e+006
1.068164978220522e+006 8.368595971460000e+006 1.068995216771951e+006 8.368595971460000e+006
1.068165208742598e+006 8.368595971460000e+006 1.068995450386572e+006 8.368595971460000e+006
1.068165419759155e+006 8.368595971460000e+006 1.068995664143054e+006 8.368595971460000e+006
1.068165613928702e+006 8.368595971460000e+006 1.068995860772074e+006 8.368595971460000e+006
1.068165793358832e+006 8.368595971460000e+006 1.068996042433189e+006 8.368595971460000e+006
1.068165959755557e+006 8.368595971460000e+006 1.068996210870328e+006 8.368595971460000e+006
1.068166114528858e+006 8.368595971460000e+006 1.068996367521690e+006 8.368595971460000e+006
1.068166258863683e+006 8.368595971460000e+006 1.068996513593785e+006 8.368595971460000e+006
1.068166393771001e+006 8.368595971460000e+006 1.068996650114520e+006 8.368595971460000e+006
1.068166520124043e+006 8.368595971460000e+006 1.068996777970799e+006 8.368595971460000e+006
1.068166638684427e+006 8.368595971460000e+006 1.068996897935547e+006 8.368595971460000e+006
1.068166750121474e+006 8.368595971460000e+006 1.068997010687638e+006 8.368595971460000e+006
1.068166855027363e+006 8.368595971460000e+006 1.068997116827437e+006 8.368595971460000e+006
1.068166953929060e+006 8.368595971460000e+006 1.068997216889020e+006 8.368595971460000e+006
1.068167047297063e+006 8.368595971460000e+006 1.068997311349124e+006 8.368595971460000e+006
1.068167135553131e+006 8.368595971460000e+006 1.068997400635036e+006 8.368595971460000e+006
1.068167219077452e+006 8.368595971460000e+006 1.068997485131862e+006 8.368595971460000e+006
1.068167298214444e+006 8.368595971460000e+006 1.068997565188462e+006 8.368595971460000e+006
1.068167373276848e+006 8.368595971460000e+006 1.068997641121696e+006 8.368595971460000e+006
1.068167444552389e+006 8.368595971460000e+006 1.068997713223352e+006 8.368595971460000e+006
1.068167512315698e+006 8.368595971460000e+006 1.068997781772706e+006 8.368595971460000e+006
1.068167576851623e+006 8.368595971460000e+006 1.068997847061655e+006 8.368595971460000e+006
1.068167638524622e+006 8.368595971460000e+006 1.068997909469268e+006 8.368595971460000e+006
1.068167697990453e+006 8.368595971460000e+006 1.068997969688582e+006 8.368595971460000e+006
and a sample of the node data (will need to be saved as nodes.txt):
0 0.26 0 -800.0
1 0.26 0 -1062.5
2 143.0 0 -800.0
3 143.0 0 -1062.5
4 0.26 0 -1150.0
5 143.0 0 -1150.0
6 1.17057404659 0 -800.0
7 2.10837283486 0 -800.0
8 3.07421037484 0 -800.0
9 4.06892500937 0 -800.0
10 5.0933801161 0 -800.0
11 6.14846485319 0 -800.0
12 7.23509501358 0 -800.0
13 8.35421377171 0 -800.0
14 9.50679247207 0 -800.0
15 10.6938315019 0 -800.0
16 11.9163611811 0 -800.0
17 13.1754426445 0 -800.0
18 14.472168711 0 -800.0
19 15.8076649025 0 -800.0
20 17.1830903981 0 -800.0
21 18.59963901 0 -800.0
22 20.0585402542 0 -800.0
23 21.5610604197 0 -800.0
24 23.1085036396 0 -800.0
25 24.7022130583 0 -800.0
26 26.3435719675 0 -800.0
27 28.0340049986 0 -800.0
28 29.7749793906 0 -800.0
29 31.5680062618 0 -800.0
30 33.4146418792 0 -800.0
31 35.3164890883 0 -800.0
32 37.2751986287 0 -800.0
33 39.2924705661 0 -800.0
34 41.3700558588 0 -800.0
35 43.5097577902 0 -800.0
36 45.7134335243 0 -800.0
37 47.9829957969 0 -800.0
38 50.3204145133 0 -800.0
39 52.7277184748 0 -800.0
40 55.2069971476 0 -800.0
41 57.7604024715 0 -800.0
42 60.3901507176 0 -800.0
43 63.0985244223 0 -800.0
44 65.8878743782 0 -800.0
45 68.7606216458 0 -800.0
46 71.7192596577 0 -800.0
47 74.7663564153 0 -800.0
48 77.904556725 0 -800.0
49 81.1365844387 0 -800.0
50 84.4652448319 0 -800.0
51 87.8934270922 0 -800.0
52 91.4241067682 0 -800.0
53 95.0603483625 0 -800.0
54 98.8053080549 0 -800.0
55 102.662236327 0 -800.0
