source: issm/trunk-jpl/src/m/shp/shpread.py@ 25065

from os import path
try:
    import shapefile
except ImportError:
    print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled")

from pairoptions import pairoptions


def shpread(filename, *args): #{{{
    '''
    SHPREAD - read a shape file and build a dict

    This routine reads a shape file .shp and builds a dict array containing the 
    fields x and y corresponding to the coordinates, one for the filename 
    of the shp file, for the density, for the nodes, and a field closed to 
    indicate if the domain is closed. If this initial shapefile is point only, 
    the fields closed and points are ommited.
    The first argument is the .shp file to be read and the second one
    (optional) indicates if the last point shall be read (1 to read it, 0 not 
    to).

    The current implementation of shpread depends on PyShp.

    Usage:
        dict = shpread(filename)

    Example:
        From underling PyShp implementation, "The shapefile format is actually 
        a collection of three files. You specify the base filename of the 
        shapefile or the complete filename of any of the shapefile component 
        files.""

        dict = shpread('domainoutline.shp')
        OR
        dict = shpread('domainoutline.dbf')
        OR
        dict = shpread('domainoutline')

        "OR any of the other 5+ formats which are potentially part of a 
        shapefile. The library does not care about file extensions". We do, 
        however, check that a file with the base filename or base filename with 
        .shp extension exists.

    See also EXPDOC, EXPWRITEASVERTICES

    Sources:
    - https://github.com/GeospatialPython/pyshp
    '''

    #recover options
    options = pairoptions(*args)

    #some checks
    if not (path.isfile(filename) or path.isfile(filename + '.shp')):
        raise RuntimeError('shpread error message: file {} or {}.shp not found!'.format(filename, filename))

    #read shapefile
    sf = shapefile.Reader(filename)

    Dicts = []

    #retrieve ID (if it exists)
    fields = sf.fields
    for i in range(1, len(fields)): # skip over first field, which is "DeletionFlag"
        if fields[i][0] == 'id':
            name = str(sf.record(i - 1)[0]) # record index is offset by one, again, because of "DeletionFlag"
            break

    shapes = sf.shapes()
    for i in range(len(shapes)):
        Dict = {}
        shape = shapes[i]
        if shape.shapeType == shapefile.POINT:
            Dict['x'] = shape.points[0][0]
            Dict['y'] = shape.points[0][1]
            Dict['density'] = 1
            Dict['Geometry'] = 'Point'
        elif shape.shapeType == shapefile.POLYLINE:
            num_points = len(shape.points)
            x = []
            y = []
            for j in range(num_points):
                point = shape.points[j]
                x.append(point[0])
                y.append(point[1])
            Dict['x'] = x
            Dict['y'] = y
            Dict['nods'] = num_points
            Dict['density'] = 1
            Dict['closed'] = 1
            Dict['BoundingBox'] = shape.bbox
            Dict['Geometry'] = 'Line'
        elif shape.shapeType == shapefile.POLYGON:
            num_points = len(shape.points)
            x = []
            y = []
            for j in range(num_points):
                point = shape.points[j]
                x.append(point[0])
                y.append(point[1])
            Dict['x'] = x
            Dict['y'] = y
            Dict['nods'] = num_points
            Dict['density'] = 1
            Dict['closed'] = 1
            Dict['BoundingBox'] = shape.bbox
            Dict['Geometry'] =  'Polygon'
        else:
            # NOTE: We could do this once before looping over shapes as all 
            #       shapes in the file must be of the same type, but we would 
            #       need to have a second check anyway in order to know how to 
            #       parse the points. So, let's just assume the file is not 
            #       malformed.
            #
            raise Exception('shpread error: geometry {} is not currently supported'.format(shape.shapeTypeName))
        name = ''
        fields = sf.fields
        for j in range(1, len(fields)): # skip over first field, which is "DeletionFlag"
            fieldname = fields[j][0]
            # 'id' field gets special treatment
            if fieldname == 'id':
                name = str(sf.record(j - 1)[0]) # record index is offset by one, again, because of "DeletionFlag"
            else:
                Dict[fieldname] = sf.record(j - 1)[0]
        Dict['name'] = name
        Dicts.append(Dict)

    invert = options.getfieldvalue('invert', 0)
    if invert:
        for i in range(len(Dicts)):
            Dicts[i].x = np.flipud(Dicts[i].x)
            Dicts[i].y = np.flipud(Dicts[i].y)

    return Dicts
#}}}