Python utilities.verify函数代码示例

def rings_equal(x, y, rtol=1.0e-6, atol=1.0e-8):
    """Compares to linear rings as numpy arrays

    Args
        * x, y: Nx2 numpy arrays

    Returns:
        * True if x == y or x' == y (up to the specified tolerance)

        where x' is x reversed in the first dimension. This corresponds to
        linear rings being seen as equal irrespective of whether they are
        organised in clock wise or counter clock wise order
    """

    x = ensure_numeric(x, numpy.float)
    y = ensure_numeric(y, numpy.float)

    msg = 'Arrays must a 2d arrays of vertices. I got %s and %s' % (x, y)
    verify(len(x.shape) == 2 and len(y.shape) == 2, msg)

    msg = 'Arrays must have two columns. I got %s and %s' % (x, y)
    verify(x.shape[1] == 2 and y.shape[1] == 2, msg)

    if (numpy.allclose(x, y, rtol=rtol, atol=atol) or
        numpy.allclose(x, y[::-1], rtol=rtol, atol=atol)):
        return True
    else:
        return False

def bboxstring2list(bbox_string):
    """Convert bounding box string to list

    Args:
        * bbox_string: String of bounding box coordinates of the form 'W,S,E,N'

    Returns:
        * bbox: List of floating point numbers with format [W, S, E, N]
    """

    msg = ('Bounding box must be a string with coordinates following the '
           'format 105.592,-7.809,110.159,-5.647\n'
           'Instead I got %s of type %s.' % (str(bbox_string),
                                             type(bbox_string)))
    verify(isinstance(bbox_string, basestring), msg)

    fields = bbox_string.split(',')
    msg = ('Bounding box string must have 4 coordinates in the form '
           '"W,S,E,N". I got bbox == "%s"' % bbox_string)
    try:
        verify(len(fields) == 4, msg)
    except VerificationError:
        raise BoundingBoxError(msg)

    for x in fields:
        try:
            float(x)
        except ValueError, e:
            msg = ('Bounding box %s contained non-numeric entry %s, '
                   'original error was "%s".' % (bbox_string, x, e))
            raise BoundingBoxError(msg)

def geotransform2resolution(geotransform, isotropic=False, rtol=1.0e-6, atol=1.0e-8):
    """Convert geotransform to resolution

    Args:
        * geotransform: GDAL geotransform (6-tuple).
                        (top left x, w-e pixel resolution, rotation,
                        top left y, rotation, n-s pixel resolution).
                        See e.g. http://www.gdal.org/gdal_tutorial.html
        * isotropic: If True, verify that dx == dy and return dx
                     If False (default) return 2-tuple (dx, dy)
        * rtol, atol: Used to control how close dx and dy must be
                      to quality for isotropic. These are passed on to
                      numpy.allclose for comparison.

    Returns:
        * resolution: grid spacing (resx, resy) in (positive) decimal
                      degrees ordered as longitude first, then latitude.
                      or resx (if isotropic is True)
    """

    resx = geotransform[1]  # w-e pixel resolution
    resy = -geotransform[5]  # n-s pixel resolution (always negative)

    if isotropic:
        msg = (
            "Resolution requested with "
            "isotropic=True, but "
            "resolutions in the horizontal and vertical "
            "are different: resx = %.12f, resy = %.12f. " % (resx, resy)
        )
        verify(numpy.allclose(resx, resy, rtol=rtol, atol=atol), msg)

        return resx
    else:
        return resx, resy

def bboxlist2string(bbox, decimals=6):
    """Convert bounding box list to comma separated string

    Args:
        * bbox: List of coordinates of the form [W, S, E, N]

    Returns:
        * bbox_string: Format 'W,S,E,N' - each will have 6 decimal points
    """

    msg = 'Got string %s, but expected bounding box as a list' % str(bbox)
    verify(not isinstance(bbox, basestring), msg)

    try:
        bbox = list(bbox)
    except:
        msg = 'Could not coerce bbox %s into a list' % str(bbox)
        raise BoundingBoxError(msg)

    msg = ('Bounding box must have 4 coordinates [W, S, E, N]. '
           'I got %s' % str(bbox))
    try:
        verify(len(bbox) == 4, msg)
    except VerificationError:
        raise BoundingBoxError(msg)

    for x in bbox:
        try:
            float(x)
        except ValueError, e:
            msg = ('Bounding box %s contained non-numeric entry %s, '
                   'original error was "%s".' % (bbox, x, e))
            raise BoundingBoxError(msg)

