Python KeywordIO.copyKeywords方法代码示例

# 需要导入模块: from safe_qgis.keyword_io import KeywordIO [as 别名]
# 或者: from safe_qgis.keyword_io.KeywordIO import copyKeywords [as 别名]

#.........这里部分代码省略.........
        Path to the output clipped layer (placed in the
        system temp dir).

    Raises:
       Exception if input layer is a density layer in projected coordinates -
       see issue #123

    """
    if not theLayer or not theExtent:
        myMessage = tr('Layer or Extent passed to clip is None.')
        raise InvalidParameterException(myMessage)

    if theLayer.type() != QgsMapLayer.RasterLayer:
        myMessage = tr('Expected a raster layer but received a %s.' %
               str(theLayer.type()))
        raise InvalidParameterException(myMessage)

    myWorkingLayer = str(theLayer.source())

    # Check for existence of keywords file
    myKeywordsPath = myWorkingLayer[:-4] + '.keywords'
    myMessage = tr('Input file to be clipped "%s" does not have the '
           'expected keywords file %s' % (myWorkingLayer,
                                          myKeywordsPath))
    verify(os.path.isfile(myKeywordsPath), myMessage)

    # Raise exception if layer is projected and refers to density (issue #123)
    # FIXME (Ole): Need to deal with it - e.g. by automatically reprojecting
    # the layer at this point and setting the native resolution accordingly
    # in its keywords.
    myKeywords = readKeywordsFromFile(myKeywordsPath)
    if 'datatype' in myKeywords and myKeywords['datatype'] == 'density':
        if theLayer.srs().epsg() != 4326:

            # This layer is not WGS84 geographic
            myMessage = ('Layer %s represents density but has spatial '
                         'reference "%s". Density layers must be given in '
                         'WGS84 geographic coordinates, so please reproject '
                         'and try again. For more information, see issue '
                         'https://github.com/AIFDR/inasafe/issues/123'
                         % (myWorkingLayer, theLayer.srs().toProj4()))
            raise InvalidProjectionException(myMessage)

    # We need to provide gdalwarp with a dataset for the clip
    # because unline gdal_translate, it does not take projwin.
    myClipKml = extentToKml(theExtent)

    # Create a filename for the clipped, resampled and reprojected layer
    myHandle, myFilename = tempfile.mkstemp('.tif', 'clip_',
                                            temp_dir())
    os.close(myHandle)
    os.remove(myFilename)

    # If no cell size is specified, we need to run gdalwarp without
    # specifying the output pixel size to ensure the raster dims
    # remain consistent.
    if theCellSize is None:
        myCommand = ('gdalwarp -q -t_srs EPSG:4326 -r near '
                     '-cutline %s -crop_to_cutline -of GTiff '
                     '"%s" "%s"' % (myClipKml,
                                    myWorkingLayer,
                                    myFilename))
    else:
        myCommand = ('gdalwarp -q -t_srs EPSG:4326 -r near -tr %f %f '
                     '-cutline %s -crop_to_cutline -of GTiff '
                     '"%s" "%s"' % (theCellSize,
                                    theCellSize,
                                    myClipKml,
                                    myWorkingLayer,
                                    myFilename))
    myExecutablePrefix = ''
    if sys.platform == 'darwin':  # Mac OS X
        # .. todo:: FIXME - softcode gdal version in this path
        myExecutablePrefix = ('/Library/Frameworks/GDAL.framework/'
                              'Versions/1.9/Programs/')
    myCommand = myExecutablePrefix + myCommand

    LOGGER.debug(myCommand)
    myResult = QProcess().execute(myCommand)

    # For QProcess exit codes see
    # http://qt-project.org/doc/qt-4.8/qprocess.html#execute
    if myResult == -2:  # cannot be started
        myMessageDetail = tr('Process could not be started.')
        myMessage = tr('<p>Error while executing the following shell command:'
                     '</p><pre>%s</pre><p>Error message: %s'
                     % (myCommand, myMessageDetail))
        raise CallGDALError(myMessage)
    elif myResult == -1:  # process crashed
        myMessageDetail = tr('Process could not be started.')
        myMessage = tr('<p>Error while executing the following shell command:'
                       '</p><pre>%s</pre><p>Error message: %s'
                       % (myCommand, myMessageDetail))
        raise CallGDALError(myMessage)

    # .. todo:: Check the result of the shell call is ok
    myKeywordIO = KeywordIO()
    myKeywordIO.copyKeywords(theLayer, myFilename,
                             theExtraKeywords=theExtraKeywords)
    return myFilename  # Filename of created file

# 需要导入模块: from safe_qgis.keyword_io import KeywordIO [as 别名]
# 或者: from safe_qgis.keyword_io.KeywordIO import copyKeywords [as 别名]

#.........这里部分代码省略.........
    if type(theExtent) is list:
        myRect = QgsRectangle(
            theExtent[0], theExtent[1],
            theExtent[2], theExtent[3])
        # noinspection PyCallByClass
        myClipPolygon = QgsGeometry.fromRect(myRect)
    elif (type(theExtent) is QgsGeometry and
          theExtent.wkbType in myAllowedClipTypes):
        myRect = theExtent.boundingBox().toRectF()
        myClipPolygon = theExtent
    else:
        raise InvalidClipGeometryError(
            tr(
                'Clip geometry must be an extent or a single part'
                'polygon based geometry.'))

    myProjectedExtent = myXForm.transformBoundingBox(myRect)

    # Get vector layer
    myProvider = theLayer.dataProvider()
    if myProvider is None:
        myMessage = tr('Could not obtain data provider from '
                       'layer "%s"' % theLayer.source())
        raise Exception(myMessage)

    # Get the layer field list, select by our extent then write to disk
    # .. todo:: FIXME - for different geometry types we should implement
    #    different clipping behaviour e.g. reject polygons that
    #    intersect the edge of the bbox. Tim
    myAttributes = myProvider.attributeIndexes()
    myFetchGeometryFlag = True
    myUseIntersectFlag = True
    myProvider.select(
        myAttributes,
        myProjectedExtent,
        myFetchGeometryFlag,
        myUseIntersectFlag)

    myFieldList = myProvider.fields()

    myWriter = QgsVectorFileWriter(
        myFilename,
        'UTF-8',
        myFieldList,
        theLayer.wkbType(),
        myGeoCrs,
        #'SQLite')  # FIXME (Ole): This works but is far too slow
        'ESRI Shapefile')
    if myWriter.hasError() != QgsVectorFileWriter.NoError:
        myMessage = tr('Error when creating shapefile: <br>Filename:'
                       '%s<br>Error: %s' %
                       (myFilename, myWriter.hasError()))
        raise Exception(myMessage)

    # Reverse the coordinate xform now so that we can convert
    # geometries from layer crs to geocrs.
    myXForm = QgsCoordinateTransform(theLayer.crs(), myGeoCrs)
    # Retrieve every feature with its geometry and attributes
    myFeature = QgsFeature()
    myCount = 0
    while myProvider.nextFeature(myFeature):
        myGeometry = myFeature.geometry()
        # Loop through the parts adding them to the output file
        # we write out single part features unless theExplodeFlag is False
        if theExplodeFlag:
            myGeometryList = explodeMultiPartGeometry(myGeometry)
        else:
            myGeometryList = [myGeometry]

        for myPart in myGeometryList:
            myPart.transform(myXForm)
            if theHardClipFlag:
                # Remove any dangling bits so only intersecting area is
                # kept.
                myPart = clipGeometry(myClipPolygon, myPart)
            if myPart is None:
                continue
            myFeature.setGeometry(myPart)
            myWriter.addFeature(myFeature)
        myCount += 1
    del myWriter  # Flush to disk

    if myCount < 1:
        myMessage = tr('No features fall within the clip extents. '
                       'Try panning / zooming to an area containing data '
                       'and then try to run your analysis again.'
                       'If hazard and exposure data doesn\'t overlap '
                       'at all, it is not possible to do an analysis.'
                       'Another possibility is that the layers do overlap '
                       'but because they may have different spatial '
                       'references, they appear to be disjoint. '
                       'If this is the case, try to turn on reproject '
                       'on-the-fly in QGIS.')
        raise NoFeaturesInExtentError(myMessage)

    myKeywordIO = KeywordIO()
    myKeywordIO.copyKeywords(
        theLayer, myFilename, theExtraKeywords=theExtraKeywords)

    return myFilename  # Filename of created file

# 需要导入模块: from safe_qgis.keyword_io import KeywordIO [as 别名]
# 或者: from safe_qgis.keyword_io.KeywordIO import copyKeywords [as 别名]

#.........这里部分代码省略.........
       None

    """
    if not theLayer or not theExtent:
        myMessage = tr('Layer or Extent passed to clip is None.')
        raise InvalidParameterException(myMessage)

    if theLayer.type() != QgsMapLayer.VectorLayer:
        myMessage = tr('Expected a vector layer but received a %s.' %
                       str(theLayer.type()))
        raise InvalidParameterException(myMessage)

    #myHandle, myFilename = tempfile.mkstemp('.sqlite', 'clip_',
    #    temp_dir())
    myHandle, myFilename = tempfile.mkstemp('.shp', 'clip_',
                                            temp_dir())

    # Ensure the file is deleted before we try to write to it
    # fixes windows specific issue where you get a message like this
    # ERROR 1: c:\temp\inasafe\clip_jpxjnt.shp is not a directory.
    # This is because mkstemp creates the file handle and leaves
    # the file open.
    os.close(myHandle)
    os.remove(myFilename)

    # Get the clip extents in the layer's native CRS
    myGeoCrs = QgsCoordinateReferenceSystem()
    myGeoCrs.createFromId(4326, QgsCoordinateReferenceSystem.EpsgCrsId)
    myXForm = QgsCoordinateTransform(myGeoCrs, theLayer.crs())
    myRect = QgsRectangle(theExtent[0], theExtent[1],
        theExtent[2], theExtent[3])
    myProjectedExtent = myXForm.transformBoundingBox(myRect)

    # Get vector layer
    myProvider = theLayer.dataProvider()
    if myProvider is None:
        myMessage = tr('Could not obtain data provider from '
                       'layer "%s"' % theLayer.source())
        raise Exception(myMessage)

    # Get the layer field list, select by our extent then write to disk
    # .. todo:: FIXME - for different geometry types we should implement
    #    different clipping behaviour e.g. reject polygons that
    #    intersect the edge of the bbox. Tim
    myAttributes = myProvider.attributeIndexes()
    myFetchGeometryFlag = True
    myUseIntersectFlag = True
    myProvider.select(myAttributes,
        myProjectedExtent,
        myFetchGeometryFlag,
        myUseIntersectFlag)

    myFieldList = myProvider.fields()

    myWriter = QgsVectorFileWriter(myFilename,
        'UTF-8',
        myFieldList,
        theLayer.wkbType(),
        myGeoCrs,
        #'SQLite')  # FIXME (Ole): This works but is far too slow
        'ESRI Shapefile')
    if myWriter.hasError() != QgsVectorFileWriter.NoError:
        myMessage = tr('Error when creating shapefile: <br>Filename:'
                       '%s<br>Error: %s' %
                       (myFilename, myWriter.hasError()))
        raise Exception(myMessage)

    # Reverse the coordinate xform now so that we can convert
    # geometries from layer crs to geocrs.
    myXForm = QgsCoordinateTransform(theLayer.crs(), myGeoCrs)
    # Retrieve every feature with its geometry and attributes
    myFeature = QgsFeature()
    myCount = 0
    while myProvider.nextFeature(myFeature):
        myGeometry = myFeature.geometry()
        # Loop through the parts adding them to the output file
        # we write out single part features unless explodeMultipart is False
        if explodeMultipart:
            myGeometryList = explodeMultiPartGeometry(myGeometry)
        else:
            myGeometryList = [myGeometry]

        for myPart in myGeometryList:
            myPart.transform(myXForm)
            myFeature.setGeometry(myPart)
            myWriter.addFeature(myFeature)
        myCount += 1
    del myWriter  # Flush to disk

    if myCount < 1:
        myMessage = tr('No features fall within the clip extents. '
                       'Try panning / zooming to an area containing data '
                       'and then try to run your analysis again.')
        raise NoFeaturesInExtentException(myMessage)

    myKeywordIO = KeywordIO()
    myKeywordIO.copyKeywords(theLayer, myFilename,
        theExtraKeywords=theExtraKeywords)

    return myFilename  # Filename of created file

