source: trunk/anuga_core/source/anuga/pmesh/mesh.py @ 7872

Last change on this file since 7872 was 7872, checked in by hudson, 13 years ago

Fixed a few unit test errors.
File size: 87.4 KB
Line 
1#!/usr/bin/env python
2#
3"""General 2D triangular classes for triangular mesh generation.
4
5   Note: A .index attribute is added to objects such as vertices and
6   segments, often when reading and writing to files.  This information
7   should not be used as persistant information.  It is not the 'index' of
8   an element in a list.
9
10   
11   Copyright 2003/2004
12   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
13   Geoscience Australia
14"""
15
16import sys
17import math
18import re
19import os
20import pickle
21
22import types
23import exceptions
24
25import numpy as num
26
27
28from anuga.coordinate_transforms.geo_reference import Geo_reference, \
29     DEFAULT_ZONE
30from  anuga.load_mesh.loadASCII import NOMAXAREA, export_mesh_file, \
31     import_mesh_file
32import anuga.alpha_shape.alpha_shape
33from anuga.geospatial_data.geospatial_data import Geospatial_data, \
34     ensure_geospatial, ensure_absolute, ensure_numeric
35from anuga.mesh_engine.mesh_engine import generate_mesh
36import anuga.utilities.log as log
37
38
39try: 
40    import kinds 
41except ImportError: 
42    # Hand-built mockup of the things we need from the kinds package, since it
43    # was recently removed from the standard numeric distro.  Some users may 
44    # not have it by default. 
45    class _bunch: 
46        pass 
47         
48    class _kinds(_bunch): 
49        default_float_kind = _bunch() 
50        default_float_kind.MIN = 2.2250738585072014e-308  #smallest +ve number
51        default_float_kind.MAX = 1.7976931348623157e+308 
52     
53    kinds = _kinds()
54   
55
56# 1st and third values must be the same
57# FIXME: maybe make this a switch that the user can change? - DSG
58initialconversions = ['', 'exterior', '']
59SEG_COLOUR = 'blue'
60
61
62class MeshObject:
63    """
64    An abstract superclass for the basic mesh objects, eg vertex, segment,
65    triangle.
66    """
67    def __init__(self):
68        pass
69   
70class Point(MeshObject): 
71    """
72    Define a point in a 2D space.
73    """
74    def __init__(self,X,Y):
75        __slots__ = ['x','y']
76        self.x=X
77        self.y=Y
78       
79    def DistanceToPoint(self, OtherPoint):
80        """
81        Returns the distance from this point to another
82        """
83        SumOfSquares = ((self.x - OtherPoint.x)**2) + \
84                       ((self.y - OtherPoint.y)**2)
85        return math.sqrt(SumOfSquares)
86       
87    def IsInsideCircle(self, Center, Radius):
88        """
89        Return 1 if this point is inside the circle,
90        0 otherwise
91        """
92       
93        if (self.DistanceToPoint(Center)<Radius):
94            return 1
95        else:
96            return 0
97       
98    def __repr__(self):
99        return "(%f,%f)" % (self.x,self.y) 
100
101    def cmp_xy(self, point):
102        if self.x < point.x:
103            return -1
104        elif self.x > point.x:
105            return 1
106        else:           
107            if self.y < point.y:
108                return -1
109            elif self.y > point.y:
110                return 1
111            else:
112                return 0
113       
114    def same_x_y(self, point):
115        if self.x == point.x and self.y == point.y:
116            return True
117        else:
118            return False
119       
120           
121
122class Vertex(Point):
123    """A point on the mesh.
124    Object attributes based on the Triangle program
125    """
126
127    VERTEXSQUARESIDELENGTH = 6
128
129    def __init__(self, X, Y, attributes=None):
130        __slots__ = ['x','y','attributes']
131
132        # we don't care what we get, as long as we can get float *value*
133        self.x = float(X)
134        self.y = float(Y)
135
136        self.attributes = [] 
137        if not attributes is None:
138            self.attributes = attributes
139   
140    def setAttributes(self,attributes):
141        """attributes is a list. """
142
143        self.attributes = attributes
144       
145    def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0,
146             yoffset =0, ):
147        x =  scale*(self.x + xoffset)
148        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
149        cornerOffset= self.VERTEXSQUARESIDELENGTH/2
150
151        # A hack to see the vert tags
152        # note: there will be many tags, since tags will not be removed
153        #when zooming
154        #canvas.create_text(x+ 2*cornerOffset,
155        #                   y+ 2*cornerOffset,
156        #                        text=tags)
157       
158        # This gives points info.  It is a mess though, since numbers are
159        # not deleted when zaooming in.
160        #canvas.create_text(x+ 2*cornerOffset,
161        #                   y+ 2*cornerOffset,
162        #                        text=str(x)+','+str(y))
163       
164        return canvas.create_rectangle(x-cornerOffset,
165                                       y-cornerOffset,
166                                       x+cornerOffset,
167                                       y+cornerOffset,
168                                       tags = tags,
169                                       outline=colour,
170                                       fill = 'white')
171   
172        #return tags
173     
174    def __repr__(self):
175        return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
176
177
178class Hole(Point):
179    """A region of the mesh were no triangles are generated.
180    Defined by a point in the hole enclosed by segments.
181    """
182
183    HOLECORNERLENGTH = 6
184   
185    def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0,
186             yoffset =0, ):
187        x =  scale*(self.x + xoffset)
188        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
189        cornerOffset= self.HOLECORNERLENGTH/2
190        return canvas.create_oval(x-cornerOffset,
191                                       y-cornerOffset,
192                                       x+cornerOffset,
193                                       y+cornerOffset,
194                                       tags = tags,
195                                       outline=colour,
196                                       fill = 'white')
197   
198class Region(Point):
199    """A region of the mesh.
200    Defined by a point in the region enclosed by segments. Used to tag areas.
201    """
202
203    CROSSLENGTH = 6
204    TAG = 0
205    MAXAREA = 1
206   
207    def __init__(self, X, Y, tag=None, maxArea=None):
208        """Precondition: tag is a string and maxArea is a real"""
209
210        self.x=X
211        self.y=Y   
212        self.attributes = [] # index 0 is the tag string
213                             # optional index 1 is the max triangle area
214                             # NOTE the size of this attribute is assumed
215                             # to be 1 or 2 in regionstrings2int
216        if tag is None:
217            self.attributes.append("")
218        else:
219            self.attributes.append(tag) #this is a string
220           
221        if maxArea is not None:
222            self.setMaxArea(maxArea) # maxArea is a number
223           
224    def getTag(self,):
225        return self.attributes[self.TAG]
226   
227    def setTag(self,tag):
228        self.attributes[self.TAG] = tag
229       
230    def getMaxArea(self):
231        """ Returns the Max Area of a Triangle or
232        None, if the Max Area has not been set.
233        """
234        if self.isMaxArea():
235            return self.attributes[1]
236        else:
237            return None
238   
239    def setMaxArea(self,MaxArea):
240        if MaxArea is not None:
241            if self.isMaxArea(): 
242                self.attributes[self.MAXAREA] = float(MaxArea)
243            else:
244                self.attributes.append( float(MaxArea) )
245   
246    def deleteMaxArea(self):
247        if self.isMaxArea():
248            self.attributes.pop(self.MAXAREA)
249           
250    def isMaxArea(self):
251        return len(self.attributes)> 1
252   
253    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0,
254             colour = "black"):
255        """
256        Draw a black cross, returning the objectID
257        """
258        x =  scale*(self.x + xoffset)
259        y = -1*scale*(self.y + yoffset) 
260        cornerOffset= self.CROSSLENGTH/2
261        return canvas.create_polygon(x,
262                                     y-cornerOffset,
263                                     x,
264                                     y,
265                                     x+cornerOffset,
266                                     y,
267                                     x,
268                                     y,
269                                     x,
270                                     y+cornerOffset,
271                                     x,
272                                     y,
273                                     x-cornerOffset,
274                                     y,
275                                     x,
276                                     y,
277                                     tags = tags,
278                                     outline = colour,fill = '')
279   
280    def __repr__(self):
281        if self.isMaxArea():
282            area = self.getMaxArea() 
283            return "(%f,%f,%s,%f)" % (self.x,self.y,
284                                      self.getTag(), self.getMaxArea())
285        else:
286            return "(%f,%f,%s)" % (self.x,self.y,
287                                   self.getTag())
288       
289class Segment(MeshObject):
290    """
291    Segments are edges whose presence in the triangulation is enforced.
292   
293    """
294    def __init__(self, vertex1, vertex2, tag = None ):
295        """
296        Each segment is specified by listing the vertices of its endpoints
297        The vertices are Vertex objects
298        """
299        assert(vertex1 != vertex2)
300        self.vertices = [vertex1,vertex2 ]
301       
302        if tag is None:
303            self.tag = self.__class__.default
304        else:
305            self.tag = tag #this is a string
306       
307    def __repr__(self):
308        return "[%s,%s]" % (self.vertices,self.tag)
309           
310       
311    def draw(self, canvas, tags,scale=1, xoffset=0,
312             yoffset=0, colour=SEG_COLOUR):
313        x=[]
314        y=[]
315        for end in self.vertices:
316            #end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
317            x.append(scale*(end.x + xoffset))
318            y.append(-1*scale*(end.y + yoffset)) # - since for a canvas - is up
319       
320        return canvas.create_line(x[0], y[0], x[1], y[1],
321                                  tags = tags,fill=colour)
322    def set_tag(self,tag):
323        self.tag = tag
324       
325    # Class methods
326    def set_default_tag(cls, default):
327        cls.default = default
328   
329    def get_default_tag(cls):
330        return cls.default
331   
332    set_default_tag = classmethod(set_default_tag) 
333    get_default_tag = classmethod(get_default_tag)
334
335Segment.set_default_tag("")       
336
337class Rigid_triangulation:
338    """
339    This is a triangulation that can't have triangles added or taken away.
340
341    It just represents the triangulation, not the mesh outline needed to
342    build the triangulation.
343
344      This is the guts of the data structure;
345        generated vertex list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
346        generated segment list: [(point1,point2),(p3,p4),...]
347            (Tuples of integers)
348        generated segment tag list: [tag,tag,...] list of strings
349
350        generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
351
352        generated triangle attribute list: [s1,s2,...] list of strings
353
354        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
355            tuple of triangles
356
357            Should I do the basic structure like triangle, general
358            mesh or something else? How about like triangle, since
359            that's where the info is from, and it'll fit easier into
360            this file..
361
362             Removing the loners is difficult, since all the vert's
363             after it must be removed.
364
365             This happens in set_triangulation.
366           
367    """
368
369    def __init__(self,
370                 triangles,
371                 segments,
372                 vertices,
373                 triangle_tags,
374                 triangle_neighbors,
375                 segment_tags,
376                 vertex_attributes,
377                 vertex_attribute_titles=None
378                 ):
379       
380        self.triangles = ensure_numeric(triangles) 
381        self.triangle_neighbors = ensure_numeric(triangle_neighbors)
382        self.triangle_tags = triangle_tags # list of strings
383        self.segments = ensure_numeric(segments) 
384       
385        self.segment_tags = segment_tags # list of strings
386        self.vertices = ensure_numeric(vertices)
387        # This is needed for __cmp__
388        if vertex_attributes is None:
389            self.vertex_attributes = []
390        else:
391            self.vertex_attributes = ensure_numeric(vertex_attributes)
392        if vertex_attribute_titles is None:
393            self.vertex_attribute_titles = []
394        else:
395            self.vertex_attribute_titles = vertex_attribute_titles
396       
397    def draw_triangulation(self, canvas, scale=1, xoffset=0, yoffset=0,
398             colour="green"):
399        """
400        Draw a triangle, returning the objectID
401        """
402       
403        # FIXME(DSG-DSG) This could be a data structure that is
404        # remembered and doesn't have any duplicates.
405        # but I wouldn't be able to use create_polygon.
406        # regard it as drawing a heap of segments
407       
408        for tri in self.triangles:
409            vertices = []
410            for v_index in range(3):
411                vertices.append(self.vertices[tri[v_index]])
412           
413            objectID = canvas.create_polygon(
414                scale*(vertices[1][0] + xoffset),
415                scale*-1*(vertices[1][1] + yoffset),
416                scale*(vertices[0][0] + xoffset),
417                scale*-1*(vertices[0][1] + yoffset),
418                scale*(vertices[2][0] + xoffset),
419                scale*-1*(vertices[2][1] + yoffset),
420                outline = colour,fill = ''
421                )
422           
423    def calc_mesh_area(self):
424        area = 0
425        for tri in self.triangles:
426            vertices = []
427            for v_index in range(3):
428                vertices.append(self.vertices[tri[v_index]])
429            ax = vertices[0][0]
430            ay = vertices[0][1]
431           
432            bx = vertices[1][0]
433            by = vertices[1][1]
434           
435            cx = vertices[2][0]
436            cy = vertices[2][1]
437           
438            area += abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
439        return area           
440       
441class Mesh:
442    """
443    Representation of a 2D triangular mesh.
444    User attributes describe the mesh region/segments/vertices/attributes
445
446    mesh attributes describe the mesh that is produced eg triangles and
447    vertices.
448    All point information is relative to the geo_reference passed in
449   
450   
451    """
452
453    def __repr__(self):
454        return """
455        mesh Triangles: %s 
456        mesh Attribute Titles: %s 
457        mesh Segments: %s 
458        mesh Vertices: %s 
459        user Segments: %s 
460        user Vertices: %s 
461        holes: %s 
462        regions: %s""" % (self.meshTriangles,
463                                self.attributeTitles,
464                                self.meshSegments,
465                                self.meshVertices,
466                                self.getUserSegments(),
467                                self.userVertices,
468                                self.holes,
469                                self.regions) 
470   
471    def __init__(self,
472                 userSegments=None,
473                 userVertices=None,
474                 holes=None,
475                 regions=None,
476                 geo_reference=None):
477        self.meshTriangles=[] 
478        self.attributeTitles=[] 
479        self.meshSegments=[]
480        self.meshVertices=[]
481
482        # Rigid
483        self.tri_mesh=None
484       
485       
486        self.visualise_graph = True
487
488        if userSegments is None:
489            self.userSegments=[]
490        else:
491            self.userSegments=userSegments
492        self.alphaUserSegments=[]
493           
494        if userVertices is None:
495            self.userVertices=[]
496        else:
497            self.userVertices=userVertices
498           
499        if holes is None:
500            self.holes=[]
501        else:
502            self.holes=holes
503           
504        if regions is None:
505            self.regions=[]
506        else:
507            self.regions=regions
508
509        if geo_reference is None:
510            self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0) 
511        else:
512            self.geo_reference = geo_reference
513
514        self.shape = None
515           
516    def __cmp__(self,other):
517       
518        # A dic for the initial m
519        dic = self.Mesh2triangList()
520        dic_mesh = self.Mesh2MeshList()
521        for element in dic_mesh.keys():
522            dic[element] = dic_mesh[element]
523        for element in dic.keys():
524            dic[element].sort()
525           
526        # A dic for the exported/imported m
527        dic_other = other.Mesh2triangList()
528        dic_mesh = other.Mesh2MeshList()
529        for element in dic_mesh.keys():
530            dic_other[element] = dic_mesh[element]
531        for element in dic.keys():
532            dic_other[element].sort()
533
534        #print "dsg************************8"
535        #print "dic ",dic
536        #print "*******8"
537        #print "mesh",dic_other
538        #print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
539        #print "dsg************************8"
540       
541        return (dic.__cmp__(dic_other))
542   
543    def addUserPoint(self, pointType, x,y):
544        if pointType == Vertex:
545            point = self.addUserVertex(x,y)
546        if pointType == Hole:
547            point = self._addHole(x,y)
548        if pointType == Region:
549            point = self._addRegion(x,y)
550        return point
551   
552    def addUserVertex(self, x,y):
553        v=Vertex(x, y)
554        self.userVertices.append(v)
555        return v
556
557    def _addHole(self, x,y):
558        h=Hole(x, y)
559        self.holes.append(h)
560        return h
561   
562    def add_hole(self, x,y, geo_reference=None):
563        """
564        adds a point, which represents a hole.
565
566        The point data can have it's own geo_refernece.
567        If geo_reference is None the data is asumed to be absolute
568        """
569        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
570                                                 points_geo_ref=geo_reference)
571        return self._addHole(x, y)
572
573    def _addRegion(self, x,y):
574        h=Region(x, y)
575        self.regions.append(h)
576        return h
577   
578    def add_region(self, x,y, geo_reference=None, tag=None):
579        """
580        adds a point, which represents a region.
581
582        The point data can have it's own geo_refernece.
583        If geo_reference is None the data is asumed to be absolute
584        """
585        #FIXME: have the user set the tag and resolution here,
586        # but still return the instance, just in case.
587        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
588                                                 points_geo_ref=geo_reference)
589        region =  self._addRegion(x, y)
590        if tag is not None:
591            region.setTag(tag)
592        return region
593
594    def build_grid(self,  vert_rows, vert_columns):
595        """
596        Build a grid with vert_rows number of vertex rows and
597        vert_columns number if vertex columns
598
599        Grid spacing of 1, the origin is the lower left hand corner.
600
601        FIXME(DSG-DSG) no test.
602        """
603
604        for i in range(vert_rows):
605            for j in range(vert_columns):
606                self.addUserVertex(j,i)
607        self.auto_segment()
608        self.generateMesh(mode = "Q", minAngle=20.0)
609       
610   
611    def add_vertices(self, point_data):
612        """
613        Add user vertices.
614
615        The point_data can be a list of (x,y) values, a numeric
616        array or a geospatial_data instance.
617        """
618        point_data = ensure_geospatial(point_data)
619        #print "point_data",point_data
620        # get points relative to the mesh geo_ref
621        points = point_data.get_data_points(geo_reference=self.geo_reference)
622   
623        for point in points:
624            v=Vertex(point[0], point[1])
625            self.userVertices.append(v)
626           
627    def add_hole_from_polygon(self, polygon, segment_tags=None,
628                              geo_reference=None):
629        """
630        Add a polygon with tags to the current mesh, as a region.
631        The maxArea of the region can be specified.
632
633        If a geo_reference of the polygon points is given, this is used.
634        If not;
635        The x,y info is assumed to be Easting and Northing, absolute,
636        for the meshes zone.
637
638        polygon a list of points, in meters that describe the polygon
639             (e.g. [[x1,y1],[x2,y2],...]
640        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
641
642        This returns the region instance, so if the user whats to modify
643        it they can.       
644        """
645        return self._add_area_from_polygon(polygon,
646                                           segment_tags=segment_tags,
647                                           hole=True,
648                                           geo_reference=geo_reference
649                                           )
650
651       
652    def add_region_from_polygon(self, polygon, segment_tags=None,
653                                max_triangle_area=None, geo_reference=None,
654                                region_tag=None):
655        """
656        Add a polygon with tags to the current mesh, as a region.
657        The maxArea of the region can be specified.
658
659        If a geo_reference of the polygon points is given, this is used.
660        If not;
661        The x,y info is assumed to be Easting and Northing, absolute,
662        for the meshes zone.
663
664        polygon a list of points, in meters that describe the polygon
665             (e.g. [[x1,y1],[x2,y2],...]
666        segment_tags (e.g.{'wall':[0,1,3],'ocean':[2]}) add tags to the
667        segements of the region
668        region_tags  - add a tag to all of the triangles in the region.
669
670        This returns the region instance (if a max_triangle_area is given),
671        so if the user whats to modify it they can.     
672        """
673        if max_triangle_area is None and region_tag is None:
674            create_region = False
675        else:
676            create_region = True
677           
678        region = self._add_area_from_polygon(polygon,
679                                             segment_tags=segment_tags,
680                                             geo_reference=geo_reference,
681                                             region=create_region)
682        if max_triangle_area is not None:
683            region.setMaxArea(max_triangle_area)
684        if region_tag is not None:
685            region.setTag(region_tag)
686       
687       
688        return region
689   
690       
691    def _add_area_from_polygon(self, polygon, segment_tags=None,
692                               geo_reference=None,
693                               hole=False,
694                               region=False):
695        """
696        Add a polygon with tags to the current mesh, as a region.
697        The maxArea of the region can be specified.
698
699        If a geo_reference of the polygon points is given, this is used.
700        If not;
701        The x,y info is assumed to be Easting and Northing, absolute,
702        for the meshes zone.
703
704        polygon a list of points, in meters that describe the polygon
705             (e.g. [[x1,y1],[x2,y2],...]
706        segment_tags (e.g.{'wall':[0,1,3],'ocean':[2]})add tags to the
707        segements of the region.
708
709        This returns the region instance, so if the user whats to modify
710        it they can.
711       
712        """
713        # Only import this if necessary.
714        # Trying to get pmesh working in an uncompiled environment
715        from anuga.geometry.polygon import point_in_polygon
716       
717        #get absolute values
718        if geo_reference is not None:
719            polygon = geo_reference.get_absolute(polygon)
720        # polygon is now absolute
721        #print "polygon  should be absolute",polygon
722       
723        #create points, segs and tags
724        region_dict = {}
725        region_dict['points'] = polygon
726       
727        #Create segments
728        #E.g. [[0,1], [1,2], [2,3], [3,0]]
729        #from polygon
730        #[0,1,2,3]
731        segments = []
732        N = len(polygon)
733        for i in range(N):
734            lo = i
735            hi = (lo + 1) % N
736            segments.append( [lo, hi] ) 
737        region_dict['segments'] = segments
738        region_dict['segment_tags'] = self._tag_dict2list(segment_tags, N)
739       
740   
741        self.addVertsSegs(region_dict) #this is passing absolute values
742
743        if region is True:
744            #get inner point - absolute values
745            inner_point = point_in_polygon(polygon)
746            inner = self.add_region(inner_point[0], inner_point[1],
747                                    geo_reference=None) 
748        elif hole is True:
749            #get inner point - absolute values
750            inner_point = point_in_polygon(polygon)
751            inner = self.add_hole(inner_point[0], inner_point[1],
752                                    geo_reference=None)
753        else:
754            inner = None
755           
756        return inner
757
758    def _tag_dict2list(self, tags, number_of_segs):
759        """
760        Convert a tag dictionary from this sort of format;
761        #{'wall':[0,3],'ocean':[2]}
762
763        To a list format
764        # ['wall', '', 'ocean', 'wall']
765
766        Note: '' is a default value of nothing
767        """
768        # FIXME (DSG-DSG): Using '' as a default isn't good.
769        #Try None.
770        # Due to this default this method is too connected to
771        # _add_area_from_polygon
772        segment_tags = ['']*number_of_segs
773        if tags is not None:
774            for key in tags:
775                indices = tags[key]
776                for i in indices:
777                    segment_tags[i] = key
778        return segment_tags
779       
780    def add_circle(self, center, radius, segment_count=100,
781                   center_geo_reference=None, tag = None,
782                   region=False, hole=False):
783        """
784        center is a point, in absulute or relative to center_geo_ref
785        radius is the radius of the circle
786        segment_count is the number of segments used to represent the circle
787        tag is a string name, given to the segments.
788        If region is True a region object is returned.
789        If hole is True a hole object is returned.
790           (Don't have them both True.)
791
792           
793        """
794        # convert center and radius to a polygon
795        cuts = []
796        factor = 2* math.pi/segment_count
797        for cut in range(segment_count):
798             cuts.append(cut*factor)
799
800        polygon = []
801        for cut in cuts:
802           
803            x = center[0] + radius * math.cos(cut)
804            y = center[1] + radius * math.sin(cut)
805            polygon.append([x,y])
806        # build the tags
807        tags = {tag:range(segment_count)}
808       
809        return self._add_area_from_polygon(polygon, segment_tags=tags,
810                                           region=region, hole=hole,
811                                           geo_reference=center_geo_reference)
812
813    def auto_set_geo_reference(self):
814        """
815        Automatically set the georeference, based in the minimum x and y
816        user vertex values.
817
818        Not to be used with the graphical interface
819       
820        Not implemented.
821        Don't implement now.  using the new georeferenced points class
822        will change this?
823        """
824        #to do
825        # find the lower left hand corner
826        [xmin, ymin, xmax, ymax] = self.boxsize()
827
828        # set all points to that lower left hand corner.
829        # use change_geo_reference
830        new_geo = Geo_reference(self.geo_reference.get_zone(), xmin, ymin)
831        self.change_geo_reference(new_geo)
832       
833    def change_geo_reference(self, new_geo_reference):
834        """
835        Change from one geo_reference to another.
836        Not to be used with the graphical interface
837        """
838        # FIXME
839        # change georeference of;
840        #self.userVertices = self.geo_reference.change_points_geo_ref( \
841        #self.userVertices)
842        #self.holes = self.geo_reference.change_points_geo_ref(self.holes)
843        #self.regions = self.geo_reference.change_points_geo_ref(self.regions)
844        # The above will not work.
845        # since userVertices (etc) is a list of point objects,
846        #not a list of lists.
847        # add a method to the points class to fix this up.
848     
849    def add_segment(self, v1, v2, tag):
850        """
851        Don't do this function.
852        what will the v's be objects?  How is the user suppost to get them?
853
854        Indexes?  If add vertstosegs or add_region_from_polygon is called
855        more than once then the actual index is not obvious.  Making this
856        function confusing.
857        """
858        pass
859
860
861    def add_points_and_segments(self, points,
862                                  segments=None, segment_tags = None):
863        """
864        Add an outline of the mesh.
865        Points is a list of points a standard representation of points.
866        Segments is a list of tuples of integers.  Each tuple defines the
867           start and end of the segment by it's point index.
868        segment_tags is an optional dictionary which is used to add tags to
869           the segments.  The key is the tag name, value is the list of segment
870           indexes the tag will apply to.
871           eg. {'wall':[0,3],'ocean':[2]}
872           
873        """
874        # make sure the points are absolute
875        # Since addVertsSegs will deal with georeferencing.
876        points = ensure_absolute(points)
877
878        if segments is None:
879            segments = []
880            for i in range(len(points)-1):
881                segments.append([i, i+1])
882           
883        #create points, segs and tags
884        region_dict = {}
885        region_dict['points'] = points
886        region_dict['segments'] = segments
887        region_dict['segment_tags'] = self._tag_dict2list(segment_tags,
888                                                          len(segments))
889        self.addVertsSegs(region_dict)
890       
891    def addVertsSegs(self, outlineDict):
892        """
893        Add out-line (user Mesh) attributes given a dictionary of the lists
894        points: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
895        segments: [(point1,point2),(p3,p4),...] (Tuples of integers)
896        segment_tags: [S1Tag, S2Tag, ...] (list of strings)
897
898        Assume the values are in Eastings and Northings, with no reference
899        point. eg absolute
900        """
901        if not outlineDict.has_key('segment_tags'):
902            outlineDict['segment_tags'] = []
903            for i in range(len(outlineDict['segments'])):
904                outlineDict['segment_tags'].append('')
905        #print "outlineDict['segment_tags']",outlineDict['segment_tags']
906        #print "outlineDict['points']",outlineDict['points']
907        #print "outlineDict['segments']",outlineDict['segments']
908       
909        i_offset = len(self.userVertices)
910        #print "self.userVertices",self.userVertices
911        #print "index_offset",index_offset
912        for point in outlineDict['points']:
913            v=Vertex(point[0]-self.geo_reference.xllcorner,
914                     point[1]-self.geo_reference.yllcorner)
915            self.userVertices.append(v)
916           
917        for seg,seg_tag in map(None,outlineDict['segments'],
918                       outlineDict['segment_tags']):
919            segObject = Segment(self.userVertices[int(seg[0])+i_offset],
920                           self.userVertices[int(seg[1])+i_offset] )
921            if not seg_tag == '':
922                segObject.set_tag(seg_tag)
923            self.userSegments.append(segObject)
924           
925       
926    def get_triangle_count(self):
927        return len(self.getTriangulation())
928       
929    def getUserVertices(self):
930        """
931        Note: The x,y values will be relative to the mesh geo_reference
932        This returns a list of vertex objects
933        """
934        return self.userVertices
935
936    def get_user_vertices(self, absolute=True):
937        """
938        This returns a list of (x,y) values
939        (maybe it should return a geospatical object?
940        It makes it a bit confusing though.)
941        """
942        pointlist=[]
943        for vertex in self.userVertices:
944            pointlist.append([vertex.x,vertex.y])
945        spat = Geospatial_data(pointlist, geo_reference=self.geo_reference)
946        return spat.get_data_points(absolute=absolute)
947   
948    def getUserSegments(self):
949        allSegments = self.userSegments + self.alphaUserSegments
950        #print "self.userSegments",self.userSegments
951        #print "self.alphaUserSegments",self.alphaUserSegments
952        #print "allSegments",allSegments
953        return allSegments
954   
955    def deleteUserSegments(self,seg):
956        if self.userSegments.count(seg) == 0:
957            self.alphaUserSegments.remove(seg)
958            pass
959        else:
960            self.userSegments.remove(seg)
961           
962    def clearUserSegments(self):
963        self.userSegments = []
964        self.alphaUserSegments = []
965
966       #FIXME see where this is used. return an array instead
967    def getTriangulation(self):
968        #return self.meshTriangles
969        return self.tri_mesh.triangles.tolist()
970   
971    def getMeshVertices(self):
972        #return self.meshVertices
973        return self.tri_mesh.vertices
974 
975    def getMeshVerticeAttributes(self):
976        #return self.meshVertices
977        return self.tri_mesh.vertex_attributes
978   
979    def getMeshSegments(self):
980        #return self.meshSegments
981        return self.tri_mesh.segments
982   
983    def getMeshSegmentTags(self):
984        #return self.meshSegments
985        return self.tri_mesh.segment_tags
986   
987    def getHoles(self):
988        return self.holes
989   
990    def getRegions(self):
991        return self.regions
992   
993    def isTriangulation(self):
994        if self.meshVertices == []:
995            return False 
996        else:
997            return True
998   
999    def addUserSegment(self, v1,v2):
1000        """
1001        PRECON: A segment between the two vertices is not already present.
1002        Check by calling isUserSegmentNew before calling this function.
1003       
1004        """
1005        s=Segment( v1,v2)
1006        self.userSegments.append(s)
1007        return s
1008       
1009    def generate_mesh(self,
1010                      maximum_triangle_area="",
1011                      minimum_triangle_angle=28.0,
1012                      verbose=True):
1013        if verbose is True:
1014            silent = ''
1015        else:
1016            silent = 'Q'
1017        self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle)
1018                                  +"a"+str(maximum_triangle_area)
1019                                  +"a")
1020        #The last a is so areas for regions will be used
1021       
1022    def generateMesh(self, mode = None, maxArea = None, minAngle= None,
1023                     isRegionalMaxAreas = True):
1024        """
1025        Based on the current user vaules, holes and regions
1026        generate a new mesh
1027        mode is a string that sets conditions on the mesh generations
1028        see triangle_instructions.txt for a definition of the commands
1029       
1030        PreCondition: maxArea is a double between 1e-20 and 1e30 or is a
1031        string.
1032        """
1033        #print "mode ",mode
1034        if mode == None:
1035            self.mode = ""
1036        else:
1037            self.mode = mode
1038       
1039        if self.mode.find('p') < 0:
1040            self.mode += 'p' #p - read a planar straight line graph.
1041            #there must be segments to use this switch
1042            # TODO throw an aception if there aren't seg's
1043            # it's more comlex than this.  eg holes
1044        if self.mode.find('z') < 0:
1045            self.mode += 'z' # z - Number all items starting from zero
1046                             # (rather than one)
1047        if self.mode.find('n'):
1048            self.mode += 'n' # n - output a list of neighboring triangles
1049        if self.mode.find('A') < 0:
1050            self.mode += 'A' # A - output region attribute list for triangles
1051
1052        if not self.mode.find('V') and not self.mode.find('Q'):
1053            self.mode += 'V' # V - output info about what Triangle is doing
1054       
1055        if self.mode.find('q') < 0 and minAngle is not None:
1056            #   print "**********8minAngle******** ",minAngle
1057            min_angle = 'q' + str(minAngle)
1058            self.mode += min_angle # z - Number all items starting from zero
1059                             # (rather than one)
1060        if maxArea != None:
1061            self.mode += 'a' + str(maxArea)
1062            try:
1063                self.mode += 'a' + '%20.20f' %maxArea
1064            except TypeError:
1065                self.mode += 'a' + str(maxArea)
1066            #print "self.mode", self.mode
1067        #FIXME (DSG-DSG) This isn't explained.
1068        if isRegionalMaxAreas:
1069            self.mode += 'a'
1070        #print "mesh#generateMesh# self.mode",self.mode 
1071        meshDict = self.Mesh2triangList()
1072
1073        #FIXME (DSG-DSG)  move below section into generate_mesh.py
1074        #                  & 4 functions eg segment_strings2ints
1075        # Actually, because of region_list.append((1.0,2.0,""))
1076        # don't move it, without careful thought
1077        #print "*************************!@!@ This is going to triangle   !@!@"
1078        #print meshDict
1079        #print "************************!@!@ This is going to triangle   !@!@"
1080
1081        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist']
1082        [meshDict['segmenttaglist'],
1083         segconverter] =  segment_strings2ints(meshDict['segmenttaglist'],
1084                                             initialconversions)
1085        #print "regionlist",meshDict['regionlist']
1086        [meshDict['regionlist'],
1087         regionconverter] =  region_strings2ints(meshDict['regionlist'])
1088        #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%regionlist",meshDict['regionlist']
1089        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist'
1090        #print "self.mode", self.mode
1091        generatedMesh = generate_mesh(
1092                              meshDict['pointlist'],
1093                              meshDict['segmentlist'],
1094                              meshDict['holelist'],
1095                              meshDict['regionlist'],
1096                              meshDict['pointattributelist'],
1097                              meshDict['segmenttaglist'],
1098                              self.mode,
1099                              meshDict['pointlist'])
1100        #print "%%%%%%%%%%%%%%%%%%%%%%%%%%%generated",generatedMesh
1101        generatedMesh['qaa'] = 1
1102        if generatedMesh['generatedsegmentmarkerlist'] is not None:
1103            generatedMesh['generatedsegmentmarkerlist'] = \
1104              segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
1105                                   segconverter)
1106        #print "processed gen",generatedMesh['generatedsegmentmarkerlist']
1107        #print "pmesh mesh generatedMesh['generatedtriangleattributelist']", generatedMesh['generatedtriangleattributelist']
1108        if generatedMesh['generatedtriangleattributelist'] is not None:
1109            generatedMesh['generatedtriangleattributelist'] = \
1110            region_ints2strings(generatedMesh['generatedtriangleattributelist'],
1111                                  regionconverter)
1112
1113        #print "pmesh mesh generatedMesh['generatedtriangleattributelist']", generatedMesh['generatedtriangleattributelist']
1114        #FIXME (DSG-DSG)  move above section into generate_mesh.py
1115       
1116        if generatedMesh['generatedpointattributelist'] is None or \
1117               generatedMesh['generatedpointattributelist'].shape[1] ==0:
1118            self.attributeTitles = []
1119        generatedMesh['generatedpointattributetitlelist']= \
1120                                            self.attributeTitles
1121        #print "################  FROM TRIANGLE"
1122        #print "generatedMesh",generatedMesh
1123        #print "##################"
1124        self.setTriangulation(generatedMesh)
1125   
1126    def clearTriangulation(self):
1127
1128        #Clear the current generated mesh values
1129        self.meshTriangles=[] 
1130        self.meshSegments=[]
1131        self.meshVertices=[]
1132
1133    def removeDuplicatedUserVertices(self):
1134        """Pre-condition: There are no user segments
1135        This function will keep the first duplicate
1136        """
1137        assert self.getUserSegments() == []
1138        self.userVertices, counter =  self.removeDuplicatedVertices(
1139            self.userVertices)
1140        return counter
1141   
1142    def removeDuplicatedVertices(self, Vertices):
1143        """
1144        This function will keep the first duplicate, remove all others
1145        Precondition: Each vertex has a dupindex, which is the list
1146        index.
1147
1148        Note: this removes vertices that have the same x,y values,
1149        not duplicate instances in the Vertices list.
1150        """
1151        remove = []
1152        index = 0
1153        for v in Vertices:
1154            v.dupindex = index
1155            index += 1
1156        t = list(Vertices)
1157        t.sort(Point.cmp_xy)
1158   
1159        length = len(t)
1160        behind = 0
1161        ahead  = 1
1162        counter = 0
1163        while ahead < length:
1164            b = t[behind]
1165            ah = t[ahead]
1166            if (b.y == ah.y and b.x == ah.x):
1167                remove.append(ah.dupindex) 
1168            behind += 1
1169            ahead += 1
1170
1171        # remove the duplicate vertices
1172        remove.sort()
1173        remove.reverse()
1174        for i in remove:
1175            Vertices.pop(i)
1176            pass
1177
1178        #Remove the attribute that this function added
1179        for v in Vertices:
1180            del v.dupindex
1181        return Vertices,counter
1182
1183    # FIXME (DSG-DSG) Move this to geospatial
1184    def thinoutVertices(self, delta):
1185        """Pre-condition: There are no user segments
1186        This function will keep the first duplicate
1187        """
1188        assert self.getUserSegments() == []
1189        #t = self.userVertices
1190        #self.userVertices =[]
1191        boxedVertices = {}
1192        thinnedUserVertices =[]
1193        delta = round(delta,1)
1194       
1195        for v in self.userVertices :
1196            # tag is the center of the boxes
1197            tag = (round(v.x/delta,0)*delta,round(v.y/delta,0)*delta)
1198            #this creates a dict of lists of faces, indexed by tag
1199            boxedVertices.setdefault(tag,[]).append(v)
1200
1201        for [tag,verts] in boxedVertices.items():
1202            min = delta
1203            bestVert = None
1204            tagVert = Vertex(tag[0],tag[1])
1205            for v in verts:
1206                dist = v.DistanceToPoint(tagVert)
1207                if (dist<min):
1208                    min = dist
1209                    bestVert = v
1210            thinnedUserVertices.append(bestVert)
1211        self.userVertices =thinnedUserVertices
1212       
1213           
1214    def isUserSegmentNew(self, v1,v2):
1215        identicalSegs= [x for x in self.getUserSegments() \
1216                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1217        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1218       
1219        return len(identicalSegs) == 0
1220
1221       
1222    def deleteSegsOfVertex(self, delVertex):
1223        """
1224        Delete this vertex and any segments that connect to it.
1225        """
1226        #Find segments that connect to delVertex
1227        deletedSegments = []
1228        for seg in self.getUserSegments():
1229            if (delVertex in seg.vertices):
1230                deletedSegments.append(seg)
1231        # Delete segments that connect to delVertex
1232        for seg in deletedSegments:
1233            self.deleteUserSegments(seg)
1234        return deletedSegments
1235
1236   
1237    def deleteMeshObject(self, MeshObject):
1238        """
1239        Returns a list of all objects that were removed
1240        """
1241        deletedObs = []
1242        if isinstance(MeshObject, Vertex ):
1243            deletedObs = self.deleteSegsOfVertex(MeshObject)
1244            deletedObs.append(MeshObject)
1245            self.userVertices.remove(MeshObject)
1246        elif isinstance(MeshObject, Segment):
1247            deletedObs.append(MeshObject)
1248            self.deleteUserSegments(MeshObject)
1249        elif isinstance(MeshObject, Hole):
1250            deletedObs.append(MeshObject)
1251            self.holes.remove(MeshObject)
1252        elif isinstance(MeshObject, Region):
1253            deletedObs.append(MeshObject)
1254            self.regions.remove(MeshObject)         
1255        return deletedObs
1256                                                 
1257    def Mesh2triangList(self, userVertices=None,
1258                        userSegments=None,
1259                        holes=None,
1260                        regions=None):
1261        """
1262        Convert the Mesh to a dictionary of the lists needed for the
1263        triang module
1264        points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1265        pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
1266        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1267        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole)
1268        regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
1269       
1270        Note, this adds an index attribute to the user Vertex objects.
1271
1272        Used to produce output to triangle
1273        """
1274        if userVertices is None:
1275            userVertices = self.getUserVertices()
1276        if userSegments is None:
1277            userSegments = self.getUserSegments()
1278        if holes is None:
1279            holes = self.getHoles()
1280        if regions is None:
1281            regions = self.getRegions()
1282           
1283        meshDict = {}
1284       
1285        pointlist=[]
1286        pointattributelist=[]
1287        index = 0
1288        for vertex in userVertices:
1289            vertex.index = index
1290            pointlist.append((vertex.x,vertex.y))
1291            pointattributelist.append((vertex.attributes))
1292           
1293            index += 1
1294        meshDict['pointlist'] = pointlist
1295        meshDict['pointattributelist'] = pointattributelist
1296
1297        segmentlist=[]
1298        segmenttaglist=[]
1299        for seg in userSegments:
1300            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1301            segmenttaglist.append(seg.tag)
1302        meshDict['segmentlist'] =segmentlist
1303        meshDict['segmenttaglist'] =segmenttaglist
1304       
1305        holelist=[]
1306        for hole in holes:
1307            holelist.append((hole.x,hole.y)) 
1308        meshDict['holelist'] = holelist
1309       
1310        regionlist=[]
1311        for region in regions:
1312            if (region.getMaxArea() != None): 
1313                regionlist.append((region.x,region.y,region.getTag(),
1314                               region.getMaxArea()))
1315            else:
1316                regionlist.append((region.x,region.y,region.getTag()))
1317        meshDict['regionlist'] = regionlist
1318        #print "*(*("
1319        #print meshDict
1320        #print meshDict['regionlist']
1321        #print "*(*("
1322        return meshDict
1323                                               
1324    def Mesh2MeshList(self):
1325        """
1326        Convert the Mesh to a dictionary of lists describing the
1327        triangulation variables;
1328
1329        This is used by __cmp__
1330        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1331        generated point attribute list: [(a11,a12,...),(a21,a22),...]
1332            (Tuples of doubles)
1333        generated point attribute title list:[A1Title, A2Title ...]
1334            (list of strings)
1335        generated segment list: [(point1,point2),(p3,p4),...]
1336            (Tuples of integers)
1337        generated segment tag list: [tag,tag,...] list of strings
1338
1339        generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
1340
1341        generated triangle attribute list: [s1,s2,...] list of strings
1342
1343        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
1344            tuple of triangles
1345       
1346        Used to produce .tsh file
1347        """
1348       
1349        meshDict = {}   
1350        #print "old meshDict['generatedpointattributetitlelist']",meshDict['generatedpointattributetitlelist']
1351        #print "self.tri_mesh", self.tri_mesh
1352        if self.tri_mesh is not None:
1353            #print "self.tri_mesh.triangles.tolist()", self.tri_mesh.triangles.tolist()
1354            meshDict['generatedtrianglelist'] = self.tri_mesh.triangles.tolist()
1355            meshDict['generatedtriangleattributelist'] = self.tri_mesh.triangle_tags
1356            meshDict['generatedtriangleneighborlist'] = self.tri_mesh.triangle_neighbors.tolist()
1357            meshDict['generatedpointlist'] = self.tri_mesh.vertices.tolist()
1358            if  self.tri_mesh.vertex_attributes == []:
1359                meshDict['generatedpointattributelist'] = []
1360            #meshDict['generatedpointattributelist'] = self.tri_mesh.vertex_attributes
1361            meshDict['generatedpointattributetitlelist'] = \
1362                       self.tri_mesh.vertex_attribute_titles
1363            meshDict['generatedsegmentlist'] = self.tri_mesh.segments.tolist()
1364            meshDict['generatedsegmenttaglist'] = self.tri_mesh.segment_tags
1365        else:
1366            meshDict['generatedtrianglelist'] = []
1367            meshDict['generatedtriangleattributelist'] = []
1368            meshDict['generatedtriangleneighborlist'] = []
1369            meshDict['generatedpointlist'] = []
1370            meshDict['generatedpointattributelist'] = []
1371            meshDict['generatedpointattributetitlelist'] = []
1372            meshDict['generatedsegmentlist'] = []
1373            meshDict['generatedsegmenttaglist'] = []
1374           
1375        #print "new meshDict['generatedpointattributetitlelist']",meshDict['generatedpointattributetitlelist']
1376        #print "mesh.Mesh2MeshList*)*)"
1377        #print meshDict
1378        #print "mesh.Mesh2MeshList*)*)"
1379
1380        return meshDict
1381
1382                               
1383    def Mesh2MeshDic(self):
1384        """
1385        Convert the user and generated info of a mesh to a dictionary
1386        structure
1387        """
1388        dic = self.Mesh2triangList()
1389        dic_mesh = self.Mesh2MeshList()
1390        for element in dic_mesh.keys():
1391            dic[element] = dic_mesh[element]
1392        return dic
1393   
1394    def setTriangulation(self, genDict):
1395        """
1396        Set the mesh attributes given a dictionary of the lists
1397        returned from the triang module       
1398        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1399        generated point attribute list:[(P1att1,P1attt2, ...),
1400            (P2att1,P2attt2,...),...]- not implemented
1401        generated point attribute title list:[A1Title, A2Title ...]
1402            (list of strings) - not implemented
1403        generated segment list: [(point1,point2),(p3,p4),...]
1404            (Tuples of integers)
1405        generated segment marker list: [S1Tag, S2Tag, ...] (list of strings)
1406        triangle list:  [(point1,point2, point3),(p5,p4, p1),...]
1407            (Tuples of integers)
1408        triangle neighbor list: [(triangle1,triangle2, triangle3),
1409            (t5,t4, t1),...] (Tuples of integers)
1410            -1 means there's no triangle neighbor
1411        triangle attribute list: [(T1att), (T2att), ...](list of strings)
1412
1413        """
1414        # Setting up the rigid triangulation
1415        self.tri_mesh = Rigid_triangulation(
1416            genDict['generatedtrianglelist']
1417            ,genDict['generatedsegmentlist']
1418            ,genDict['generatedpointlist']
1419            ,genDict['generatedtriangleattributelist']
1420            ,genDict['generatedtriangleneighborlist']
1421            ,genDict['generatedsegmentmarkerlist']
1422            ,genDict['generatedpointattributelist']
1423            ,genDict['generatedpointattributetitlelist']
1424            )
1425           
1426    def setMesh(self, genDict):
1427        """
1428        Set the user Mesh attributes given a dictionary of the lists
1429        point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1430        point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1431        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1432        segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1433        region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1434        region attribute list: ["","reservoir",""] list of strings
1435        region max area list:[real, None, Real,...] list of None and reals
1436       
1437        mesh is an instance of a mesh object
1438        """
1439        #Clear the current user mesh values
1440        self.clearUserSegments()
1441        self.userVertices=[]
1442        self.Holes=[]
1443        self.Regions=[]
1444
1445        #print "mesh.setMesh@#@#@#"
1446        #print genDict
1447        #print "@#@#@#"
1448       
1449        #index = 0
1450        for point in genDict['pointlist']:
1451            v=Vertex(point[0], point[1])
1452            #v.index = index
1453            #index +=1
1454            self.userVertices.append(v)
1455
1456        #index = 0
1457        for seg,tag in map(None,genDict['segmentlist'],
1458                           genDict['segmenttaglist']):
1459            segObject = Segment( self.userVertices[int(seg[0])],
1460                           self.userVertices[int(seg[1])], tag = tag )
1461            #segObject.index = index
1462            #index +=1
1463            self.userSegments.append(segObject)
1464
1465# Remove the loading of attribute info.
1466# Have attribute info added using least_squares in pyvolution
1467#         index = 0
1468#         for att in genDict['pointattributelist']:
1469#             if att == None:
1470#                 self.userVertices[index].setAttributes([])
1471#             else:
1472#                 self.userVertices[index].setAttributes(att)
1473#            index += 1
1474       
1475        #index = 0
1476        for point in genDict['holelist']:
1477            h=Hole(point[0], point[1])
1478            #h.index = index
1479            #index +=1
1480            self.holes.append(h)
1481
1482        #index = 0
1483        for reg,att,maxArea in map(None,
1484                                   genDict['regionlist'],
1485                                   genDict['regionattributelist'],
1486                                   genDict['regionmaxarealist']):
1487            Object = Region( reg[0],
1488                             reg[1],
1489                             tag = att,
1490                             maxArea = maxArea)
1491            #Object.index = index
1492            #index +=1
1493            self.regions.append(Object)
1494           
1495    def Testauto_segment(self):
1496        newsegs = []
1497        s1 = Segment(self.userVertices[0],
1498                               self.userVertices[1])
1499        s2 = Segment(self.userVertices[0],
1500                               self.userVertices[2])
1501        s3 = Segment(self.userVertices[2],
1502                               self.userVertices[1])
1503        if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1504            newsegs.append(s1)
1505        if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1506            newsegs.append(s2)
1507        if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1508            newsegs.append(s3)
1509        #DSG!!!
1510        self.userSegments.extend(newsegs)
1511        return newsegs
1512
1513   
1514    def savePickle(self, currentName):
1515        fd = open(currentName, 'w')
1516        pickle.dump(self,fd)
1517        fd.close()
1518
1519    def auto_segmentFilter(self,raw_boundary=True,
1520                    remove_holes=False,
1521                    smooth_indents=False,
1522                    expand_pinch=False):
1523        """
1524        Change the filters applied on the alpha shape boundary
1525        """
1526        if self.shape is None:
1527            return [],[],0.0
1528        return self._boundary2mesh(raw_boundary=raw_boundary,
1529                    remove_holes=remove_holes,
1530                    smooth_indents=smooth_indents,
1531                    expand_pinch=expand_pinch)
1532   
1533       
1534   
1535    def auto_segment(self, alpha = None,
1536                    raw_boundary=True,
1537                    remove_holes=False,
1538                    smooth_indents=False,
1539                    expand_pinch=False): 
1540        """
1541        Precon: There must be 3 or more vertices in the userVertices structure
1542       
1543        This returns alpha_segs_no_user_segs, segs2delete, optimum_alpha
1544        """
1545        self._createBoundary(alpha=alpha)
1546        return self._boundary2mesh(raw_boundary=raw_boundary,
1547                    remove_holes=remove_holes,
1548                    smooth_indents=smooth_indents,
1549                    expand_pinch=expand_pinch)
1550
1551    def _createBoundary(self,alpha=None):
1552        """
1553        """
1554        points=[]
1555        for vertex in self.getUserVertices():
1556            points.append((vertex.x,vertex.y))
1557        self.shape = anuga.alpha_shape.alpha_shape.Alpha_Shape(points,
1558                                                               alpha=alpha)
1559
1560
1561    def _boundary2mesh(self, raw_boundary=True,
1562                    remove_holes=False,
1563                    smooth_indents=False,
1564                    expand_pinch=False):
1565        """
1566        Precon there must be a shape object.
1567        """
1568        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1569                                 remove_holes=remove_holes,
1570                                 smooth_indents=smooth_indents,
1571                                 expand_pinch=expand_pinch)
1572        boundary_segs = self.shape.get_boundary()
1573        #print "boundary_segs",boundary_segs
1574        segs2delete = self.alphaUserSegments
1575        #FIXME(DSG-DSG) this algorithm needs comments
1576        new_segs = {}
1577        #alpha_segs = []
1578        #user_segs = []
1579        for seg in boundary_segs:
1580            v1 = self.userVertices[int(seg[0])]
1581            v2 = self.userVertices[int(seg[1])]
1582            boundary_seg = Segment(v1, v2)
1583            new_segs[(v1,v2)] = boundary_seg
1584
1585        for user_seg in self.userSegments:
1586            if new_segs.has_key((user_seg.vertices[0],
1587                                user_seg.vertices[1])):
1588                del new_segs[user_seg.vertices[0],
1589                                user_seg.vertices[1]]
1590            elif new_segs.has_key((user_seg.vertices[1],
1591                                user_seg.vertices[0])):
1592                del new_segs[user_seg.vertices[1],
1593                                user_seg.vertices[0]]
1594               
1595        optimum_alpha = self.shape.get_alpha()
1596        alpha_segs_no_user_segs  = new_segs.values()
1597        self.alphaUserSegments = alpha_segs_no_user_segs
1598        return alpha_segs_no_user_segs, segs2delete, optimum_alpha
1599   
1600   
1601    def representedAlphaUserSegment(self, v1,v2):
1602        identicalSegs= [x for x in self.alphaUserSegments \
1603                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1604        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1605
1606        if identicalSegs == []:
1607            return None
1608        else:
1609            # Only return the first one.
1610            return identicalSegs[0]
1611   
1612    def representedUserSegment(self, v1,v2):
1613        identicalSegs= [x for x in self.userSegments \
1614                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1615        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1616
1617        if identicalSegs == []:
1618            return None
1619        else:
1620            # Only return the first one.
1621            return identicalSegs[0]
1622       
1623    def joinVertices(self):
1624        """
1625        Return list of segments connecting all userVertices
1626        in the order they were given
1627       
1628        Precon: There must be 3 or more vertices in the userVertices structure
1629        """
1630
1631        newsegs = []
1632       
1633        v1 = self.userVertices[0]
1634        for v2 in self.userVertices[1:]:
1635            if self.isUserSegmentNew(v1,v2):           
1636                newseg = Segment(v1, v2)
1637                newsegs.append(newseg)
1638            v1 = v2
1639
1640        #Connect last point to the first
1641        v2 = self.userVertices[0]       
1642        if self.isUserSegmentNew(v1,v2):           
1643                newseg = Segment(v1, v2)
1644                newsegs.append(newseg)
1645           
1646
1647        #Update list of user segments       
1648        #DSG!!!
1649        self.userSegments.extend(newsegs)               
1650        return newsegs
1651   
1652    def normaliseMesh(self,scale, offset, height_scale):
1653        [xmin, ymin, xmax, ymax] = self.boxsize()
1654        [attmin0, attmax0] = self.maxMinVertAtt(0)
1655        #print "[attmin0, attmax0]" ,[attmin0, attmax0]
1656        [attmin1, attmax1] = self.maxMinVertAtt(1)
1657        #print [xmin, ymin, xmax, ymax]
1658        xrange = xmax - xmin
1659        yrange = ymax - ymin
1660        if xrange > yrange:
1661            min,max = xmin, xmax
1662        else:
1663            min,max = ymin, ymax
1664           
1665        for obj in self.getUserVertices():
1666            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1667            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1668            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1669                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1670                                    (attmax0-attmin0)*height_scale
1671            if len(obj.attributes) > 1 and attmin1 != attmax1:
1672                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1673                                    (attmax1-attmin1)*height_scale
1674           
1675        for obj in self.getMeshVertices():
1676            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1677            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1678            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1679                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1680                                    (attmax0-attmin0)*height_scale
1681            if len(obj.attributes) > 1 and attmin1 != attmax1:
1682                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1683                                    (attmax1-attmin1)*height_scale
1684               
1685        for obj in self.getHoles():
1686            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1687            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1688        for obj in self.getRegions():
1689            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1690            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1691        [xmin, ymin, xmax, ymax] = self.boxsize()
1692        #print [xmin, ymin, xmax, ymax]
1693     
1694    def boxsizeVerts(self):
1695        """
1696        Returns a list of verts denoting a box or triangle that contains
1697        verts on the xmin, ymin, xmax and ymax axis.
1698        Structure: list of verts
1699        """
1700       
1701        large = kinds.default_float_kind.MAX
1702        xmin= large
1703        xmax=-large
1704        ymin= large
1705        ymax=-large
1706        for vertex in self.userVertices:
1707            if vertex.x < xmin:
1708                xmin = vertex.x
1709                xminVert = vertex
1710            if vertex.x > xmax:
1711                xmax = vertex.x
1712                xmaxVert = vertex
1713               
1714            if vertex.y < ymin:
1715                ymin = vertex.y
1716                yminVert = vertex
1717            if vertex.y > ymax:
1718                ymax = vertex.y
1719                ymaxVert = vertex
1720        verts, count = self.removeDuplicatedVertices([xminVert,
1721                                                      xmaxVert,
1722                                                      yminVert,
1723                                                      ymaxVert])
1724         
1725        return verts
1726   
1727    def boxsize(self):
1728        """
1729        Returns a list denoting a box that contains the entire
1730        structure of vertices
1731        Structure: [xmin, ymin, xmax, ymax]
1732        """
1733     
1734        large = kinds.default_float_kind.MAX
1735        xmin= large
1736        xmax=-large
1737        ymin= large
1738        ymax=-large
1739        for vertex in self.userVertices:
1740            if vertex.x < xmin:
1741                xmin = vertex.x
1742            if vertex.x > xmax:
1743                xmax = vertex.x
1744               
1745            if vertex.y < ymin:
1746                ymin = vertex.y
1747            if vertex.y > ymax:
1748                ymax = vertex.y
1749        return [xmin, ymin, xmax, ymax]
1750 
1751    def maxMinVertAtt(self, iatt):
1752        """
1753        Returns a list denoting a box that contains the entire structure
1754        of vertices
1755        Structure: [xmin, ymin, xmax, ymax]
1756        """
1757       
1758        large = kinds.default_float_kind.MAX
1759        min= large
1760        max=-large
1761        for vertex in self.userVertices:
1762            if len(vertex.attributes) > iatt:
1763                if vertex.attributes[iatt] < min:
1764                    min = vertex.attributes[iatt]
1765                if vertex.attributes[iatt] > max:
1766                    max = vertex.attributes[iatt] 
1767        for vertex in self.meshVertices:
1768            if len(vertex.attributes) > iatt:
1769                if vertex.attributes[iatt] < min:
1770                    min = vertex.attributes[iatt]
1771                if vertex.attributes[iatt] > max:
1772                    max = vertex.attributes[iatt] 
1773        return [min, max]
1774   
1775    def scaleoffset(self, WIDTH, HEIGHT):
1776        """
1777        Returns a list denoting the scale and offset terms that need to be
1778        applied when converting  mesh co-ordinates onto grid co-ordinates
1779        Structure: [scale, xoffset, yoffset]
1780        """   
1781        OFFSET = 0.05*min([WIDTH, HEIGHT])
1782        [xmin, ymin, xmax, ymax] = self.boxsize()
1783        SCALE = min([0.9*WIDTH, 0.9*HEIGHT])/max([xmax-xmin, ymax-ymin])
1784       
1785        if SCALE*xmin < OFFSET:
1786            xoffset = abs(SCALE*xmin) + OFFSET
1787        if SCALE*xmax > WIDTH - OFFSET:
1788            xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1789        if SCALE*ymin < OFFSET:
1790            b = abs(SCALE*ymin)+OFFSET
1791        if SCALE*ymax > HEIGHT-OFFSET:
1792            b = -(SCALE*ymax - HEIGHT + OFFSET)
1793        yoffset = HEIGHT - b
1794        return [SCALE, xoffset, yoffset]
1795
1796         
1797    def exportASCIIobj(self,ofile):
1798        """
1799        export a file, ofile, with the format
1800         lines:  v <x> <y> <first attribute>
1801        f <vertex #>  <vertex #> <vertex #> (of the triangles)
1802        """
1803        fd = open(ofile,'w')
1804        self.writeASCIIobj(fd)   
1805        fd.close()
1806
1807
1808    def writeASCIIobj(self,fd):
1809        fd.write(" # Triangulation as an obj file\n")
1810        numVert = str(len(self.meshVertices))
1811       
1812        index1 = 1 
1813        for vert in self.meshVertices:
1814            vert.index1 = index1
1815            index1 += 1
1816           
1817            fd.write("v "
1818                     + str(vert.x) + " "
1819                     + str(vert.y) + " "
1820                     + str(vert.attributes[0]) + "\n")
1821           
1822        for tri in self.meshTriangles:
1823            fd.write("f "
1824                     + str(tri.vertices[0].index1) + " " 
1825                     + str(tri.vertices[1].index1) + " " 
1826                     + str(tri.vertices[2].index1) + "\n")
1827           
1828    def exportASCIIsegmentoutlinefile(self,ofile):
1829        """Write the boundary user mesh info, eg
1830         vertices that are connected to segments,
1831         segments
1832        """
1833       
1834        verts = {}
1835        for seg in self.getUserSegments():
1836            verts[seg.vertices[0]] = seg.vertices[0]
1837            verts[seg.vertices[1]] = seg.vertices[1]
1838        meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values())
1839        export_mesh_file(ofile,meshDict)
1840       
1841        # exportASCIImeshfile   - this function is used
1842    def export_mesh_file(self,ofile):
1843        """
1844        export a file, ofile, with the format
1845        """
1846       
1847        dict = self.Mesh2IODict()
1848        export_mesh_file(ofile,dict)
1849
1850    # FIXME(DSG-DSG):Break this into two functions.
1851    #One for the outline points.
1852    #One for the mesh points.
1853    # Note: this function is not in the gui
1854    def exportPointsFile(self,ofile):
1855        """
1856        export a points file, ofile.
1857       
1858        """
1859       
1860        mesh_dict = self.Mesh2IODict()
1861        #point_dict = {}
1862        #point_dict['attributelist'] = {} #this will need to be expanded..
1863                                         # if attributes are brought back in.
1864        #point_dict['geo_reference'] = self.geo_reference
1865        if mesh_dict['vertices'] == []:
1866            #point_dict['pointlist'] = mesh_dict['points']
1867            geo = Geospatial_data(mesh_dict['points'],
1868                                  geo_reference=self.geo_reference)
1869        else:
1870            #point_dict['pointlist'] = mesh_dict['vertices']
1871            geo = Geospatial_data(mesh_dict['vertices'],
1872                                  geo_reference=self.geo_reference)
1873
1874        geo.export_points_file(ofile, absolute=True)
1875       
1876
1877
1878    def import_ungenerate_file(self,ofile, tag=None, region_tag=None):
1879        """
1880        Imports an ungenerate file, from arcGIS into mesh.
1881        This file describes many polygons.
1882
1883        ofile is the name of the ungenerated file.
1884        Tag is a string name to be taggged on each segment.
1885       
1886        region_tag is the tag applied to the polygon regions.
1887        if it is a string the one value will be assigned to all regions
1888        if it is a list the first value in the list will be applied to the first polygon etc.
1889        WARNING: size of list and number of polygons isn't checked
1890
1891        WARNING values are assumed to be absolute.
1892        geo-refs are not taken into account..
1893        """
1894   
1895        dict = importUngenerateFile(ofile)
1896        default_tag = Segment.get_default_tag()
1897        if tag is not None:
1898            Segment.set_default_tag(str(tag))
1899       
1900        if region_tag is None:
1901            self.addVertsSegs(dict)
1902        else:           
1903            if not isinstance(region_tag, list):
1904                region_tag = [region_tag]*len(dict['polygons'])
1905            for a_tag,polygon in map(None, region_tag, dict['polygons']): 
1906                segment_tags = {tag:range(len(polygon))}
1907                self.add_region_from_polygon(polygon,segment_tags=segment_tags,
1908                                region_tag=a_tag)
1909           
1910       
1911        Segment.set_default_tag(default_tag)
1912
1913        # change the tag back to it's default
1914   
1915       
1916########### IO CONVERTERS ##################
1917        """
1918        The dict fromat for IO with .tsh files is;
1919        (the triangulation)
1920        vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
1921        vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1922        vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
1923        segments: [[v1,v2],[v3,v4],...] (lists of integers)
1924        segment_tags : [tag,tag,...] list of strings
1925        triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
1926        triangle tags: [s1,s2,...] A list of strings
1927        triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
1928       
1929        (the outline)   
1930        points: [[x1,y1],[x2,y2],...] (lists of doubles)
1931        point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1932        outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
1933        outline_segment_tags : [tag1,tag2,...] list of strings
1934        holes : [[x1,y1],...](List of doubles, one inside each hole region)
1935        regions : [ [x1,y1],...] (List of 4 doubles)
1936        region_tags : [tag1,tag2,...] (list of strings)
1937        region_max_areas: [ma1,ma2,...] (A list of doubles)
1938        {Convension: A -ve max area means no max area}
1939       
1940        """
1941     
1942
1943                               
1944    def Mesh2IODict(self):
1945        """
1946        Convert the triangulation and outline info of a mesh to a dictionary
1947        structure
1948        """
1949        dict = self.Mesh2IOTriangulationDict()
1950        dict_mesh = self.Mesh2IOOutlineDict()
1951        for element in dict_mesh.keys():
1952            dict[element] = dict_mesh[element]
1953
1954        # add the geo reference
1955        dict['geo_reference'] = self.geo_reference
1956        return dict
1957   
1958    def Mesh2IOTriangulationDict(self):
1959        """
1960        Convert the Mesh to a dictionary of lists describing the
1961        triangulation variables;
1962       
1963        Used to produce .tsh file
1964        """
1965        meshDict = {}
1966        if self.tri_mesh is not None:
1967            meshDict['triangles'] = self.tri_mesh.triangles
1968            meshDict['triangle_tags'] = self.tri_mesh.triangle_tags
1969            #print "mesh meshDict['triangle_tags']", meshDict['triangle_tags']
1970            meshDict['triangle_neighbors'] = self.tri_mesh.triangle_neighbors
1971            meshDict['vertices'] = self.tri_mesh.vertices
1972            meshDict['vertex_attributes'] = self.tri_mesh.vertex_attributes
1973            meshDict['vertex_attribute_titles'] = \
1974                               self.tri_mesh.vertex_attribute_titles
1975            meshDict['segments'] = self.tri_mesh.segments
1976            meshDict['segment_tags'] = self.tri_mesh.segment_tags
1977        else:
1978            meshDict['triangles'] = []
1979            meshDict['triangle_tags'] = []
1980            meshDict['triangle_neighbors'] = []
1981            meshDict['vertices'] = []
1982            meshDict['vertex_attributes'] = []
1983            meshDict['vertex_attribute_titles'] = []
1984            meshDict['segments'] = []
1985            meshDict['segment_tags'] = []
1986        #print "mesh.Mesh2IOTriangulationDict*)*)"
1987        #print meshDict
1988        #print "mesh.Mesh2IOTriangulationDict*)*)"
1989
1990        return meshDict
1991
1992                                                     
1993    def Mesh2IOOutlineDict(self, userVertices=None,
1994                        userSegments=None,
1995                        holes=None,
1996                        regions=None):
1997        """
1998        Convert the mesh outline to a dictionary of the lists needed for the
1999        triang module;
2000       
2001        Note, this adds an index attribute to the user Vertex objects.
2002
2003        Used to produce .tsh file and output to triangle
2004        """
2005       
2006        if userVertices is None:
2007            userVertices = self.getUserVertices()
2008        if userSegments is None:
2009            userSegments = self.getUserSegments()
2010        if holes is None:
2011            holes = self.getHoles()
2012        if regions is None:
2013            regions = self.getRegions()
2014           
2015        meshDict = {}
2016        #print "userVertices",userVertices
2017        #print "userSegments",userSegments
2018        pointlist=[]
2019        pointattributelist=[]
2020        index = 0
2021        for vertex in userVertices:
2022            vertex.index = index
2023            pointlist.append([vertex.x,vertex.y])
2024            pointattributelist.append(vertex.attributes)
2025           
2026            index += 1
2027        meshDict['points'] = pointlist
2028        meshDict['point_attributes'] = pointattributelist
2029
2030        segmentlist=[]
2031        segmenttaglist=[]
2032        for seg in userSegments:
2033            segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
2034            segmenttaglist.append(seg.tag)
2035        meshDict['outline_segments'] =segmentlist
2036        meshDict['outline_segment_tags'] =segmenttaglist
2037       
2038        holelist=[]
2039        for hole in holes:
2040            holelist.append([hole.x,hole.y]) 
2041        meshDict['holes'] = holelist
2042       
2043        regionlist=[]
2044        regiontaglist = []
2045        regionmaxarealist = []
2046        for region in regions:
2047            regionlist.append([region.x,region.y])
2048            regiontaglist.append(region.getTag())
2049           
2050            if (region.getMaxArea() != None): 
2051                regionmaxarealist.append(region.getMaxArea())
2052            else:
2053                regionmaxarealist.append(NOMAXAREA)
2054        meshDict['regions'] = regionlist
2055        meshDict['region_tags'] = regiontaglist
2056        meshDict['region_max_areas'] = regionmaxarealist
2057        #print "*(*("
2058        #print meshDict
2059        #print meshDict['regionlist']
2060        #print "*(*("
2061        return meshDict
2062
2063    def IOTriangulation2Mesh(self, genDict):
2064        """
2065        Set the mesh attributes given an tsh IO dictionary
2066        """
2067        #Clear the current generated mesh values
2068        self.tri_mesh = None
2069       
2070        self.tri_mesh = Rigid_triangulation(
2071            genDict['triangles']
2072            ,genDict['segments']
2073            ,genDict['vertices']
2074            ,genDict['triangle_tags']
2075            ,genDict['triangle_neighbors']
2076            ,genDict['segment_tags']
2077            ,genDict['vertex_attributes']
2078            ,genDict['vertex_attribute_titles']
2079            )
2080        self.attributeTitles = genDict['vertex_attribute_titles']
2081        self.maxVertexIndex = len(genDict['vertices'])
2082        #print "self.maxVertexIndex ", self.maxVertexIndex
2083       
2084    def IOOutline2Mesh(self, genDict):
2085        """
2086        Set the outline (user Mesh attributes) given a IO tsh dictionary
2087       
2088        mesh is an instance of a mesh object
2089        """
2090        #Clear the current user mesh values
2091        self.clearUserSegments()
2092        self.userVertices=[]
2093        self.Holes=[]
2094        self.Regions=[]
2095
2096        #print "mesh.IOOutline2Mesh@#@#@#"
2097        #print "genDict",genDict
2098        #print "@#@#@#"
2099       
2100        #index = 0
2101        for point in genDict['points']:
2102            v=Vertex(point[0], point[1])
2103            #v.index = index
2104            #index +=1
2105            self.userVertices.append(v)
2106
2107        #index = 0
2108        for seg,tag in map(None,genDict['outline_segments'],
2109                           genDict['outline_segment_tags']):
2110
2111            segObject = Segment( self.userVertices[int(seg[0])],
2112                           self.userVertices[int(seg[1])], tag = tag )
2113            #segObject.index = index
2114            #index +=1
2115            self.userSegments.append(segObject)
2116
2117# Remove the loading of attribute info.
2118# Have attribute info added using least_squares in pyvolution
2119#         index = 0
2120#         for att in genDict['point_attributes']:
2121#             if att == None:
2122#                 self.userVertices[index].setAttributes([])
2123#             else:
2124#                 self.userVertices[index].setAttributes(att)
2125#            index += 1
2126       
2127        #index = 0
2128        for point in genDict['holes']:
2129            h=Hole(point[0], point[1])
2130            #h.index = index
2131            #index +=1
2132            self.holes.append(h)
2133
2134        #index = 0
2135        for reg,att,maxArea in map(None,
2136                                   genDict['regions'],
2137                                   genDict['region_tags'],
2138                                   genDict['region_max_areas']):
2139            if maxArea > 0:  # maybe I should ref NOMAXAREA? Prob' not though
2140                Object = Region( reg[0],
2141                                 reg[1],
2142                                 tag = att,
2143                                 maxArea = maxArea)
2144            else:
2145                Object = Region( reg[0],
2146                                 reg[1],
2147                                 tag = att)
2148               
2149            #Object.index = index
2150            #index +=1
2151            self.regions.append(Object)
2152 
2153def importUngenerateFile(ofile):
2154    """
2155    import a file, ofile, with the format
2156    [poly]
2157    poly format:
2158    First line:  <# of vertices> <x centroid> <y centroid>
2159    Following lines: <x> <y>
2160    last line:  "END"
2161
2162    Note: These are clockwise.
2163    """
2164    fd = open(ofile,'r')
2165    Dict = readUngenerateFile(fd)
2166    fd.close()
2167    return Dict
2168
2169def readUngenerateFile(fd):
2170    """
2171    import a file, ofile, with the format
2172    [poly]
2173    poly format:
2174    First line:  <# of polynomial> <x centroid> <y centroid>
2175    Following lines: <x> <y>
2176    last line:  "END"
2177    """
2178   
2179    END_DELIMITER = 'END'
2180   
2181    points = []
2182    segments = []
2183    polygons = []
2184   
2185    isEnd = False
2186    line = fd.readline() #not used <# of polynomial> <x> <y>
2187    while not isEnd:
2188        poly = []
2189        line = fd.readline()
2190        fragments = line.split()
2191        x = float(fragments.pop(0))
2192        y = float(fragments.pop(0))
2193        points.append([x,y])
2194        poly.append([x,y])
2195        PreviousVertIndex = len(points)-1
2196        firstVertIndex = PreviousVertIndex
2197       
2198        line = fd.readline() #Read the next line
2199        while not line.startswith(END_DELIMITER): 
2200            #print "line >" + line + "<"
2201            fragments = line.split()
2202            x = float(fragments.pop(0))
2203            y = float(fragments.pop(0))
2204            points.append([x,y])
2205            poly.append([x,y])
2206            thisVertIndex = len(points)-1
2207            segment = [PreviousVertIndex,thisVertIndex]
2208            segments.append(segment)
2209            PreviousVertIndex = thisVertIndex
2210            line = fd.readline() #Read the next line
2211        # If the last and first segments are the same,
2212        # Remove the last segment and the last vertex
2213        # then add a segment from the second last vert to the 1st vert
2214        thisVertIndex = len(points)-1
2215        firstVert = points[firstVertIndex]
2216        thisVert = points[thisVertIndex]
2217        #print "firstVert",firstVert
2218        #print "thisVert",thisVert
2219        if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
2220            points.pop()
2221            segments.pop()
2222            poly.pop()
2223            thisVertIndex = len(points)-1
2224        segments.append([thisVertIndex, firstVertIndex])
2225       
2226        line = fd.readline() # read <# of polynomial> <x> <y> OR END
2227        #print "line >>" + line + "<<"
2228        # do poly stuff here
2229        polygons.append(poly)
2230        if line.startswith(END_DELIMITER):
2231            isEnd = True
2232   
2233    #print "points", points       
2234    #print "segments", segments
2235    ungenerated_dict = {}
2236    ungenerated_dict['points'] = points
2237    ungenerated_dict['segments'] = segments
2238    ungenerated_dict['polygons'] = polygons
2239    return ungenerated_dict
2240
2241def importMeshFromFile(ofile):
2242    """returns a mesh object, made from a points file or .tsh/.msh file
2243    Often raises IOError,RuntimeError
2244    """
2245    newmesh = None
2246    if (ofile[-4:]== ".pts" or ofile[-4:]== ".txt" or \
2247        ofile[-4:]== ".csv"):
2248        geospatial = Geospatial_data(ofile)
2249        dict = {}
2250        dict['points'] = geospatial.get_data_points(absolute=False)
2251        dict['outline_segments'] = []
2252        dict['outline_segment_tags'] = []
2253        dict['regions'] = []
2254        dict['region_tags'] = []
2255        dict['region_max_areas'] = []
2256        dict['holes'] = [] 
2257        newmesh= Mesh(geo_reference = geospatial.geo_reference)
2258        newmesh.IOOutline2Mesh(dict)
2259        counter = newmesh.removeDuplicatedUserVertices()
2260        if (counter >0):
2261            log.critical("%i duplicate vertices removed from dataset" % counter)
2262    elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
2263        dict = import_mesh_file(ofile)
2264        #print "********"
2265        #print "zq mesh.dict",dict
2266        #print "********"
2267        newmesh= Mesh()
2268        newmesh.IOOutline2Mesh(dict)
2269        newmesh.IOTriangulation2Mesh(dict)
2270    else:
2271        raise RuntimeError
2272   
2273    if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
2274        newmesh.geo_reference = dict['geo_reference']
2275    return newmesh
2276
2277def loadPickle(currentName):
2278    fd = open(currentName)
2279    mesh = pickle.load(fd)
2280    fd.close()
2281    return mesh
2282   
2283def square_outline(side_length = 1,up = "top", left = "left", right = "right",
2284                   down = "bottom", regions = False):
2285   
2286        a = Vertex (0,0)
2287        b = Vertex (0,side_length)
2288        c = Vertex (side_length,0)
2289        d = Vertex (side_length,side_length)
2290     
2291        s2 = Segment(b,d, tag = up)
2292        s3 = Segment(b,a, tag = left)
2293        s4 = Segment(d,c, tag = right)
2294        s5 = Segment(a,c, tag = down)
2295
2296        if regions:
2297            e =  Vertex (side_length/2,side_length/2)
2298            s6 = Segment(a,e, tag = down + left)
2299            s7 = Segment(b,e, tag = up + left)
2300            s8 = Segment(c,e, tag = down + right)
2301            s9 = Segment(d,e, tag = up + right)
2302            r1 = Region(side_length/2,3.*side_length/4, tag = up)
2303            r2 = Region(1.*side_length/4,side_length/2, tag = left)
2304            r3 = Region(3.*side_length/4,side_length/2, tag = right)
2305            r4 = Region(side_length/2,1.*side_length/4, tag = down)
2306            mesh = Mesh(userVertices=[a,b,c,d,e],
2307                        userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
2308                        regions = [r1,r2,r3,r4])
2309        else:
2310            mesh = Mesh(userVertices=[a,b,c,d],
2311                 userSegments=[s2,s3,s4,s5])
2312     
2313        return mesh
2314
2315   
2316
2317def region_strings2ints(region_list):
2318    """Given a list of (x_int,y_int,tag_string) lists it returns a list of
2319    (x_int,y_int,tag_int) and a list to convert the tag_int's back to
2320    the tag_strings 
2321    """
2322    # Make sure "" has an index of 0
2323    region_list.reverse()
2324    region_list.append((1.0,2.0,""))
2325    region_list.reverse()
2326    convertint2string = []
2327    for i in xrange(len(region_list)):
2328        convertint2string.append(region_list[i][2])
2329        if len(region_list[i]) == 4: # there's an area value
2330            region_list[i] = (region_list[i][0],
2331                              region_list[i][1],i,region_list[i][3])
2332        elif len(region_list[i]) == 3: # no area value
2333            region_list[i] = (region_list[i][0],region_list[i][1],i)
2334        else:
2335            log.critical("The region list has a bad size")
2336            # raise an error ..
2337            raise Error
2338
2339    #remove "" from the region_list
2340    region_list.pop(0)
2341   
2342    return [region_list, convertint2string]
2343
2344def region_ints2strings(region_list,convertint2string):
2345    """Reverses the transformation of region_strings2ints
2346    """
2347    #print 'region_ints2strings region_list', region_list
2348   
2349    returned_region_list = []
2350    # may not need (not region_list[0] == [])
2351    # or region_list[0] == [0.0]
2352    if (not region_list[0] == []): # or region_list[0] == [0.0]:
2353        #print "in loop"
2354        for i in xrange(len(region_list)):
2355            temp = region_list[i]
2356            returned_region_list.append(convertint2string[int(temp[0])])
2357    return returned_region_list
2358
2359def segment_ints2strings(intlist, convertint2string):
2360    """Reverses the transformation of segment_strings2ints """
2361    stringlist = []
2362    for x in intlist:
2363        stringlist.append(convertint2string[int(x)])
2364    return stringlist
2365
2366def segment_strings2ints(stringlist, preset):
2367    """Given a list of strings return a list of 0 to n ints which represent
2368    the strings and a converting list of the strings, indexed by 0 to n ints.
2369    Also, input an initial converting list of the strings
2370    Note, the converting list starts off with
2371    ["internal boundary", "external boundary", "internal boundary"]
2372    example input and output
2373    input ["a","b","a","c"],["c"]
2374    output [[2, 1, 2, 0], ['c', 'b', 'a']]
2375
2376    the first element in the converting list will be
2377    overwritten with "".
2378    ?This will become the third item in the converting list?
2379   
2380    # note, order the initial converting list is important,
2381     since the index = the int tag
2382    """
2383    nodups = unique(stringlist)
2384    # note, order is important, the index = the int tag
2385    #preset = ["internal boundary", "external boundary"]
2386    #Remove the preset tags from the list with no duplicates
2387    nodups = [x for x in nodups if x not in preset]
2388
2389    try:
2390        nodups.remove("") # this has to go to zero
2391    except ValueError:
2392        pass
2393   
2394    # Add the preset list at the beginning of no duplicates
2395    preset.reverse()
2396    nodups.extend(preset)
2397    nodups.reverse()
2398
2399       
2400    convertstring2int = {}
2401    convertint2string = []
2402    index = 0
2403    for x in nodups:
2404        convertstring2int[x] = index
2405        convertint2string.append(x)
2406        index += 1
2407    convertstring2int[""] = 0
2408
2409    intlist = []
2410    for x in stringlist:
2411        intlist.append(convertstring2int[x])
2412    return [intlist, convertint2string]
2413
2414
2415def unique(s):
2416    """Return a list of the elements in s, but without duplicates.
2417
2418    For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
2419    unique("abcabc") some permutation of ["a", "b", "c"], and
2420    unique(([1, 2], [2, 3], [1, 2])) some permutation of
2421    [[2, 3], [1, 2]].
2422
2423    For best speed, all sequence elements should be hashable.  Then
2424    unique() will usually work in linear time.
2425
2426    If not possible, the sequence elements should enjoy a total
2427    ordering, and if list(s).sort() doesn't raise TypeError it's
2428    assumed that they do enjoy a total ordering.  Then unique() will
2429    usually work in O(N*log2(N)) time.
2430
2431    If that's not possible either, the sequence elements must support
2432    equality-testing.  Then unique() will usually work in quadratic
2433    time.
2434    """
2435
2436    n = len(s)
2437    if n == 0:
2438        return []
2439
2440    # Try using a dict first, as that's the fastest and will usually
2441    # work.  If it doesn't work, it will usually fail quickly, so it
2442    # usually doesn't cost much to *try* it.  It requires that all the
2443    # sequence elements be hashable, and support equality comparison.
2444    u = {}
2445    try:
2446        for x in s:
2447            u[x] = 1
2448    except TypeError:
2449        del u  # move on to the next method
2450    else:
2451        return u.keys()
2452
2453    # We can't hash all the elements.  Second fastest is to sort,
2454    # which brings the equal elements together; then duplicates are
2455    # easy to weed out in a single pass.
2456    # NOTE:  Python's list.sort() was designed to be efficient in the
2457    # presence of many duplicate elements.  This isn't true of all
2458    # sort functions in all languages or libraries, so this approach
2459    # is more effective in Python than it may be elsewhere.
2460    try:
2461        t = list(s)
2462        t.sort()
2463    except TypeError:
2464        del t  # move on to the next method
2465    else:
2466        assert n > 0
2467        last = t[0]
2468        lasti = i = 1
2469        while i < n:
2470            if t[i] != last:
2471                t[lasti] = last = t[i]
2472                lasti += 1
2473            i += 1
2474        return t[:lasti]
2475
2476    # Brute force is all that's left.
2477    u = []
2478    for x in s:
2479        if x not in u:
2480            u.append(x)
2481    return u
2482
2483
2484# FIXME (DSG-DSG)
2485# To do-
2486# Create a clear interface. eg
2487# have the interface methods more at the top of this file and add comments
2488# for the interface functions/methods, use function_name (not functionName),
2489
2490#Currently
2491#function_name methods assume absolute values.  Geo-refs can be passed in.
2492#
2493
2494# instead of functionName
2495if __name__ == "__main__":
2496    pass
Note: See TracBrowser for help on using the repository browser.