Source code for safe.impact_statistics.test.test_zonal_stats

# coding=utf-8
"""
InaSAFE Disaster risk assessment tool developed by AusAid -
 **Zonal Tests.**

Contact : [email protected]

.. note:: This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation; either version 2 of the License, or
     (at your option) any later version.

"""
__author__ = '[email protected]'
__date__ = '17/10/2013'
__copyright__ = ('Copyright 2012, Australia Indonesia Facility for '
                 'Disaster Reduction')
import unittest
import os

from qgis.core import QgsRectangle

from safe.impact_statistics.zonal_stats import (
    calculate_zonal_stats, intersection_box)
from safe.test.utilities import (
    load_layer,
    test_data_path,
    get_qgis_app)

QGIS_APP, CANVAS, IFACE, PARENT = get_qgis_app()


[docs]class ZonalStatsTest(unittest.TestCase): """Tests for zonal related functions. # Note there is a weird issue on OSX where # geometry comes back as none all the time when testing # running the test again will randomly succeed in fetching geometries # properly. """ # noinspection PyPep8Naming
[docs] def setUp(self): os.environ['LANG'] = 'en' # noinspection PyPep8Naming
[docs] def tearDown(self): pass
[docs] def test_zonal(self): """Test that zonal stats returns the expected output.""" raster_layer, _ = load_layer( test_data_path('other', 'tenbytenraster.asc')) vector_layer, _ = load_layer( test_data_path('other', 'zonal_polygons.shp')) result = calculate_zonal_stats( raster_layer=raster_layer, polygon_layer=vector_layer) expected_result = { 0L: {'count': 4, 'sum': 34.0, 'mean': 8.5}, # BR polygon 1L: {'count': 9, 'sum': 36.0, 'mean': 4.0}, # center polygon 2L: {'count': 4, 'sum': 2.0, 'mean': 0.5}, # TL polygon 3L: {'count': 4, 'sum': 2.0, 'mean': 0.5}, # BL Polygon 4L: {'count': 4, 'sum': 34.0, 'mean': 8.5}} # TR polygon # noinspection PyPep8Naming self.maxDiff = None self.assertDictEqual(expected_result, result)
[docs] def test_zonal_with_exact_cell_boundaries(self): """Test that zonal stats returns the expected output.""" raster_layer, _ = load_layer( test_data_path('other', 'tenbytenraster.asc')) # Note this is a matrix of 11x11 polygons - one per cell # and one poly extending beyond to the right of each row # and one poly extending beyond the bottom of each col vector_layer, _ = load_layer( test_data_path('other', 'ten_by_ten_raster_as_polys.shp')) result = calculate_zonal_stats( raster_layer=raster_layer, polygon_layer=vector_layer) expected_result = { 0L: {'count': 1.0, 'sum': 0.0, 'mean': 0.0}, # TL polygon 9L: {'count': 1.0, 'sum': 9.0, 'mean': 9.0}, # TR polygon 25L: {'count': 1.0, 'sum': 3.0, 'mean': 3.0}, # Central polygon 88L: {'count': 1.0, 'sum': 0.0, 'mean': 0.0}, # BL polygon 108L: {'count': 1.0, 'sum': 9.0, 'mean': 9.0}} # BR polygon # We will just check TL, TR, Middle, BL and BR cells result = { 0L: result[0L], 9L: result[9L], 25L: result[25L], 88L: result[88L], 108L: result[108L]} # noinspection PyPep8Naming self.maxDiff = None self.assertDictEqual(expected_result, result)
[docs] def test_cell_info_for_bbox(self): """Test that cell info for bbox returns expected values.""" raster_box = QgsRectangle(1535375.0, 5083255.0, 1535475.0, 5083355.0) feature_box = QgsRectangle(1535455.0, 5083345.0, 1535465.0, 5083355.0) cell_size_x = 10 cell_size_y = 10 offset_x = 8 offset_y = 0 cells_x = cells_y = 1 expected_result = offset_x, offset_y, cells_x, cells_y result = intersection_box( raster_box, feature_box, cell_size_x, cell_size_y) self.assertTupleEqual(expected_result, result) # Now try with a poly bbox slightly larger than 1 px (1535470) cells_x = 2 expected_result = offset_x, offset_y, cells_x, cells_y feature_box = QgsRectangle(1535455.0, 5083345.0, 1535470.0, 5083355.0) result = intersection_box( raster_box, feature_box, cell_size_x, cell_size_y) self.assertTupleEqual(expected_result, result)