def convert_polygons_to_centroids(V):
    """Convert polygon vector data to point vector data

    Args:
        * V: Vector layer with polygon data

    Returns:
        * Vector layer with point data and the same attributes as V
    """

    msg = 'Input data %s must be polygon vector data' % V
    verify(V.is_polygon_data, msg)

    geometry = V.get_geometry()
    N = len(V)

    # Calculate points for each polygon
    centroids = []
    for i in range(N):
        c = calculate_polygon_centroid(geometry[i])
        centroids.append(c)

    # Create new point vector layer with same attributes and return
    V = Vector(data=V.get_data(),
               projection=V.get_projection(),
               geometry=centroids,
               name='%s_centroid_data' % V.get_name(),
               keywords=V.get_keywords())
    return V

def _keywords_to_string(keywords, sublayer=None):
    """Create a string from a keywords dict.

    Args:
        * keywords: A required dictionary containing the keywords to stringify.
        * sublayer: str optional group marker for a sub layer.

    Returns:
        str: a String containing the rendered keywords list

    Raises:
        Any exceptions are propogated.

    .. note: Only simple keyword dicts should be passed here, not multilayer
       dicts.

    For example you pass a dict like this::

        {'datatype': 'osm',
         'category': 'exposure',
         'title': 'buildings_osm_4326',
         'subcategory': 'building',
         'purpose': 'dki'}

    and the following string would be returned:

        datatype: osm
        category: exposure
        title: buildings_osm_4326
        subcategory: building
        purpose: dki

    If sublayer is provided e.g. _keywords_to_string(keywords, sublayer='foo'),
    the following:

        [foo]
        datatype: osm
        category: exposure
        title: buildings_osm_4326
        subcategory: building
        purpose: dki
    """

    # Write
    result = get_unicode("")
    if sublayer is not None:
        result = "[%s]\n" % sublayer
    for k, value in keywords.items():
        # Create key
        msg = "Key in keywords dictionary must be a string. " "I got %s with type %s" % (k, str(type(k))[1:-1])
        verify(isinstance(k, basestring), msg)

        key = k
        msg = 'Key in keywords dictionary must not contain the ":" ' 'character. I got "%s"' % key
        verify(":" not in key, msg)

        # Store
        result += "%s: %s\n" % (key, value)
    return result

def bbox_intersection(*args):
    """Compute intersection between two or more bounding boxes

    Args:
        * args: two or more bounding boxes.
              Each is assumed to be a list or a tuple with
              four coordinates (W, S, E, N)

    Returns:
        * result: The minimal common bounding box

    """

    msg = 'Function bbox_intersection must take at least 2 arguments.'
    verify(len(args) > 1, msg)

    result = [-180, -90, 180, 90]
    for a in args:
        if a is None:
            continue

        msg = ('Bounding box expected to be a list of the '
               'form [W, S, E, N]. '
               'Instead i got "%s"' % str(a))

        try:
            box = list(a)
        except:
            raise Exception(msg)

        if not len(box) == 4:
            raise BoundingBoxError(msg)

        msg = ('Western boundary must be less than or equal to eastern. '
               'I got %s' % box)
        if not box[0] <= box[2]:
            raise BoundingBoxError(msg)

        msg = ('Southern boundary must be less than or equal to northern. '
               'I got %s' % box)
        if not box[1] <= box[3]:
            raise BoundingBoxError(msg)

        # Compute intersection

        # West and South
        for i in [0, 1]:
            result[i] = max(result[i], box[i])

        # East and North
        for i in [2, 3]:
            result[i] = min(result[i], box[i])

    # Check validity and return
    if result[0] <= result[2] and result[1] <= result[3]:
        return result
    else:
        return None