# 需要导入模块: from safe_qgis.keyword_io import KeywordIO [as 别名]
# 或者: from safe_qgis.keyword_io.KeywordIO import copyKeywords [as 别名]
        myMessage = tr('<p>Error while executing the following shell command:'
                     '</p><pre>%s</pre><p>Error message: %s'
                     % (myCommand, str(e)))
        # shameless hack - see https://github.com/AIFDR/inasafe/issues/141
        if sys.platform == 'darwin':  # Mac OS X
            if 'Errno 4' in str(e):
                # continue as the error seems to be non critical
                pass
            else:
                raise Exception(myMessage)
        else:
            raise Exception(myMessage)

    # .. todo:: Check the result of the shell call is ok
    myKeywordIO = KeywordIO()
    myKeywordIO.copyKeywords(theLayer, myFilename,
                             theExtraKeywords=theExtraKeywords)
    return myFilename  # Filename of created file


def extentToKml(theExtent):
    """A helper to get a little kml doc for an extent so that
    we can use it with gdal warp for clipping."""

    myBottomLeftCorner = '%f,%f' % (theExtent[0], theExtent[1])
    myTopLeftCorner = '%f,%f' % (theExtent[0], theExtent[3])
    myTopRightCorner = '%f,%f' % (theExtent[2], theExtent[3])
    myBottomRightCorner = '%f,%f' % (theExtent[2], theExtent[1])
    myKml = ("""<?xml version="1.0" encoding="utf-8" ?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Document>
    <Folder>