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

i have set of data represents set of nodes, each node associated value (represented color in image). want achieve selectively changing values.

the mesh represents porous system (say rock example) model. pressure in system specified @ nodes. input contains pressure attributed each node want able re-assign initial conditions pressure @ specific nodes (the ones located inside polygon). weights of nodes pressure @ node.

i want define polygon, , attribute value each vertex (think of weight), , using weight of vertex , distance vertices each node inside polygon correct value node.

this output like:

i working on algorithm takes in set of values of form [x,y,z] , [value,value,value,value]. both have same number of rows. e.i, rows in first input location of node, , rows in second values associated node.

i made algorithm takes in set of points forming polygon , set of weights corresponding each vertex of polygon.

i scan merged inputs , replace value of node found inside polygon. value defined algorithm in post. if node not inside polygon it's value preserved.

i write new_values file.

my concern not able handle large inputs of few million nodes without memoryerror. @ moment handling 9239 rows of inputs, takes me 9 seconds.

this code:

# pip problem import shapely.geometry shapely # plot import matplotlib.pyplot plt # handling data import csv import itertools # 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]) 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]) p in points)     return poly_shapely,points_shapely       def isbetween(a, b, c): #is c between , 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):     # edges     edges = []     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:         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 shapely poly , 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 = [],[],[]     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 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 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 }  data in zip(polys.itervalues(),points.itervalues()):     plot(data[0],data[1], true)  #================================================================================ #                           weighting function #================================================================================  def add_weights(poly, weights):     poly.weights = [float(w) w in weights]+[weights[0]] #need add first weight                                                           # @ end account                                                           # first point being added close 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) n,p in enumerate(poly.points_shapely) if poly.weights[n] != 'nan'])  def get_weighted_dist(poly, point):     return sum([1/distance(point,p) 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 data inside polygone #==============================================================================  '''function definitions'''  def data_extraction(filename,start_line,node_num,span_start,span_end):     open(filename, "r") myfile:         file_= csv.reader(myfile, delimiter=' ')  #extracts data .txt lines         return (x x in [filter(lambda a: != '', row[span_start:span_end]) \         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 data points     points_init = ([float(pair[0][0]),float(pair[0][2])]for pair in data)     #convert shapely poly , points     poly,points = convert_to_shapely_points_and_poly(poly_init,points_init)     add_weights(poly,weights)     #fliter out points in polygon     new_pair = []     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 x in new_pair)   def make_file(filename,data):     open(filename, "w") f:         f.writelines('\t'.join(i) + '\n' in data)         f.close()  def run(directory_path,poly_list,weight_list):     ''' directory , 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 code 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 sample data play (you need save values.txt in 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 sample of node data (will need saved 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 

i'm trying answer, not sure if suffice:

you can save list creation replacing part in edit_values:

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 x in new_pair) 

by

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] 

you save creation of new_pair list , associated memory.