def interpolate_raster_vector(source,
                              target,
                              layer_name=None,
                              attribute_name=None,
                              mode='linear'):
    """Interpolate from raster layer to vector data

    Args:
        * source: Raster data set (grid)
        * target: Vector data set (points or polygons)
        * layer_name: Optional name of returned interpolated layer.
              If None the name of V is used for the returned layer.
        * attribute_name: Name for new attribute.
              If None (default) the name of R is used

    Returns:
        I: Vector data set; points located as target with values
           interpolated from source

    Note: If target geometry is polygon, data will be interpolated to
    its centroids and the output is a point data set.
    """

    # Input checks
    verify(source.is_raster)
    verify(target.is_vector)

    if target.is_point_data:
        # Interpolate from raster to point data
        R = interpolate_raster_vector_points(
            source,
            target,
            layer_name=layer_name,
            attribute_name=attribute_name,
            mode=mode)
    #elif target.is_line_data:
    # TBA - issue https://github.com/AIFDR/inasafe/issues/36
    #
    elif target.is_polygon_data:
        # Use centroids, in case of polygons
        P = convert_polygons_to_centroids(target)
        R = interpolate_raster_vector_points(
            source, P, layer_name=layer_name, attribute_name=attribute_name)
        # In case of polygon data, restore the polygon geometry
        # Do this setting the geometry of the returned set to
        # that of the original polygon
        R = Vector(
            data=R.get_data(),
            projection=R.get_projection(),
            geometry=target.get_geometry(),
            name=R.get_name())
    else:
        msg = ('Unknown datatype for raster2vector interpolation: '
               'I got %s' % str(target))
        raise InaSAFEError(msg)

    # Return interpolated vector layer
    return R

    def copy_keywords(
            self,
            source_layer,
            destination_file,
            extra_keywords=None):
        """Helper to copy the keywords file from a source to a target dataset.

        e.g.::

            copyKeywords('foo.shp', 'bar.shp')

        Will result in the foo.keywords file being copied to bar.keyword.

        Optional argument extraKeywords is a dictionary with additional
        keywords that will be added to the destination file e.g::

            copyKeywords('foo.shp', 'bar.shp', {'resolution': 0.01})

        :param source_layer: A QGIS QgsMapLayer instance.
        :type source_layer: QgsMapLayer

        :param destination_file: The output filename that should be used
            to store the keywords in. It can be a .shp or a .keywords for
            example since the suffix will always be replaced with .keywords.
        :type destination_file: str

        :param extra_keywords: A dict containing all the extra keywords
            to be written for the layer. The written keywords will consist of
            any original keywords from the source layer's keywords file and
            and the extra keywords (which will replace the source layers
            keywords if the key is identical).
        :type extra_keywords: dict

        """
        keywords = self.read_keywords(source_layer)
        if extra_keywords is None:
            extra_keywords = {}
        message = self.tr(
            'Expected extraKeywords to be a dictionary. Got '
            '%s' % str(type(extra_keywords))[1:-1])
        verify(isinstance(extra_keywords, dict), message)
        # compute the output keywords file name
        destination_base = os.path.splitext(destination_file)[0]
        new_destination = destination_base + '.keywords'
        # write the extra keywords into the source dict
        try:
            for key in extra_keywords:
                keywords[key] = extra_keywords[key]
            write_keywords_to_file(new_destination, keywords)
        except Exception, e:
            message = self.tr(
                'Failed to copy keywords file from : \n%s\nto\n%s: %s' % (
                    source_layer.source(), new_destination, str(e)))
            raise Exception(message)

