Python clipper.clip_layer函数代码示例

    def test_clip_raster_small(self):
        """Raster layers can be clipped in small and precise size. For #710."""

        # Create a raster layer
        layer_name = 'shake'
        raster_layer = QgsRasterLayer(RASTERPATH, layer_name)

        message = (
            'Did not find layer "%s" in path "%s"' % (layer_name, RASTERPATH))
        assert raster_layer is not None, message

        # Create a bounding box
        bounding_box = [97, -3, 104, 1]

        # Clip the vector to the bbox
        result = clip_layer(raster_layer, bounding_box)

        # Check the output is valid
        assert os.path.exists(result.source())

        # Clip and give a desired resolution for the output

        # small pixel size and high precision

        # based on pixel size of Flood_Current_Depth_Jakarta_geographic.asc
        size = 0.00045228819716
        result = clip_layer(raster_layer, bounding_box, size)
        new_raster_layer = QgsRasterLayer(result.source(), layer_name)
        assert new_raster_layer.isValid(), 'Resampled raster is not valid'
        message = (
            'Resampled raster has incorrect pixel size. Expected: %.14f, '
            'Actual: %.14f' % (
                size, new_raster_layer.rasterUnitsPerPixelX()))
        result_size = new_raster_layer.rasterUnitsPerPixelX()
        self.assertAlmostEqual(result_size, size, places=13, msg=message)

    def test_clipRaster(self):
        """Raster layers can be clipped
        """

        # Create a raster layer
        myName = 'shake'
        myRasterLayer = QgsRasterLayer(RASTERPATH, myName)

        myMessage = 'Did not find layer "%s" in path "%s"' % \
                    (myName, RASTERPATH)
        assert myRasterLayer is not None, myMessage

        # Create a bounding box
        myRect = [97, -3, 104, 1]

        # Clip the vector to the bbox
        myResult = clip_layer(myRasterLayer, myRect)

        # Check the output is valid
        assert os.path.exists(myResult.source())

        # Clip and give a desired resolution for the output
        mySize = 0.05
        myResult = clip_layer(myRasterLayer, myRect, mySize)
        myNewRasterLayer = QgsRasterLayer(myResult.source(), myName)
        assert myNewRasterLayer.isValid(), 'Resampled raster is not valid'

        myMessage = ('Resampled raster has incorrect pixel size.'
                     'Expected: %f, Actual: %f' %
                     (mySize, myNewRasterLayer.rasterUnitsPerPixelX()))
        assert myNewRasterLayer.rasterUnitsPerPixelX() == mySize, myMessage

    def test_clip_raster(self):
        """Raster layers can be clipped."""

        # Create a raster layer
        layer_name = 'shake'
        raster_layer = QgsRasterLayer(RASTERPATH, layer_name)

        message = (
            'Did not find layer "%s" in path "%s"' % (layer_name, RASTERPATH))
        assert raster_layer is not None, message

        # Create a bounding box
        bounding_box = [97, -3, 104, 1]

        # Clip the vector to the bbox
        result = clip_layer(raster_layer, bounding_box)

        # Check the output is valid
        assert os.path.exists(result.source())

        # Clip and give a desired resolution for the output
        # big pixel size
        size = 0.05
        result = clip_layer(raster_layer, bounding_box, size)
        new_raster_layer = QgsRasterLayer(result.source(), layer_name)
        assert new_raster_layer.isValid(), 'Resampled raster is not valid'
        message = (
            'Resampled raster has incorrect pixel size. Expected: %5f, '
            'Actual: %5f' % (size, new_raster_layer.rasterUnitsPerPixelX()))
        assert new_raster_layer.rasterUnitsPerPixelX() == size, message

    def test_raster_scaling_projected(self):
        """Attempt to scale projected density raster layers raise exception.

        Automatic scaling when resampling density data
        does not currently work for projected layers. See issue #123.

        For the time being this test checks that an exception is raised
        when scaling is attempted on projected layers.
        When we resolve issue #123, this test should be rewritten.
        """

        test_filename = 'Population_Jakarta_UTM48N.tif'
        raster_path = ('%s/%s' % (TESTDATA, test_filename))

        # Get reference values
        safe_layer = read_safe_layer(raster_path)
        min_value, max_value = safe_layer.get_extrema()
        native_resolution = safe_layer.get_resolution()

        print min_value, max_value
        print native_resolution

        # Define bounding box in EPSG:4326
        bounding_box = [106.61, -6.38, 107.05, -6.07]

        resolutions = [0.02, 0.01, 0.005, 0.002, 0.001]
        # Test for a range of resolutions
        for resolution in resolutions:

            # Clip the raster to the bbox
            extra_keywords = {'resolution': native_resolution}
            raster_layer = QgsRasterLayer(raster_path, 'xxx')
            try:
                clip_layer(
                    raster_layer,
                    bounding_box,
                    resolution,
                    extra_keywords=extra_keywords
                )
            except InvalidProjectionError:
                pass
            else:
                message = 'Should have raised InvalidProjectionError'
                raise Exception(message)

    def testRasterScaling_projected(self):
        """Attempt to scale projected density raster layers raise exception

        Automatic scaling when resampling density data
        does not currently work for projected layers. See issue #123.

        For the time being this test checks that an exception is raised
        when scaling is attempted on projected layers.
        When we resolve issue #123, this test should be rewritten.
        """

        myTestFilename = 'Population_Jakarta_UTM48N.tif'
        myRasterPath = ('%s/%s' % (TESTDATA, myTestFilename))

        # Get reference values
        mySafeLayer = readSafeLayer(myRasterPath)
        myMinimum, myMaximum = mySafeLayer.get_extrema()
        myNativeResolution = mySafeLayer.get_resolution()

        print
        print myMinimum, myMaximum
        print myNativeResolution

        # Define bounding box in EPSG:4326
        myBoundingBox = [106.61, -6.38, 107.05, -6.07]

        # Test for a range of resolutions
        for myResolution in [0.02, 0.01, 0.005, 0.002, 0.001]:

            # Clip the raster to the bbox
            myExtraKeywords = {'resolution': myNativeResolution}
            myRasterLayer = QgsRasterLayer(myRasterPath, 'xxx')
            try:
                clip_layer(myRasterLayer,
                           myBoundingBox,
                           myResolution,
                           extra_keywords=myExtraKeywords)
            except InvalidProjectionError:
                pass
            else:
                myMessage = 'Should have raised InvalidProjectionError'
                raise Exception(myMessage)

    def test_clipOnePixel(self):
        # Create a raster layer
        myRasterPath = os.path.join(EXPDATA, 'glp10ag.asc')
        myTitle = 'people'
        myRasterLayer = QgsRasterLayer(myRasterPath, myTitle)

        # Create very small bounding box
        mySmallRect = [100.3430, -0.9089, 100.3588, -0.9022]

        # Clip the raster to the bbox
        try:
            _ = clip_layer(myRasterLayer, mySmallRect)
        except CallGDALError:
            pass
        except Exception, e:
            raise Exception('Exception is not expected, %s' % e)

 def test_vector_projections(self):
     """Test that vector input data is reprojected properly during clip."""
     # Input data is OSM in GOOGLE CRS
     # We are reprojecting to GEO and expecting the output shp to be in GEO
     # see https://github.com/AIFDR/inasafe/issues/119
     # and https://github.com/AIFDR/inasafe/issues/95
     vector_layer = QgsVectorLayer(
         VECTOR_PATH2, 'OSM Buildings', 'ogr')
     message = 'Failed to load osm buildings'
     assert vector_layer is not None, message
     assert vector_layer.isValid()
     set_canvas_crs(GEOCRS, True)
     set_jakarta_extent()
     clip_rectangle = [106.52, -6.38, 107.14, -6.07]
     # Clip the vector to the bbox
     result = clip_layer(vector_layer, clip_rectangle)
     assert(os.path.exists(result.source()))

    def test_clip_one_pixel(self):
        """Test for clipping less than one pixel."""
        # Create a raster layer
        raster_path = os.path.join(EXPDATA, 'glp10ag.asc')
        title = 'people'
        raster_layer = QgsRasterLayer(raster_path, title)

        # Create very small bounding box
        small_bounding_box = [100.3430, -0.9089, 100.3588, -0.9022]

        # Clip the raster to the bbox
        try:
            _ = clip_layer(raster_layer, small_bounding_box)
        except CallGDALError:
            pass
        except Exception, e:
            raise Exception('Exception is not expected, %s' % e)

    def test_clip_vector_hard(self):
        """Vector layers can be hard clipped.

        Hard clipping will remove any dangling, non intersecting elements.
        """
        vector_layer = QgsVectorLayer(
            VECTOR_PATH3, 'OSM Buildings', 'ogr')
        message = 'Failed to load osm buildings'
        assert vector_layer is not None, message
        assert vector_layer.isValid()
        set_canvas_crs(GEOCRS, True)
        set_jakarta_extent()
        clip_rectangle = [106.8218, -6.1842, 106.8232, -6.1830]

        # Clip the vector to the bbox
        result = clip_layer(vector_layer, clip_rectangle, hard_clip_flag=True)

        # Check the output is valid
        assert(os.path.exists(result.source()))

    def test_clip_vector(self):
        """Vector layers can be clipped."""

        # Create a vector layer
        layer_name = 'padang'
        vector_layer = QgsVectorLayer(VECTOR_PATH, layer_name, 'ogr')

        message = (
            'Did not find layer "%s" in path "%s"' % (layer_name, VECTOR_PATH))
        assert vector_layer is not None, message
        assert vector_layer.isValid()
        # Create a bounding box
        bounding_box = [100.03, -1.14, 100.81, -0.73]

        # Clip the vector to the bbox
        result = clip_layer(vector_layer, bounding_box)

        # Check the output is valid
        assert(os.path.exists(result.source()))

    def test_clipVector(self):
        """Vector layers can be clipped.
        """

        # Create a vector layer
        myName = 'padang'
        myVectorLayer = QgsVectorLayer(VECTOR_PATH, myName, 'ogr')

        myMessage = 'Did not find layer "%s" in path "%s"' % \
                    (myName, VECTOR_PATH)
        assert myVectorLayer is not None, myMessage
        assert myVectorLayer.isValid()
        # Create a bounding box
        myRect = [100.03, -1.14, 100.81, -0.73]

        # Clip the vector to the bbox
        myResult = clip_layer(myVectorLayer, myRect)

        # Check the output is valid
        assert(os.path.exists(myResult.source()))

    def _aggregate(self, myImpactLayer, myExpectedResults):
        myAggregationLayer = QgsVectorLayer(
            os.path.join(BOUNDDATA, 'kabupaten_jakarta.shp'),
            'test aggregation',
            'ogr')
        # create a copy of aggregation layer
        myGeoExtent = extent_to_geo_array(
            myAggregationLayer.extent(),
            myAggregationLayer.crs())

        myAggrAttribute = self.keywordIO.read_keywords(
            myAggregationLayer, self.defaults['AGGR_ATTR_KEY'])
        # noinspection PyArgumentEqualDefault
        myAggregationLayer = clip_layer(
            layer=myAggregationLayer,
            extent=myGeoExtent,
            explode_flag=True,
            explode_attribute=myAggrAttribute)

        myAggregator = Aggregator(None, myAggregationLayer)
        # setting up
        myAggregator.isValid = True
        myAggregator.layer = myAggregationLayer
        myAggregator.safeLayer = safe_read_layer(
            str(myAggregator.layer.source()))
        myAggregator.aoiMode = False
        myAggregator.aggregate(myImpactLayer)

        myProvider = myAggregator.layer.dataProvider()
        myProvider.select(myProvider.attributeIndexes())
        myFeature = QgsFeature()
        myResults = []

        while myProvider.nextFeature(myFeature):
            myFeatureResults = {}
            myAtMap = myFeature.attributeMap()
            for (k, attr) in myAtMap.iteritems():
                myFeatureResults[k] = str(attr.toString())
            myResults.append(myFeatureResults)

        self.assertEqual(myExpectedResults, myResults)

    def test_clip_raster_with_no_extension(self):
        """Test we can clip a raster with no extension - see #659."""
        # Create a raster layer
        source_file = os.path.join(UNITDATA, 'other', 'tenbytenraster.asc')
        test_file = unique_filename(prefix='tenbytenraster-')
        shutil.copyfile(source_file, test_file)

        # Create a keywords file
        source_file = os.path.join(
            UNITDATA, 'other', 'tenbytenraster.keywords')
        keywords_file = test_file + '.keywords'
        shutil.copyfile(source_file, keywords_file)

        # Test the raster layer
        raster_layer = QgsRasterLayer(test_file, 'ten by ten')
        # Create a bounding box
        rectangle = [1535380, 5083260, 1535380 + 40, 5083260 + 40]
        # Clip the vector to the bounding box
        result = clip_layer(raster_layer, rectangle)
        # Check the output is valid
        assert os.path.exists(result.source())

    def _prepareLayer(self):
        """Prepare the aggregation layer to match analysis extents."""
        myMessage = m.Message(
            m.Heading(
                self.tr('Preparing aggregation layer'),
                **PROGRESS_UPDATE_STYLE),
            m.Paragraph(self.tr(
                'We are clipping the aggregation layer to match the '
                'intersection of the hazard and exposure layer extents.')))
        self._sendMessage(myMessage)

        # This is used to hold an *in memory copy* of the aggregation layer
        # or a in memory layer with the clip extents as a feature.
        if self.aoiMode:
            self.layer = self._extentsToLayer()
            # Area Of Interest (AOI) mode flag
        else:
            # we use only the exposure extent, because both exposure and hazard
            # have the same extent at this point.
            myGeoExtent = extent_to_geo_array(
                self.exposureLayer.extent(),
                self.exposureLayer.crs())

            myAggrAttribute = self.keywordIO.read_keywords(
                self.layer, self.defaults['AGGR_ATTR_KEY'])

            myClippedLayer = clip_layer(
                layer=self.layer,
                extent=myGeoExtent,
                explode_flag=True,
                explode_attribute=myAggrAttribute)

            myName = '%s %s' % (self.layer.name(), self.tr('aggregation'))
            self.layer = myClippedLayer
            self.layer.setLayerName(myName)
            if self.showIntermediateLayers:
                self.keywordIO.update_keywords(self.layer, {'title': myName})
                QgsMapLayerRegistry.instance().addMapLayer(self.layer)

    def testRasterScaling(self):
        """Raster layers can be scaled when resampled

        This is a test for ticket #52

        Native test .asc data has

        Population_Jakarta_geographic.asc
        ncols         638
        nrows         649
        cellsize      0.00045228819716044

        Population_2010.asc
        ncols         5525
        nrows         2050
        cellsize      0.0083333333333333

        Scaling is necessary for raster data that represents density
        such as population per km^2
        """

        for myFilename in ['Population_Jakarta_geographic.asc',
                           'Population_2010.asc']:

            myRasterPath = ('%s/%s' % (TESTDATA, myFilename))

            # Get reference values
            mySafeLayer = readSafeLayer(myRasterPath)
            myMinimum, myMaximum = mySafeLayer.get_extrema()
            del myMaximum
            del myMinimum
            myNativeResolution = mySafeLayer.get_resolution()

            # Get the Hazard extents as an array in EPSG:4326
            myBoundingBox = mySafeLayer.get_bounding_box()

            # Test for a range of resolutions
            for myResolution in [0.02,
                                 0.01,
                                 0.005,
                                 0.002,
                                 0.001,
                                 0.0005,   # Coarser
                                 0.0002]:  # Finer

                # To save time only do two resolutions for the
                # large population set
                if myFilename.startswith('Population_2010'):
                    if myResolution > 0.01 or myResolution < 0.005:
                        break

                # Clip the raster to the bbox
                myExtraKeywords = {'resolution': myNativeResolution}
                myRasterLayer = QgsRasterLayer(myRasterPath, 'xxx')
                myResult = clip_layer(myRasterLayer,
                                      myBoundingBox,
                                      myResolution,
                                      extra_keywords=myExtraKeywords)

                mySafeLayer = readSafeLayer(myResult.source())
                myNativeData = mySafeLayer.get_data(scaling=False)
                myScaledData = mySafeLayer.get_data(scaling=True)

                mySigma = (mySafeLayer.get_resolution()[0] /
                           myNativeResolution[0]) ** 2

                # Compare extrema
                myExpectedScaledMax = mySigma * numpy.nanmax(myNativeData)
                myMessage = ('Resampled raster was not rescaled correctly: '
                             'max(myScaledData) was %f but expected %f'
                             % (numpy.nanmax(myScaledData),
                                myExpectedScaledMax))

                # FIXME (Ole): The rtol used to be 1.0e-8 -
                #              now it has to be 1.0e-6, otherwise we get
                #              max(myScaledData) was 12083021.000000 but
                #              expected 12083020.414316
                #              Is something being rounded to the nearest
                #              integer?
                assert numpy.allclose(myExpectedScaledMax,
                                      numpy.nanmax(myScaledData),
                                      rtol=1.0e-6, atol=1.0e-8), myMessage

                myExpectedScaledMin = mySigma * numpy.nanmin(myNativeData)
                myMessage = ('Resampled raster was not rescaled correctly: '
                             'min(myScaledData) was %f but expected %f'
                             % (numpy.nanmin(myScaledData),
                                myExpectedScaledMin))
                assert numpy.allclose(myExpectedScaledMin,
                                      numpy.nanmin(myScaledData),
                                      rtol=1.0e-8, atol=1.0e-12), myMessage

                # Compare elementwise
                myMessage = 'Resampled raster was not rescaled correctly'
                assert nan_allclose(myNativeData * mySigma, myScaledData,
                                   rtol=1.0e-8, atol=1.0e-8), myMessage

                # Check that it also works with manual scaling
                myManualData = mySafeLayer.get_data(scaling=mySigma)
                myMessage = 'Resampled raster was not rescaled correctly'
#.........这里部分代码省略.........