Changeset 6074


Ignore:
Timestamp:
Dec 12, 2008, 1:22:18 PM (15 years ago)
Author:
rwilson
Message:

Looked at NetCDF code, also did a pep8.

Location:
anuga_core/source/anuga/load_mesh
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/load_mesh/loadASCII.py

    r4903 r6074  
    11"""
    2 
    32The format for a Points dictionary is:
    43
    … …  
    1211    dic['attributelist']['elevation'] = [[7.0,5.0]
    1312
    14    
     13
    1514    The dict format for IO with mesh files is;
    1615    (the triangulation)
    … …  
    1817    vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
    1918    vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
    20     segments: [[v1,v2],[v3,v4],...] (lists of integers) 
     19    segments: [[v1,v2],[v3,v4],...] (lists of integers)
    2120    segment_tags : [tag,tag,...] list of strings
    2221    triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
    23     triangle_tags: [s1,s2,...] A list of strings 
     22    triangle_tags: [s1,s2,...] A list of strings
    2423    triangle_neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
    25        
    26     (the outline)   
     24
     25    (the outline)
    2726    points: [[x1,y1],[x2,y2],...] (lists of doubles)
    2827    point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
    29     outline_segments: [[point1,point2],[p3,p4],...] (lists of integers) 
     28    outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
    3029    outline_segment_tags : [tag1,tag2,...] list of strings
    3130    holes : [[x1,y1],...](List of doubles, one inside each hole region)
    … …  
    4443
    4544    The format for a .tsh file is (FIXME update this)
    46    
     45
    4746    First line:  <# of vertices> <# of attributes>
    4847    Following lines:  <vertex #> <x> <y> [attributes]
    49     One line:  <# of triangles>
    50     Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
    51     One line:  <# of segments>
     48    One line:  <# of triangles>
     49    Following lines: <triangle #> <vertex #>  <vertex #> <vertex #>
     50                         <neigbouring triangle #> <neigbouring triangle #>
     51                         <neigbouring triangle #> [attribute of region]
     52    One line:  <# of segments>
    5253    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
    5354"""
    … …  
    5758
    5859from string import  find, rfind
    59 from Numeric import array, Float, Int16, Int32, Character,reshape, concatenate, take
     60from Numeric import array, Float, Int16, Int32, Character,reshape, \
     61                    concatenate, take
    6062from os.path import splitext
    6163
    62 from anuga.coordinate_transforms.geo_reference import Geo_reference,TITLE, TitleError
     64from anuga.coordinate_transforms.geo_reference import Geo_reference, TITLE, \
     65                                                      TitleError
    6366
    6467from Scientific.IO.NetCDF import NetCDFFile
    65    
     68
    6669import exceptions
    6770class TitleAmountError(exceptions.Exception): pass
    6871
     72
    6973NOMAXAREA=-999
    7074
    71 # This is used by pmesh
     75
     76##
     77# @brief Read a mesh file (.tsh or .msh) and return a dictionary.
     78# @param ofile Path to the file to read.
     79# @return A dictionary of data from the input file.
     80# @note Throws IOError if can't read file.
     81# @note This is used by pmesh.
    7282def import_mesh_file(ofile):
    7383    """
    7484    read a mesh file, either .tsh or .msh
    75    
     85
    7686    Note: will throw an IOError if it can't load the file.
    7787    Catch these!
    … …  
    8292        elif ofile[-4:]== ".msh":
    8393            dict = _read_msh_file(ofile)
    84         else:     
    85             msg = 'Extension .%s is unknown' %ofile[-4:]
     94        else:
     95            msg = 'Extension .%s is unknown' % ofile[-4:]
    8696            raise IOError, msg
    87     except (TitleError,SyntaxError,IndexError, ValueError): #FIXME No test for ValueError
     97    #FIXME No test for ValueError
     98    except (TitleError, SyntaxError, IndexError, ValueError):
    8899            msg = 'File could not be opened'
    89             raise IOError, msg 
     100            raise IOError, msg
    90101    return dict
    91102
    92103
    93 def export_mesh_file(ofile,mesh_dict):
     104##
     105# @brief Write a mesh file from a dictionary.
     106# @param ofile The path of the mesh file to write.
     107# @param mesh_dict Dictionary containing data to write to output mesh file.
     108# @note Raises IOError if file extension is unrecognized.
     109def export_mesh_file(ofile, mesh_dict):
    94110    """
    95111    write a file, ofile, with the format
    96    
    97     First line:  <# of vertices> <# of attributes>
     112
     113    First line:       <# of vertices> <# of attributes>
    98114    Following lines:  <vertex #> <x> <y> [attributes]
    99     One line:  <# of triangles>
    100     Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
    101     One line:  <# of segments>
     115    One line:         <# of triangles>
     116    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #>
     117                          <neigbouring triangle #> <neigbouring triangle #>
     118                          <neigbouring triangle #> [attribute of region]
     119    One line:         <# of segments>
    102120    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
    103121    """
    104     #FIXME DSG-anyone: automate
    105     if not mesh_dict.has_key('points'):
    106         mesh_dict['points'] = []
    107     if not mesh_dict.has_key('point_attributes'):
    108         mesh_dict['point_attributes'] = []
    109     if not mesh_dict.has_key('outline_segments'):
    110         mesh_dict['outline_segments'] = []
    111     if not mesh_dict.has_key('outline_segment_tags'):
    112         mesh_dict['outline_segment_tags'] = []
    113     if not mesh_dict.has_key('holes'):
    114         mesh_dict['holes'] = []
    115     if not mesh_dict.has_key('regions'):
    116         mesh_dict['regions'] = []
    117        
    118     if not mesh_dict.has_key('region_tags'):
    119         mesh_dict['region_tags'] = []
    120     if not mesh_dict.has_key('region_max_areas'):
    121         mesh_dict['region_max_areas'] = []
    122     if not mesh_dict.has_key('vertices'):
    123         mesh_dict['vertices'] = []
    124     if not mesh_dict.has_key('vertex_attributes'):
    125         mesh_dict['vertex_attributes'] = []
    126     if not mesh_dict.has_key('vertex_attribute_titles'):
    127         mesh_dict['vertex_attribute_titles'] = []
    128     if not mesh_dict.has_key('segments'):
    129         mesh_dict['segments'] = []
    130     if not mesh_dict.has_key('segment_tags'):
    131         mesh_dict['segment_tags'] = []
    132     if not mesh_dict.has_key('triangles'):
    133         mesh_dict['triangles'] = []
    134     if not mesh_dict.has_key('triangle_tags'):
    135         mesh_dict['triangle_tags'] = []
    136     if not mesh_dict.has_key('triangle_neighbors'):
    137         mesh_dict['triangle_neighbors'] = []
    138     #print "DSG************"
    139     #print "mesh_dict",mesh_dict
    140     #print "DSG************"
    141    
     122
     123    # Ensure that if required key isn't present, we add it with [] value.
     124    # .setdefault() for dictionaries was added in python 2.0.
     125    # The only other neat way to do this is to override the dictionary
     126    # .__getitem__() method to return [] if key not found.
     127    reqd_keys = ['points', 'point_attributes', 'outline_segments',
     128                 'outline_segment_tags', 'holes', 'regions', 'region_tags',
     129                 'region_max_areas', 'vertices', 'vertex_attributes',
     130                 'vertex_attribute_titles', 'segments', 'segment_tags',
     131                 'triangles', 'triangle_tags', 'triangle_neighbors']
     132
     133    for k in reqd_keys:
     134        mesh_dict.setdefault(k, [])
     135
     136    # hand-off to appropriate function
    142137    if (ofile[-4:] == ".tsh"):
    143138        _write_tsh_file(ofile,mesh_dict)
    … …  
    146141    else:
    147142        msg = 'Unknown file type %s ' %ofile
    148         raise IOError, msg
    149 
     143        raise IOError, msg
     144
     145
     146##
     147# @brief Read .tsh file into a dictionary.
     148# @param ofile Path of the file to read.
     149# @return Data dictionary from the .tsh file.
    150150def _read_tsh_file(ofile):
    151151    """
    152152    Read the text file format for meshes
    153153    """
    154     fd = open(ofile,'r')
     154
     155    fd = open(ofile, 'r')
    155156    dict = _read_triangulation(fd)
    156157    dict_mesh = _read_outline(fd)
    … …  
    158159        dict[element] = dict_mesh[element]
    159160    fd.close()
     161
    160162    return dict
    161163
    162    
     164
     165##
     166# @brief ??
     167# @param fd An open descriptor of the file to read.
     168# @return A dictionary ...
    163169def _read_triangulation(fd):
    164170    """
    165     Read the generated triangulation, NOT the outline.
    166     """
     171    Read the generated triangulation, NOT the outline.
     172    """
     173
    167174    delimiter = " "
    168     ######### loading the point info
     175
     176    # loading the point info
    169177    line = fd.readline()
    170     #print line
    171178    fragments = line.split()
    172     #for fragment in fragments:
    173     #print fragment
    174179    if fragments ==[]:
    175180        NumOfVertices = 0
    … …  
    180185    points = []
    181186    pointattributes = []
    182     for index in range(int(NumOfVertices)):       
    183         #print index
     187    for index in range(int(NumOfVertices)):
    184188        fragments = fd.readline().split()
    185         #print fragments
    186189        fragments.pop(0) #pop off the index
    187        
    188         # pop the x y off so we're left with a list of attributes       
    189         vert = [float(fragments.pop(0)),float(fragments.pop(0))]
     190
     191        # pop the x y off so we're left with a list of attributes
     192        vert = [float(fragments.pop(0)), float(fragments.pop(0))]
    190193        points.append(vert)
    191194        apointattributes  = []
    192         #print fragments
    193195        for fragment in fragments:
    194196            apointattributes.append(float(fragment))
    195197        if apointattributes != []:
    196198            pointattributes.append(apointattributes)
    197        
    198     ######### loading the point title info
     199
     200    # loading the point title info
    199201    line = fd.readline()
    200     #print "point title comments",line
    201202    vertTitle = []
    202     for index in range(int(NumOfVertAttributes)):       
    203         #print index
    204         fragments = fd.readline().strip()
     203    for index in range(int(NumOfVertAttributes)):
     204        fragments = fd.readline().strip()
    205205        vertTitle.append(fragments)
    206         #print vertTitle
    207        
    208     ######### loading the triangle info
     206
     207    # loading the triangle info
    209208    line = fd.readline()
    210     #print "triangle comments",line
    211209    fragments = line.split()
    212     #for fragment in fragments:
    213     #    print fragment
    214210    NumOfTriangles = fragments[0]
    215211    triangles = []
    216212    triangleattributes = []
    217213    triangleneighbors = []
    218     for index in range(int(NumOfTriangles)):       
    219         #print index
     214    for index in range(int(NumOfTriangles)):
    220215        line = fd.readline()
    221216        line.strip() # so we can get the region string
    222217        fragments = line.split()
    223         #print "triangle info", fragments
    224218        fragments.pop(0) #pop off the index
    225        
    226         tri = [int(fragments[0]),int(fragments[1]),int(fragments[2])]
     219
     220        tri = [int(fragments[0]), int(fragments[1]), int(fragments[2])]
    227221        triangles.append(tri)
    228         neighbors = [int(fragments[3]),int(fragments[4]),int(fragments[5])]
     222        neighbors = [int(fragments[3]), int(fragments[4]), int(fragments[5])]
    229223        triangleneighbors.append(neighbors)
    230224        for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
    231             line = line[find(line,delimiter):] # remove index
     225            line = line[find(line, delimiter):] # remove index
    232226            line = line.lstrip()
    233227        stringtag = line.strip()
    234228        triangleattributes.append(stringtag)
    235        
    236     ######### loading the segment info
     229
     230    # loading the segment info
    237231    line = fd.readline()
    238     #print "seg comment line",line
    239232    fragments = line.split()
    240     #for fragment in fragments:
    241     #    print fragment
    242233    NumOfSegments = fragments[0]
    243234    segments = []
    244235    segmenttags = []
    245     for index in range(int(NumOfSegments)):       
    246         #print index
     236    for index in range(int(NumOfSegments)):
    247237        line = fd.readline()
    248238        line.strip() # to get the segment string
    249239        fragments = line.split()
    250         #print fragments
    251240        fragments.pop(0) #pop off the index
    252         seg = [int(fragments[0]),int(fragments[1])]
     241        seg = [int(fragments[0]), int(fragments[1])]
    253242        segments.append(seg)
    254         line = line[find(line,delimiter):] # remove index
     243        line = line[find(line, delimiter):] # remove index
    255244        line = line.lstrip()
    256         line = line[find(line,delimiter):] # remove x
     245        line = line[find(line, delimiter):] # remove x
    257246        line = line.lstrip()
    258         line = line[find(line,delimiter):] # remove y
     247        line = line[find(line, delimiter):] # remove y
    259248        stringtag = line.strip()
    260249        segmenttags.append(stringtag)
    261250
    262            
    263251    meshDict = {}
    264252    meshDict['vertices'] = points
    … …  
    269257    meshDict['triangles'] = triangles
    270258    meshDict['triangle_tags'] = triangleattributes
    271     meshDict['triangle_neighbors'] = triangleneighbors 
    272     meshDict['segments'] = segments 
     259    meshDict['triangle_neighbors'] = triangleneighbors
     260    meshDict['segments'] = segments
    273261    meshDict['segment_tags'] = segmenttags
    274262    meshDict['vertex_attribute_titles'] = vertTitle
    275    
     263
    276264    return meshDict
    277    
     265
     266
     267##
     268# @brief
     269# @param fd An open descriptor of the file to read.
     270# @return A dictionary ...
     271# @note If file has no mesh info, an empty dict will be returned.
    278272def _read_outline(fd):
    279273    """
    … …  
    281275    returned will be 'empty'.
    282276    """
     277
    283278    delimiter = " " # warning: split() calls are using default whitespace
    284    
    285     ######### loading the point info
     279
     280    # loading the point info
    286281    line = fd.readline()
    287     #print 'read_outline - line',line
    288282    fragments = line.split()
    289     #for fragment in fragments:
    290283    if fragments ==[]:
    291284        NumOfVertices = 0
    … …  
    296289    points = []
    297290    pointattributes = []
    298     for index in range(int(NumOfVertices)):       
    299         #print index
    300         fragments = fd.readline().split()
    301         #print fragments
     291    for index in range(int(NumOfVertices)):
     292        fragments = fd.readline().split()
    302293        fragments.pop(0) #pop off the index
    303294        # pop the x y off so we're left with a list of attributes
    304         vert = [float(fragments.pop(0)),float(fragments.pop(0))]
     295        vert = [float(fragments.pop(0)), float(fragments.pop(0))]
    305296        points.append(vert)
    306         apointattributes  = []
    307         #print fragments
     297        apointattributes = []
    308298        for fragment in fragments:
    309299            apointattributes.append(float(fragment))
    310300        pointattributes.append(apointattributes)
    311        
    312 
    313     ######### loading the segment info
     301
     302    # loading the segment info
    314303    line = fd.readline()
    315304    fragments = line.split()
    316     #for fragment in fragments:
    317     #    print fragment
    318305    if fragments ==[]:
    319306        NumOfSegments = 0
    … …  
    322309    segments = []
    323310    segmenttags = []
    324     for index in range(int(NumOfSegments)): 
    325         #print index
     311    for index in range(int(NumOfSegments)):
    326312        line = fd.readline()
    327313        fragments = line.split()
    328         #print fragments
    329314        fragments.pop(0) #pop off the index
    330        
    331         seg = [int(fragments[0]),int(fragments[1])]
     315        seg = [int(fragments[0]), int(fragments[1])]
    332316        segments.append(seg)
    333         line = line[find(line,delimiter):] # remove index
     317        line = line[find(line, delimiter):] # remove index
    334318        line = line.lstrip()
    335         line = line[find(line,delimiter):] # remove x
     319        line = line[find(line, delimiter):] # remove x
    336320        line = line.lstrip()
    337         line = line[find(line,delimiter):] # remove y
     321        line = line[find(line, delimiter):] # remove y
    338322        stringtag = line.strip()
    339         segmenttags.append(stringtag) 
    340 
    341     ######### loading the hole info
     323        segmenttags.append(stringtag)
     324
     325    # loading the hole info
    342326    line = fd.readline()
    343     #print line
    344327    fragments = line.split()
    345     #for fragment in fragments:
    346     #    print fragment
    347     if fragments ==[]:
     328    if fragments == []:
    348329        numOfHoles = 0
    349330    else:
    350331        numOfHoles = fragments[0]
    351332    holes = []
    352     for index in range(int(numOfHoles)):       
    353         #print index
     333    for index in range(int(numOfHoles)):
    354334        fragments = fd.readline().split()
    355         #print fragments
    356335        fragments.pop(0) #pop off the index
    357         hole = [float(fragments[0]),float(fragments[1])]
     336        hole = [float(fragments[0]), float(fragments[1])]
    358337        holes.append(hole)
    359338
    360    
    361     ######### loading the region info
     339    # loading the region info
    362340    line = fd.readline()
    363     #print line
    364341    fragments = line.split()
    365     #for fragment in fragments:
    366     #    print fragment
    367     if fragments ==[]:
     342    if fragments == []:
    368343        numOfRegions = 0
    369344    else:
    … …  
    374349        line = fd.readline()
    375350        fragments = line.split()
    376         #print fragments
    377351        fragments.pop(0) #pop off the index
    378         region = [float(fragments[0]),float(fragments[1])]
     352        region = [float(fragments[0]), float(fragments[1])]
    379353        regions.append(region)
    380354
    381         line = line[find(line,delimiter):] # remove index
     355        line = line[find(line, delimiter):] # remove index
    382356        line = line.lstrip()
    383         line = line[find(line,delimiter):] # remove x
     357        line = line[find(line, delimiter):] # remove x
    384358        line = line.lstrip()
    385         line = line[find(line,delimiter):] # remove y
     359        line = line[find(line, delimiter):] # remove y
    386360        stringtag = line.strip()
    387361        regionattributes.append(stringtag)
     362
    388363    regionmaxareas = []
    389364    line = fd.readline()
    … …  
    391366        line = fd.readline()
    392367        fragments = line.split()
    393         #print fragments
    394368        # The try is here for format compatibility
    395369        try:
    … …  
    405379        geo_reference = Geo_reference(ASCIIFile=fd)
    406380    except:
    407         #geo_ref not compulsory 
     381        #geo_ref not compulsory
    408382        geo_reference = None
    409        
     383
    410384    meshDict = {}
    411385    meshDict['points'] = points
    412386    meshDict['point_attributes'] = pointattributes
    413     meshDict['outline_segments'] = segments 
     387    meshDict['outline_segments'] = segments
    414388    meshDict['outline_segment_tags'] = segmenttags
    415389    meshDict['holes'] = holes
    … …  
    418392    meshDict['region_max_areas'] = regionmaxareas
    419393    meshDict['geo_reference'] = geo_reference
    420    
    421     #print "meshDict",meshDict
     394
    422395    return meshDict
    423396
    424 def clean_line(line,delimiter):     
     397
     398##
     399# @brief Clean up a line with delimited fields.
     400# @param line The string to be cleaned.
     401# @param delimiter String used a delimiter when splitting 'line'.
     402# @return A list of delimited fields with white-space removed from ends.
     403# @note Will also remove any field that was originally ''.
     404def clean_line(line, delimiter):
    425405    """Remove whitespace
    426406    """
    427     #print ">%s" %line
     407
    428408    line = line.strip()
    429     #print "stripped>%s" %line
    430409    numbers = line.split(delimiter)
    431410    i = len(numbers) - 1
    … …  
    435414        else:
    436415            numbers[i] = numbers[i].strip()
    437        
     416
    438417        i += -1
    439     #for num in numbers:
    440     #    print "num>%s<" %num
     418
    441419    return numbers
    442420
    443 def _write_ASCII_triangulation(fd,
    444                                gen_dict):
     421
     422##
     423# @brief Write an ASCII triangulation file.
     424# @param fd An open descriptor to the file to write.
     425# @param gen_dict Dictionary containing data to write.
     426def _write_ASCII_triangulation(fd, gen_dict):
    445427    vertices = gen_dict['vertices']
    446428    vertices_attributes = gen_dict['vertex_attributes']
    447429    try:
    448         vertices_attribute_titles = gen_dict['vertex_attribute_titles']   
     430        vertices_attribute_titles = gen_dict['vertex_attribute_titles']
    449431    except KeyError, e:
    450         #FIXME is this the best way? 
     432        #FIXME is this the best way?
    451433        if vertices_attributes == [] or vertices_attributes[0] == []:
    452434             vertices_attribute_titles = []
    453435        else:
    454             raise KeyError, e 
    455        
     436            raise KeyError, e
     437
    456438    triangles = gen_dict['triangles']
    457439    triangles_attributes = gen_dict['triangle_tags']
    458440    triangle_neighbors = gen_dict['triangle_neighbors']
    459        
     441
    460442    segments = gen_dict['segments']
    461443    segment_tags = gen_dict['segment_tags']
    462      
     444
    463445    numVert = str(len(vertices))
    464     #print "load vertices_attributes", vertices_attributes
    465446    # Don't understand why we have to do vertices_attributes[0] is None,
    466447    # but it doesn't work otherwise...
    467     if vertices_attributes == None or \
    468            (numVert == "0") or \
    469            len(vertices_attributes) == 0:
     448    if vertices_attributes == None \
     449       or (numVert == "0") \
     450       or len(vertices_attributes) == 0:
    470451        numVertAttrib = "0"
    471452    else:
    472453        numVertAttrib = str(len(vertices_attributes[0]))
    473     fd.write(numVert + " " + numVertAttrib + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices..." + "\n")
    474 
    475     #<vertex #> <x> <y> [attributes]   
     454
     455    fd.write(numVert + " " + numVertAttrib +
     456             " # <# of verts> <# of vert attributes>, next lines <vertex #> "
     457             "<x> <y> [attributes] ...Triangulation Vertices...\n")
     458
     459    #<vertex #> <x> <y> [attributes]
    476460    index = 0
    477461    for vert in vertices:
    478462        attlist = ""
    479        
    480         if vertices_attributes == None or \
    481                vertices_attributes == []:
     463
     464        if vertices_attributes == None or vertices_attributes == []:
    482465            attlist = ""
    483466        else:
    484467            for att in vertices_attributes[index]:
    485                 attlist = attlist + str(att)+" "
     468                attlist = attlist + str(att) + " "
    486469        attlist.strip()
    487         fd.write(str(index) + " "
    488                  + str(vert[0]) + " "
    489                  + str(vert[1]) + " "
    490                  + attlist + "\n")
     470        fd.write(str(index) + " " + str(vert[0]) + " " + str(vert[1])
     471                 + " " + attlist + "\n")
    491472        index += 1
    492473
    493     # write comments for title 
    494     fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
     474    # write comments for title
     475    fd.write("# attribute column titles ...Triangulation Vertex Titles...\n")
    495476    for title in vertices_attribute_titles:
    496477        fd.write(title + "\n")
    497        
     478
    498479    #<# of triangles>
    499480    n = len(triangles)
    500     fd.write(str(n) + " # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
    501        
    502     # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
     481    fd.write(str(n)
     482             + " # <# of triangles>, next lines <triangle #> [<vertex #>] "
     483               "[<neigbouring triangle #>] [attribute of region] ..."
     484               "Triangulation Triangles...\n")
     485
     486    # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #>
     487    #    <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
    503488    for index in range(n):
    504489        neighbors = ""
    … …  
    521506        # with triangle, and it seems to have the option of returning
    522507        # more than one value for triangle attributex
    523         if triangles_attributes == None or \
    524                triangles_attributes == [] or \
    525                triangles_attributes[index] == ['']:
     508        if triangles_attributes == None \
     509           or triangles_attributes == [] \
     510           or triangles_attributes[index] == ['']:
    526511            att = ""
    527512        else:
    528513            att = str(triangles_attributes[index])
    529514        fd.write(str(index) + " "
    530                  + str(tri[0]) + " " 
    531                  + str(tri[1]) + " " 
    532                  + str(tri[2]) + " " 
     515                 + str(tri[0]) + " "
     516                 + str(tri[1]) + " "
     517                 + str(tri[2]) + " "
    533518                 + neighbors + " "
    534519                 + att + "\n")
    535            
    536     #One line:  <# of segments>
    537     fd.write(str(len(segments)) +
    538              " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments..." + "\n")
    539        
     520
     521    #One line:  <# of segments>
     522    fd.write(str(len(segments)) +
     523             " # <# of segments>, next lines <segment #> <vertex #>  "
     524             "<vertex #> [boundary tag] ...Triangulation Segments...\n")
     525
    540526    #Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
    541527    for i in range(len(segments)):
    542528        seg = segments[i]
    543529        fd.write(str(i) + " "
    544                  + str(seg[0]) + " " 
    545                  + str(seg[1]) + " " 
     530                 + str(seg[0]) + " "
     531                 + str(seg[1]) + " "
    546532                 + str(segment_tags[i]) + "\n")
    547533
    548534
    549 def _write_tsh_file(ofile,mesh_dict):
    550             fd = open(ofile,'w')
    551             _write_ASCII_triangulation(fd,mesh_dict)
    552             _write_ASCII_outline(fd,mesh_dict)   
    553             fd.close()
    554 
    555 def _write_ASCII_outline(fd,
    556                      dict):
     535##
     536# @brief Write triangulation and outline to a .tsh file.
     537# @param ofile Path of the file to write.
     538# @mesh_dict Dictionary of data to write.
     539def _write_tsh_file(ofile, mesh_dict):
     540    fd = open(ofile, 'w')
     541    _write_ASCII_triangulation(fd, mesh_dict)
     542    _write_ASCII_outline(fd, mesh_dict)
     543    fd.close()
     544
     545
     546##
     547# @brief Write an ASCII outline to a file.
     548# @param fd An open descriptor of the file to write.
     549# @param dict Dictionary holding data to write.
     550def _write_ASCII_outline(fd, dict):
    557551    points = dict['points']
    558552    point_attributes = dict['point_attributes']
    … …  
    563557    region_tags = dict['region_tags']
    564558    region_max_areas = dict['region_max_areas']
    565        
     559
    566560    num_points = str(len(points))
    567561    if (num_points == "0"):
    … …  
    569563    else:
    570564        num_point_atts = str(len(point_attributes[0]))
    571     fd.write(num_points + " " + num_point_atts + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
    572        
     565
     566    fd.write(num_points + " " + num_point_atts +
     567             " # <# of verts> <# of vert attributes>, next lines <vertex #> "
     568             "<x> <y> [attributes] ...Mesh Vertices...\n")
     569
    573570    # <x> <y> [attributes]
    574571    for i,point in enumerate(points):
    575572        attlist = ""
    576573        for att in point_attributes[i]:
    577             attlist = attlist + str(att)+" "
     574            attlist = attlist + str(att) + " "
    578575        attlist.strip()
    579576        fd.write(str(i) + " "
    580                  +str(point[0]) + " "
     577                 + str(point[0]) + " "
    581578                 + str(point[1]) + " "
    582579                 + attlist + "\n")
    583            
    584     #One line:  <# of segments>
    585     fd.write(str(len(segments)) +
    586              " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments..." + "\n")
    587        
     580
     581    #One line:  <# of segments>
     582    fd.write(str(len(segments)) +
     583             " # <# of segments>, next lines <segment #> <vertex #>  "
     584             "<vertex #> [boundary tag] ...Mesh Segments...\n")
     585
    588586    #Following lines:  <vertex #>  <vertex #> [boundary tag]
    589587    for i,seg in enumerate(segments):
    590588        fd.write(str(i) + " "
    591                  + str(seg[0]) + " " 
    592                  + str(seg[1]) + " " 
     589                 + str(seg[0]) + " "
     590                 + str(seg[1]) + " "
    593591                 + str(segment_tags[i]) + "\n")
    594592
    595     #One line:  <# of holes> 
    596     fd.write(str(len(holes)) + 
    597              " # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes..." + "\n")   
     593    #One line:  <# of holes>
     594    fd.write(str(len(holes)) +
     595             " # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes...\n")
    598596    # <x> <y>
    599597    for i,h in enumerate(holes):
    … …  
    601599                 + str(h[0]) + " "
    602600                 + str(h[1]) + "\n")
    603        
    604     #One line:  <# of regions>
    605     fd.write(str(len(regions)) +
    606              " # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")   
     601
     602    #One line:  <# of regions>
     603    fd.write(str(len(regions)) +
     604             " # <# of regions>, next lines <Region #> <x> <y> <tag>"
     605             "...Mesh Regions...\n")
     606
    607607    # <index> <x> <y> <tag>
    608     #print "regions",regions
    609608    for i,r in enumerate(regions):
    610         #print "r",r
    611         fd.write(str(i) + " "
     609        fd.write(str(i) + " "
    612610                 + str(r[0]) + " "
    613611                 + str(r[1])+ " "
    614612                 + str(region_tags[i]) + "\n")
    615        
     613
    616614    # <index> [<MaxArea>|""]
    617        
    618     #One line:  <# of regions>
    619     fd.write(str(len(regions)) +
    620              " # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions..." + "\n")
     615
     616    #One line:  <# of regions>
     617    fd.write(str(len(regions)) +
     618             " # <# of regions>, next lines <Region #> [Max Area] "
     619             "...Mesh Regions...\n")
    621620    for i,r in enumerate(regions):
    622621        area = str(region_max_areas[i])
    623            
     622
    624623        fd.write(str(i) + " " + area + "\n")
    625624
    … …  
    628627        dict['geo_reference'].write_ASCII(fd)
    629628
     629
     630##
     631# @brief Write a .msh file.
     632# @param file Path to the file to write.
     633# @param mesh Dictionary holding data to write.
    630634def _write_msh_file(file_name, mesh):
    631     """Write .msh NetCDF file   
     635    """Write .msh NetCDF file
    632636
    633637    WARNING: This function mangles the mesh data structure
    634638    """
    635  
     639
    636640    # FIXME(Ole and John): We ran into a problem on Bogong (64 bit)
    637641    # where integers appeared as arrays.  This may be similar to
    … …  
    643647    # created with the type Int64. Need to look at the NetCDF library
    644648    # in more detail.
    645    
     649
    646650    IntType = Int32
    647651    #IntType = Int
    648    
    649     #
     652
    650653    #the triangulation
    651654    mesh['vertices'] = array(mesh['vertices']).astype(Float)
    652655    if mesh['vertex_attributes'] != None:
    653656        mesh['vertex_attributes'] = \
    654                array(mesh['vertex_attributes']).astype(Float)
    655     mesh['vertex_attribute_titles'] = array(mesh['vertex_attribute_titles']).astype(Character)
     657            array(mesh['vertex_attributes']).astype(Float)
     658    mesh['vertex_attribute_titles'] = \
     659        array(mesh['vertex_attribute_titles']).astype(Character)
    656660    mesh['segments'] = array(mesh['segments']).astype(IntType)
    657661    mesh['segment_tags'] = array(mesh['segment_tags']).astype(Character)
    658662    mesh['triangles'] = array(mesh['triangles']).astype(IntType)
    659663    mesh['triangle_tags'] = array(mesh['triangle_tags']) #.astype(Character)
    660     mesh['triangle_neighbors'] = array(mesh['triangle_neighbors']).astype(IntType)
    661    
    662     #the outline
     664    mesh['triangle_neighbors'] = \
     665        array(mesh['triangle_neighbors']).astype(IntType)
     666
     667    #the outline
    663668    mesh['points'] = array(mesh['points']).astype(Float)
    664669    mesh['point_attributes'] = array(mesh['point_attributes']).astype(Float)
    665670    mesh['outline_segments'] = array(mesh['outline_segments']).astype(IntType)
    666     mesh['outline_segment_tags'] = array(mesh['outline_segment_tags']).astype(Character)
     671    mesh['outline_segment_tags'] = \
     672        array(mesh['outline_segment_tags']).astype(Character)
    667673    mesh['holes'] = array(mesh['holes']).astype(Float)
    668674    mesh['regions'] = array(mesh['regions']).astype(Float)
    669675    mesh['region_tags'] = array(mesh['region_tags']).astype(Character)
    670676    mesh['region_max_areas'] = array(mesh['region_max_areas']).astype(Float)
    671    
     677
    672678    #mesh = mesh_dict2array(mesh)
    673679    #print "new_mesh",new_mesh
    674     #print "mesh",mesh 
    675    
     680    #print "mesh",mesh
     681
    676682    # NetCDF file definition
    677683    try:
    678684        outfile = NetCDFFile(file_name, 'w')
    679685    except IOError:
    680         msg = 'File %s could not be created' %file_name
     686        msg = 'File %s could not be created' % file_name
    681687        raise msg
    682        
     688
    683689    #Create new file
    684690    outfile.institution = 'Geoscience Australia'
    685691    outfile.description = 'NetCDF format for compact and portable storage ' +\
    686                       'of spatial point data'
    687    
     692                          'of spatial point data'
     693
    688694    # dimension definitions - fixed
    689695    outfile.createDimension('num_of_dimensions', 2) #This is 2d data
    … …  
    694700
    695701    # Create dimensions, variables and set the variables
    696    
     702
    697703    # trianglulation
    698704    # vertices
    … …  
    702708                                                   'num_of_dimensions'))
    703709        outfile.variables['vertices'][:] = mesh['vertices']
    704         if mesh['vertex_attributes']  != None and \
    705                (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0):
     710        if mesh['vertex_attributes']  != None \
     711           and (mesh['vertex_attributes'].shape[0] > 0 \
     712           and mesh['vertex_attributes'].shape[1] > 0):
    706713            outfile.createDimension('num_of_vertex_attributes',
    707714                                    mesh['vertex_attributes'].shape[1])
    … …  
    711718                                   Float,
    712719                                   ('num_of_vertices',
    713                                     'num_of_vertex_attributes'))   
     720                                    'num_of_vertex_attributes'))
    714721            outfile.createVariable('vertex_attribute_titles',
    715722                                   Character,
    716                                    ( 'num_of_vertex_attributes',
    717                                      'num_of_vertex_attribute_title_chars' ))
     723                                   ('num_of_vertex_attributes',
     724                                    'num_of_vertex_attribute_title_chars'))
    718725            outfile.variables['vertex_attributes'][:] = \
    719                                                       mesh['vertex_attributes']
     726                                     mesh['vertex_attributes']
    720727            outfile.variables['vertex_attribute_titles'][:] = \
    721728                                     mesh['vertex_attribute_titles']
     729
    722730    # segments
    723     if (mesh['segments'].shape[0] > 0):       
    724         outfile.createDimension('num_of_segments',
    725                                 mesh['segments'].shape[0])
    726         outfile.createVariable('segments',
    727                                IntType,
     731    if (mesh['segments'].shape[0] > 0):
     732        outfile.createDimension('num_of_segments', mesh['segments'].shape[0])
     733        outfile.createVariable('segments', IntType,
    728734                               ('num_of_segments', 'num_of_segment_ends'))
    729735        outfile.variables['segments'][:] = mesh['segments']
    … …  
    736742                                    'num_of_segment_tag_chars'))
    737743            outfile.variables['segment_tags'][:] = mesh['segment_tags']
    738     # triangles   
     744
     745    # triangles
    739746    if (mesh['triangles'].shape[0] > 0):
    740         outfile.createDimension('num_of_triangles',
    741                                 mesh['triangles'].shape[0])
     747        outfile.createDimension('num_of_triangles', mesh['triangles'].shape[0])
    742748        outfile.createVariable('triangles',
    743749                               IntType,
    744                                ('num_of_triangles',
    745                                 'num_of_triangle_vertices'))
     750                               ('num_of_triangles', 'num_of_triangle_vertices'))
    746751        outfile.createVariable('triangle_neighbors',
    747752                               IntType,
    748                                ('num_of_triangles',
    749                                 'num_of_triangle_faces'))
     753                               ('num_of_triangles', 'num_of_triangle_faces'))
    750754        outfile.variables['triangles'][:] = mesh['triangles']
    751755        outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
    752         if mesh['triangle_tags'] != None and \
    753                (mesh['triangle_tags'].shape[1] > 0):
     756        if mesh['triangle_tags'] != None \
     757           and (mesh['triangle_tags'].shape[1] > 0):
    754758            outfile.createDimension('num_of_triangle_tag_chars',
    755                                     mesh['triangle_tags'].shape[1]) 
     759                                    mesh['triangle_tags'].shape[1])
    756760            outfile.createVariable('triangle_tags',
    757761                                   Character,
    … …  
    759763                                    'num_of_triangle_tag_chars'))
    760764            outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
    761    
    762765
    763766    # outline
    … …  
    767770        outfile.createVariable('points', Float, ('num_of_points',
    768771                                                 'num_of_dimensions'))
    769         outfile.variables['points'][:] = mesh['points']
    770         if (mesh['point_attributes'].shape[0] > 0  and mesh['point_attributes'].shape[1] > 0):
     772        outfile.variables['points'][:] = mesh['points']
     773        if mesh['point_attributes'].shape[0] > 0  \
     774           and mesh['point_attributes'].shape[1] > 0:
    771775            outfile.createDimension('num_of_point_attributes',
    772776                                    mesh['point_attributes'].shape[1])
    773777            outfile.createVariable('point_attributes',
    774778                                   Float,
    775                                    ('num_of_points',
    776                                     'num_of_point_attributes'))
     779                                   ('num_of_points', 'num_of_point_attributes'))
    777780            outfile.variables['point_attributes'][:] = mesh['point_attributes']
    778     # outline_segments
    779     if (mesh['outline_segments'].shape[0] > 0):
     781
     782    # outline_segments
     783    if mesh['outline_segments'].shape[0] > 0:
    780784        outfile.createDimension('num_of_outline_segments',
    781785                                mesh['outline_segments'].shape[0])
    … …  
    792796                                   ('num_of_outline_segments',
    793797                                    'num_of_outline_segment_tag_chars'))
    794             outfile.variables['outline_segment_tags'][:] = mesh['outline_segment_tags']
     798            outfile.variables['outline_segment_tags'][:] = \
     799                mesh['outline_segment_tags']
     800
    795801    # holes
    796802    if (mesh['holes'].shape[0] > 0):
    … …  
    799805                                                'num_of_dimensions'))
    800806        outfile.variables['holes'][:] = mesh['holes']
     807
    801808    # regions
    802809    if (mesh['regions'].shape[0] > 0):
    … …  
    804811        outfile.createVariable('regions', Float, ('num_of_regions',
    805812                                                  'num_of_dimensions'))
    806         outfile.createVariable('region_max_areas',
    807                                Float,
    808                                ('num_of_regions',))
     813        outfile.createVariable('region_max_areas', Float, ('num_of_regions',))
    809814        outfile.variables['regions'][:] = mesh['regions']
    810815        outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
    … …  
    812817            outfile.createDimension('num_of_region_tag_chars',
    813818                                    mesh['region_tags'].shape[1])
    814             outfile.createVariable('region_tags',
    815                                    Character,
     819            outfile.createVariable('region_tags', Character,
    816820                                   ('num_of_regions',
    817821                                    'num_of_region_tag_chars'))
    … …  
    821825    if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
    822826        mesh['geo_reference'].write_NetCDF(outfile)
    823                                                  
     827
    824828    outfile.close()
    825829
    826830
    827    
     831##
     832# @brief Read a mesh file.
     833# @param file_name Path to the file to read.
     834# @return A dictionary containing the .msh file data.
    828835def _read_msh_file(file_name):
    829     """
    830     Read in an msh file.
    831 
    832     """
    833                
    834            
     836    """ Read in an msh file.
     837    """
     838
    835839    #Check contents
    836840    #Get NetCDF
    837    
    838841    # see if the file is there.  Throw a QUIET IO error if it isn't
    839     fd = open(file_name,'r')
     842    fd = open(file_name, 'r')
    840843    fd.close()
    841844
    842845    #throws prints to screen if file not present
    843     fid = NetCDFFile(file_name, 'r') 
     846    fid = NetCDFFile(file_name, 'r')
    844847    mesh = {}
     848
    845849    # Get the variables
    846850    # the triangulation
    … …  
    849853    except KeyError:
    850854        mesh['vertices'] = array([])
     855
    851856    try:
    852857        mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
    … …  
    855860        #for ob in mesh['vertices']:
    856861            #mesh['vertex_attributes'].append([])
    857    
     862
    858863    mesh['vertex_attribute_titles'] = []
    859864    try:
    … …  
    862867            mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
    863868    except KeyError:
    864         pass 
     869        pass
     870
    865871    try:
    866872        mesh['segments'] = fid.variables['segments'][:]
    867873    except KeyError:
    868874        mesh['segments'] = array([])
     875
    869876    mesh['segment_tags'] =[]
    870877    try:
    … …  
    875882        for ob in mesh['segments']:
    876883            mesh['segment_tags'].append('')
     884
    877885    try:
    878886        mesh['triangles'] = fid.variables['triangles'][:]
    … …  
    881889        mesh['triangles'] = array([])
    882890        mesh['triangle_neighbors'] = array([])
     891
    883892    mesh['triangle_tags'] =[]
    884893    try:
    … …  
    889898        for ob in mesh['triangles']:
    890899            mesh['triangle_tags'].append('')
    891    
    892     #the outline 
     900
     901    #the outline
    893902    try:
    894903        mesh['points'] = fid.variables['points'][:]
    895904    except KeyError:
    896905        mesh['points'] = []
     906
    897907    try:
    898908        mesh['point_attributes'] = fid.variables['point_attributes'][:]
    … …  
    901911        for point in mesh['points']:
    902912            mesh['point_attributes'].append([])
     913
    903914    try:
    904915        mesh['outline_segments'] = fid.variables['outline_segments'][:]
    905916    except KeyError:
    906917        mesh['outline_segments'] = array([])
     918
    907919    mesh['outline_segment_tags'] =[]
    908920    try:
    … …  
    913925        for ob in mesh['outline_segments']:
    914926            mesh['outline_segment_tags'].append('')
     927
    915928    try:
    916929        mesh['holes'] = fid.variables['holes'][:]
    917930    except KeyError:
    918931        mesh['holes'] = array([])
     932
    919933    try:
    920934        mesh['regions'] = fid.variables['regions'][:]
    921935    except KeyError:
    922936        mesh['regions'] = array([])
     937
    923938    mesh['region_tags'] =[]
    924939    try:
    … …  
    929944        for ob in mesh['regions']:
    930945            mesh['region_tags'].append('')
     946
    931947    try:
    932948        mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
    … …  
    934950        mesh['region_max_areas'] = array([])
    935951    #mesh[''] = fid.variables[''][:]
    936      
     952
    937953    try:
    938954        geo_reference = Geo_reference(NetCDFObject=fid)
    939955        mesh['geo_reference'] = geo_reference
    940956    except AttributeError, e:
    941         #geo_ref not compulsory 
     957        #geo_ref not compulsory
    942958        mesh['geo_reference'] = None
    943    
     959
    944960    fid.close()
    945      
     961
    946962    return mesh
    947            
    948 
    949 ## used by alpha shapes
    950    
    951 def export_boundary_file( file_name, points, title, delimiter = ','):
     963
     964
     965##
     966# @brief Export a boundary file.
     967# @param file_name Path to the file to write.
     968# @param points List of point index pairs.
     969# @paran title Title string to write into file (first line).
     970# @param delimiter Delimiter string to write between points.
     971# @note Used by alpha shapes
     972def export_boundary_file(file_name, points, title, delimiter=','):
    952973    """
    953974    export a file, ofile, with the format
    954    
     975
    955976    First line: Title variable
    956     Following lines:  [point index][delimiter][point index] 
     977    Following lines:  [point index][delimiter][point index]
    957978
    958979    file_name - the name of the new file
    959     points - List of point index pairs [[p1, p2],[p3, p4]..] 
     980    points - List of point index pairs [[p1, p2],[p3, p4]..]
    960981    title - info to write in the first line
    961982    """
    962    
    963     fd = open(file_name,'w')
    964    
    965     fd.write(title+"\n")
    966     #[point index][delimiter][point index]
     983
     984    fd = open(file_name, 'w')
     985
     986    fd.write(title + "\n")
     987
     988    #[point index][delimiter][point index]
    967989    for point in points:
    968         fd.write( str(point[0]) + delimiter
    969                   + str(point[1]) + "\n")
     990        fd.write(str(point[0]) + delimiter + str(point[1]) + "\n")
     991
    970992    fd.close()
    971993
    … …  
    975997###
    976998
     999##
     1000# @brief
     1001# @param point_atts
    9771002def extent_point_atts(point_atts):
    9781003    """
    … …  
    9821007    point_atts['pointlist'] = extent(point_atts['pointlist'])
    9831008    point_atts['attributelist'] = {}
     1009
    9841010    return point_atts
    985    
     1011
     1012
     1013##
     1014# @brief Get an extent array for a set of points.
     1015# @param points A set of points.
     1016# @return An extent array of form [[min_x, min_y], [max_x, min_y],
     1017#                                  [max_x, max_y], [min_x, max_y]]
    9861018def extent(points):
    9871019    points = array(points).astype(Float)
     1020
    9881021    max_x = min_x = points[0][0]
    9891022    max_y = min_y = points[0][1]
     1023
    9901024    for point in points[1:]:
    991         x = point[0] 
     1025        x = point[0]
    9921026        if x > max_x: max_x = x
    993         if x < min_x: min_x = x           
    994         y = point[1]
     1027        if x < min_x: min_x = x
     1028
     1029        y = point[1]
    9951030        if y > max_y: max_y = y
    9961031        if y < min_y: min_y = y
    997     extent = array([[min_x,min_y],[max_x,min_y],[max_x,max_y],[min_x,max_y]])
    998     #print "extent",extent
     1032
     1033    extent = array([[min_x, min_y],
     1034                    [max_x, min_y],
     1035                    [max_x, max_y],
     1036                    [min_x, max_y]])
     1037
    9991038    return extent
    1000    
     1039
     1040
     1041##
     1042# @brief Reduce a points file until number of points is less than limit.
     1043# @param infile Path to input file to thin.
     1044# @param outfile Path to output thinned file.
     1045# @param max_points Max number of points in output file.
     1046# @param verbose True if this function is to be verbose.
    10011047def reduce_pts(infile, outfile, max_points, verbose = False):
    10021048    """
    1003     reduces a points file by removong every second point until the # of points
     1049    reduces a points file by removing every second point until the # of points
    10041050    is less than max_points.
    10051051    """
     1052
    10061053    # check out pts2rectangular in least squares, and the use of reduction.
    10071054    # Maybe it does the same sort of thing?
    10081055    point_atts = _read_pts_file(infile)
     1056
    10091057    while point_atts['pointlist'].shape[0] > max_points:
    10101058        if verbose: print "point_atts['pointlist'].shape[0]"
    10111059        point_atts = half_pts(point_atts)
     1060
    10121061    export_points_file(outfile, point_atts)
    1013  
    1014 def produce_half_point_files(infile, max_points, delimiter, verbose = False):
     1062
     1063
     1064##
     1065# @brief
     1066# @param infile
     1067# @param max_points
     1068# @param delimiter
     1069# @param verbose True if this function is to be verbose.
     1070# @return A list of
     1071def produce_half_point_files(infile, max_points, delimiter, verbose=False):
    10151072    point_atts = _read_pts_file(infile)
    10161073    root, ext = splitext(infile)
    10171074    outfiles = []
     1075
    10181076    if verbose: print "# of points", point_atts['pointlist'].shape[0]
     1077
    10191078    while point_atts['pointlist'].shape[0] > max_points:
    10201079        point_atts = half_pts(point_atts)
     1080
    10211081        if verbose: print "# of points", point_atts['pointlist'].shape[0]
     1082
    10221083        outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
    10231084        outfiles.append(outfile)
    1024         export_points_file(outfile,
    1025                   point_atts)
     1085        export_points_file(outfile, point_atts)
     1086
    10261087    return outfiles
    1027    
     1088
     1089
     1090##
     1091# @brief
     1092# @param point_atts Object with attribute of 'points list'.
     1093# @return
    10281094def point_atts2array(point_atts):
     1095    # convert attribute list to array of floats
    10291096    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
    1030    
     1097
    10311098    for key in point_atts['attributelist'].keys():
    1032         point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
     1099        point_atts['attributelist'][key] = \
     1100            array(point_atts['attributelist'][key]).astype(Float)
     1101
    10331102    return point_atts
    1034    
     1103
     1104
     1105##
     1106# @brief ??
     1107# @param point_atts ??
     1108# @return ??
    10351109def half_pts(point_atts):
    10361110    point_atts2array(point_atts)
    10371111    point_atts['pointlist'] = point_atts['pointlist'][::2]
    1038    
     1112
    10391113    for key in point_atts['attributelist'].keys():
    1040         point_atts['attributelist'][key]=  point_atts['attributelist'][key] [::2]
     1114        point_atts['attributelist'][key] = point_atts['attributelist'][key][::2]
     1115
    10411116    return point_atts
    10421117
     1118##
     1119# @brief ??
     1120# @param dic ??
     1121# @return ??
    10431122def concatinate_attributelist(dic):
    10441123    """
    10451124    giving a dic[attribute title] = attribute
    1046     return list of attribute titles, array of attributes
    1047     """
     1125    return list of attribute titles, array of attributes
     1126    """
     1127
    10481128    point_attributes = array([]).astype(Float)
    10491129    keys = dic.keys()
    … …  
    10541134        reshaped = reshape(dic[key],(dic[key].shape[0],1))
    10551135        point_attributes = concatenate([point_attributes, reshaped], axis=1)
     1136
    10561137    return dic.keys(), point_attributes
    10571138
     1139
     1140##
     1141# @brief
     1142# @param dict
     1143# @param indices_to_keep
     1144# @return
    10581145#FIXME(dsg), turn this dict plus methods into a class?
    1059 def take_points(dict,indices_to_keep):
     1146def take_points(dict, indices_to_keep):
    10601147    dict = point_atts2array(dict)
    10611148    #FIXME maybe the points data structure should become a class?
    … …  
    10651152        dict['attributelist'][key]= take(dict['attributelist'][key],
    10661153                                         indices_to_keep)
     1154
    10671155    return dict
    1068    
     1156
     1157##
     1158# @brief ??
     1159# @param dict1 ??
     1160# @param dict2 ??
     1161# @return ??
    10691162def add_point_dictionaries (dict1, dict2):
    10701163    """
    10711164    """
     1165
    10721166    dict1 = point_atts2array(dict1)
    10731167    dict2 = point_atts2array(dict2)
    1074    
    1075     combined = {} 
     1168
     1169    combined = {}
    10761170    combined['pointlist'] = concatenate((dict2['pointlist'],
    1077                                          dict1['pointlist']),axis=0)   
    1078     atts = {}   
     1171                                         dict1['pointlist']),axis=0)
     1172
     1173    atts = {}
    10791174    for key in dict2['attributelist'].keys():
    10801175        atts[key]= concatenate((dict2['attributelist'][key],
    … …  
    10821177    combined['attributelist']=atts
    10831178    combined['geo_reference'] = dict1['geo_reference']
     1179
    10841180    return combined
    1085      
     1181
     1182
    10861183if __name__ == "__main__":
    10871184    pass

  • anuga_core/source/anuga/load_mesh/test_loadASCII.py

    r5875 r6074  
    543543                        'bad msh file did not raise error!')       
    544544        os.remove(fileName)         
     545
     546######
     547# Test the clean_line() utility function.
     548######
     549
     550    def test_clean_line_01(self):
     551        test_line = 'abc, ,,xyz,123'
     552        result = clean_line(test_line, ',')
     553        self.failUnless(result == ['abc', '', 'xyz', '123'])
     554             
     555    def test_clean_line_02(self):
     556        test_line = ' abc , ,, xyz  , 123  '
     557        result = clean_line(test_line, ',')
     558        self.failUnless(result == ['abc', '', 'xyz', '123'])
     559             
     560    def test_clean_line_03(self):
     561        test_line = '1||||2'
     562        result = clean_line(test_line, '|')
     563        self.failUnless(result == ['1', '2'])
    545564             
    546565#-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.