def tag_polygons_by_grid(polygons, grid, threshold=0, tag='affected'):
    """Tag polygons by raster values

    Args:
        * polygons: Polygon layer
        * grid: Raster layer
        * threshold: Threshold for grid value to tag polygon
        * tag: Name of new tag

    Returns:
        Polygon layer: Same as input polygon but with extra attribute tag
                       set according to grid values

    """

    verify(polygons.is_polygon_data)
    verify(grid.is_raster)

    polygon_attributes = polygons.get_data()
    polygon_geometry = polygons.get_geometry(as_geometry_objects=True)

    # Separate grid points by polygon
    res, _ = clip_grid_by_polygons(
        grid.get_data(),
        grid.get_geotransform(),
        polygon_geometry)

    # Create new polygon layer with tag set according to grid values
    # and threshold
    new_attributes = []
    for i, (_, values) in enumerate(res):
        # For each polygon check if any grid value in it exceeds the threshold
        affected = False
        for val in values:
            # Check each grid value in this polygon
            if val > threshold:
                affected = True

        # Existing attributes for this polygon
        attr = polygon_attributes[i].copy()

        # Create tagged polygon feature
        if affected:
            attr[tag] = True
        else:
            attr[tag] = False

        new_attributes.append(attr)

    R = Vector(data=new_attributes,
               projection=polygons.get_projection(),
               geometry=polygon_geometry,
               name='%s_tagged_by_%s' % (polygons.name, grid.name))
    return R

def dom2object(node):
    """Convert DOM representation to XML_object hierarchy.
    """

    value = []
    textnode_encountered = None
    for n in node.childNodes:

        if n.nodeType == 3:
            # Child is a text element - omit the dom tag #text and
            # go straight to the text value.

            # Note - only the last text value will be recorded

            msg = 'Text element has child nodes - this shouldn\'t happen'
            verify(len(n.childNodes) == 0, msg)


            x = n.nodeValue.strip()
            if len(x) == 0:
                # Skip empty text children
                continue

            textnode_encountered = value = x
        else:
            # XML element


            if textnode_encountered is not None:
                msg = 'A text node was followed by a non-text tag. This is not allowed.\n'
                msg += 'Offending text node: "%s" ' %str(textnode_encountered)
                msg += 'was followed by node named: "<%s>"' %str(n.nodeName)
                raise Exception, msg


            value.append(dom2object(n))


    # Deal with empty elements
    if len(value) == 0: value = ''


    if node.nodeType == 9:
        # Root node (document)
        tag = None
    else:
        # Normal XML node
        tag = node.nodeName


    X = XML_element(tag=tag,
                    value=value)

    return X

def sigab2bnpb(E, target_attribute='VCLASS'):
    """Map SIGAB point data to BNPB vulnerability classes

    Input
        E: Vector object representing the OSM data
        target_attribute: Optional name of the attribute containing
                          the mapped vulnerability class. Default
                          value is 'VCLASS'

    Output:
        Vector object like E, but with one new attribute (e.g. 'VCLASS')
        representing the vulnerability class used in the guidelines
    """

    # Input check
    required = ['Struktur_B', 'Lantai', 'Atap', 'Dinding', 'Tingkat']
    actual = E.get_attribute_names()

    msg = ('Input data to sigab2bnpb must have attributes %s. '
           'It has %s' % (str(required), str(actual)))
    for attribute in required:
        verify(attribute in actual, msg)

    # Start mapping
    N = len(E)
    attributes = E.get_data()
    for i in range(N):
        levels = E.get_data('Tingkat', i).lower()
        structure = E.get_data('Struktur_B', i).lower()
        roof_type = E.get_data('Atap', i).lower()
        wall_type = E.get_data('Dinding', i).lower()
        floor_type = E.get_data('Lantai', i).lower()
        if levels == 'none' or structure == 'none':
            vulnerability_class = 'URM'
        elif structure.startswith('beton') or structure.startswith('kayu'):
            vulnerability_class = 'RM'
        else:
            if int(levels) >= 2:
                vulnerability_class = 'RM'
            else:
                vulnerability_class = 'URM'

        # Store new attribute value
        attributes[i][target_attribute] = vulnerability_class

    # Create new vector instance and return
    V = Vector(data=attributes,
               projection=E.get_projection(),
               geometry=E.get_geometry(),
               name=E.get_name() + ' mapped to BNPB vulnerability classes',
               keywords=E.get_keywords())
    return V

