# Replacing nodal values in a mesh with >1e6 inputs selectively using a polygon

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
#================================================================================

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 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)
#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'
dir_path = directory_path
msh_path = dir_path+msh_file

'''Running the Code with your data'''
mesh_data = data_extraction(msh_path,0,54,1,4)
new_values = edit_value(data,poly_list,weight_list)

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

• Could you come up with a less generic title that actually explains what this code does? Also, you've told us a lot about how you're manipulating the data, but you've explained nothing about what you are really trying to accomplish. Commented Sep 16, 2016 at 16:51
• I have edited the question with a visual representation of what I want to achieve. I also modified the title. Commented Sep 16, 2016 at 17:02
• @Sorade: Can you say more about the motivation for this problem? What do these nodes represent? What does the polygon represent? Why do you want to change the nodes inside the polygon? What is the significance of the weighted distance to the vertices of the polygon? Commented Sep 18, 2016 at 9:10
• I've just added some information in the question. The mesh is a model of a porous system (say a rock) and the node weight is the pressure at the node. Commented Sep 18, 2016 at 9:22

I find that I get a 7× speedup simply by rewriting the inPoly function as follows:
def inPoly(poly, point, inclusive=False):

See the shapely documentation for BaseGeometry.boundary for an explanation of how this works.