Replacing nodal values in a mesh with >1e6 inputs selectively using a polygon
<p>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.</p>

<p>The mesh represents a porous system (say a rock for example) model. The pressure in my system is specified at the nodes. My input already contains a pressure attributed to each node but I want to be able to re-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. </p>

<p>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.</p>

<p>This is what my output look like:</p>

<p><a href="http://i.stack.imgur.com/B9YBP.png" rel="nofollow"><img src="http://i.stack.imgur.com/B9YBP.png" alt="enter image description here"></a></p>

<p>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.</p>

<p>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.</p>

<p>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 <a href="http://codereview.stackexchange.com/questions/141546/algorithm-to-find-the-value-of-a-point-inside-a-polygone-using-weighted-vertices">post</a>. If the node is not inside the polygon it's value is preserved.</p>

<p>I then write all the new_values to a file.</p>

<p>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.</p>

<p>This is my code:</p>

<pre><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) &gt; 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 &lt; 0 : return False

    squaredlengthba = (b.x - a.x)*(b.x - a.x) + (b.y - a.y)*(b.y - a.y)
    if dotproduct &gt; 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
</code></pre>

<p>And this is sample data for you to play with (you will need to save it as a <code>values.txt</code> in the folder of interest):</p>

<pre><code>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 
</code></pre>

<p>and a sample of the node data (will need to be saved as <code>nodes.txt</code>):</p>

<pre><code>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
</code></pre>

<p>I'm trying to answer, but not sure if this will suffice:</p>

<p>You can save a list creation by replacing this part in <code>edit_values</code>:</p>

<pre><code>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)
</code></pre>

<p>by</p>

<pre><code>for n,point in enumerate(points):
    if inPoly(poly, point, True):
        value = str(get_point_weighted_value(poly, point))
        yield [value,value,value,value]
    else:
        yield data[n][1]
</code></pre>

<p>you save the creation of the <code>new_pair</code> list and the associated memory.</p>

python performance python-2.7 optimization geometry