def calculate_polygon_area(polygon, signed=False):
    """Calculate the signed area of non-self-intersecting polygon

    :param polygon: Numeric array of points (longitude, latitude). It is
        assumed to be closed, i.e. first and last points are identical
    :type polygon: numpy.ndarray

    :param signed: Optional flag deciding whether returned area retains its
     sign:

            If points are ordered counter clockwise, the signed area
            will be positive.

            If points are ordered clockwise, it will be negative
            Default is False which means that the area is always
            positive.

    :type signed: bool

    :returns: area: Area of polygon (subject to the value of argument signed)
    :rtype: numpy.ndarray

    Note:
        Sources
            http://paulbourke.net/geometry/polyarea/
            http://en.wikipedia.org/wiki/Centroid
    """

    # Make sure it is numeric
    P = numpy.array(polygon)

    msg = ('Polygon is assumed to consist of coordinate pairs. '
           'I got second dimension %i instead of 2' % P.shape[1])
    verify(P.shape[1] == 2, msg)

    x = P[:, 0]
    y = P[:, 1]

    # Calculate 0.5 sum_{i=0}^{N-1} (x_i y_{i+1} - x_{i+1} y_i)
    a = x[:-1] * y[1:]
    b = y[:-1] * x[1:]

    A = numpy.sum(a - b) / 2.

    if signed:
        return A
    else:
        return abs(A)

def raster_geometry_to_geotransform(longitudes, latitudes):
    """Convert vectors of longitudes and latitudes to geotransform

    Note:
        This is the inverse operation of Raster.get_geometry().

    :param longitudes: Vectors of geographic coordinates
    :type longitudes:

    :param latitudes: Vectors of geographic coordinates
    :type latitudes:

    :returns: geotransform: 6-tuple (top left x, w-e pixel resolution,
        rotation, top left y, rotation, n-s pixel resolution)
    """

    nx = len(longitudes)
    ny = len(latitudes)

    msg = ('You must specify more than 1 longitude to make geotransform: '
           'I got %s' % str(longitudes))
    verify(nx > 1, msg)

    msg = ('You must specify more than 1 latitude to make geotransform: '
           'I got %s' % str(latitudes))
    verify(ny > 1, msg)

    dx = float(longitudes[1] - longitudes[0])  # Longitudinal resolution
    dy = float(latitudes[0] - latitudes[1])  # Latitudinal resolution (neg)

    # Define pixel centers along each directions
    # This is to achieve pixel registration rather
    # than gridline registration
    dx2 = dx / 2
    dy2 = dy / 2

    geotransform = (longitudes[0] - dx2,  # Longitude of upper left corner
                    dx,                   # w-e pixel resolution
                    0,                    # rotation
                    latitudes[-1] - dy2,  # Latitude of upper left corner
                    0,                    # rotation
                    dy)                   # n-s pixel resolution

    return geotransform

    def write_to_file(self, filename):
        """Save raster data to file

        Args:
            * filename: filename with extension .tif
        """

        # Check file format
        basename, extension = os.path.splitext(filename)

        msg = ('Invalid file type for file %s. Only extension '
               'tif allowed.' % filename)
        verify(extension in ['.tif', '.asc'], msg)
        file_format = DRIVER_MAP[extension]

        # Get raster data
        A = self.get_data()

        # Get Dimensions. Note numpy and Gdal swap order
        N, M = A.shape

        # Create empty file.
        # FIXME (Ole): It appears that this is created as single
        #              precision even though Float64 is specified
        #              - see issue #17
        driver = gdal.GetDriverByName(file_format)
        fid = driver.Create(filename, M, N, 1, gdal.GDT_Float64)
        if fid is None:
            msg = ('Gdal could not create filename %s using '
                   'format %s' % (filename, file_format))
            raise WriteLayerError(msg)

        self.filename = filename

        # Write metada
        fid.SetProjection(str(self.projection))
        fid.SetGeoTransform(self.geotransform)

        # Write data
        fid.GetRasterBand(1).WriteArray(A)

        # Write keywords if any
        write_keywords(self.keywords, basename + '.keywords')