Source code for safe_qgis.test_clipper
"""
InaSAFE Disaster risk assessment tool developed by AusAid -
**ISClipper test suite.**
Contact : ole.moller.nielsen@gmail.com
.. 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.
"""
from unittest import expectedFailure
__author__ = 'tim@linfiniti.com'
__date__ = '20/01/2011'
__copyright__ = ('Copyright 2012, Australia Indonesia Facility for '
'Disaster Reduction')
import unittest
import sys
import os
# Add PARENT directory to path to make test aware of other modules
pardir = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
sys.path.append(pardir)
import numpy
from qgis.core import (QgsVectorLayer,
QgsRasterLayer,
QgsGeometry,
QgsPoint)
from safe_qgis.safe_interface import readSafeLayer
from safe_qgis.safe_interface import getOptimalExtent
from safe_qgis.exceptions import InvalidProjectionError, CallGDALError
from safe_qgis.clipper import (clipLayer,
extentToKml,
explodeMultiPartGeometry,
clipGeometry)
from safe_qgis.utilities import qgisVersion
from safe_qgis.utilities_test import (getQgisTestApp,
setCanvasCrs,
RedirectStdStreams,
DEVNULL,
GEOCRS,
setJakartaGeoExtent)
from safe.common.testing import HAZDATA, TESTDATA, EXPDATA, UNITDATA
from safe.common.exceptions import GetDataError
from safe.api import nanallclose
# Setup pathnames for test data sets
VECTOR_PATH = os.path.join(TESTDATA, 'Padang_WGS84.shp')
VECTOR_PATH2 = os.path.join(TESTDATA, 'OSM_subset_google_mercator.shp')
VECTOR_PATH3 = os.path.join(UNITDATA, 'exposure', 'buildings_osm_4326.shp')
RASTERPATH = os.path.join(HAZDATA, 'Shakemap_Padang_2009.asc')
RASTERPATH2 = os.path.join(TESTDATA, 'population_padang_1.asc')
# Handle to common QGis test app
QGISAPP, CANVAS, IFACE, PARENT = getQgisTestApp()
[docs]class ClipperTest(unittest.TestCase):
"""Test the InaSAFE clipper"""
def setUp(self):
pass
def tearDown(self):
pass
[docs] 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 = clipLayer(myVectorLayer, myRect)
# Check the output is valid
assert(os.path.exists(myResult))
[docs] 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 = clipLayer(myRasterLayer, myRect)
# Check the output is valid
assert os.path.exists(myResult)
# Clip and give a desired resolution for the output
mySize = 0.05
myResult = clipLayer(myRasterLayer, myRect, mySize)
myNewRasterLayer = QgsRasterLayer(myResult, myName)
assert myNewRasterLayer.isValid(), 'Resampled raster is not valid'
myMessage = ('Resampled raster has incorrect pixel size.'
'Expected: %f, Actual: %f' %
(mySize, myNewRasterLayer.rasterUnitsPerPixel()))
assert myNewRasterLayer.rasterUnitsPerPixel() == mySize, myMessage
# See issue #349
@expectedFailure
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:
myResult = clipLayer(myRasterLayer, mySmallRect)
except CallGDALError:
pass
except Exception, e:
raise Exception('Exception is not expected, %s' % e)
else:
myMsg = "Failed, does not raise exception"
raise Exception(myMsg)
# Create not very small bounding box
myNotSmallRect = [100.3430, -0.9089, 101.3588, -1.9022]
# Clip the raster to the bbox
myResult = clipLayer(myRasterLayer, myNotSmallRect)
assert myResult is not None, 'Should be a success clipping'
[docs] def test_invalidFilenamesCaught(self):
"""Invalid filenames raise appropriate exceptions
Wrote this test because test_clipBoth raised the wrong error
when file was missing. Instead of reporting that, it gave
Western boundary must be less than eastern. I got [0.0, 0.0, 0.0, 0.0]
See issue #170
"""
# Try to create a vector layer from non-existing filename
myName = 'stnhaoeu_78oeukqjkrcgA'
myPath = 'OEk_tnshoeu_439_kstnhoe'
with RedirectStdStreams(stdout=DEVNULL, stderr=DEVNULL):
myVectorLayer = QgsVectorLayer(myPath, myName, 'ogr')
myMessage = ('QgsVectorLayer reported "valid" for non '
'existent path "%s" and name "%s".'
% (myPath, myName))
assert not myVectorLayer.isValid(), myMessage
# Create a raster layer
with RedirectStdStreams(stdout=DEVNULL, stderr=DEVNULL):
myRasterLayer = QgsRasterLayer(myPath, myName)
myMessage = ('QgsRasterLayer reported "valid" for non '
'existent path "%s" and name "%s".'
% (myPath, myName))
assert not myRasterLayer.isValid(), myMessage
[docs] def test_clipBoth(self):
"""Raster and 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.isValid(), myMessage
# Create a raster layer
myName = 'shake'
myRasterLayer = QgsRasterLayer(RASTERPATH, myName)
myMessage = 'Did not find layer "%s" in path "%s"' % (myName,
RASTERPATH)
assert myRasterLayer.isValid(), myMessage
# Create a bounding box
myViewportGeoExtent = [99.53, -1.22, 101.20, -0.36]
# Get the Hazard extents as an array in EPSG:4326
myHazardGeoExtent = [myRasterLayer.extent().xMinimum(),
myRasterLayer.extent().yMinimum(),
myRasterLayer.extent().xMaximum(),
myRasterLayer.extent().yMaximum()]
# Get the Exposure extents as an array in EPSG:4326
myExposureGeoExtent = [myVectorLayer.extent().xMinimum(),
myVectorLayer.extent().yMinimum(),
myVectorLayer.extent().xMaximum(),
myVectorLayer.extent().yMaximum()]
# Now work out the optimal extent between the two layers and
# the current view extent. The optimal extent is the intersection
# between the two layers and the viewport.
# Extent is returned as an array [xmin,ymin,xmax,ymax]
myGeoExtent = getOptimalExtent(myHazardGeoExtent,
myExposureGeoExtent,
myViewportGeoExtent)
# Clip the vector to the bbox
myResult = clipLayer(myVectorLayer, myGeoExtent)
# Check the output is valid
assert os.path.exists(myResult)
readSafeLayer(myResult)
# Clip the raster to the bbox
myResult = clipLayer(myRasterLayer, myGeoExtent)
# Check the output is valid
assert os.path.exists(myResult)
readSafeLayer(myResult)
# -------------------------------
# Check the extra keywords option
# -------------------------------
# Clip the vector to the bbox
myResult = clipLayer(myVectorLayer, myGeoExtent,
theExtraKeywords={'kermit': 'piggy'})
# Check the output is valid
assert os.path.exists(myResult)
L = readSafeLayer(myResult)
kwds = L.get_keywords()
myMessage = 'Extra keyword was not found in %s: %s' % (myResult, kwds)
assert kwds['kermit'] == 'piggy'
# Clip the raster to the bbox
myResult = clipLayer(myRasterLayer, myGeoExtent,
theExtraKeywords={'zoot': 'animal'})
# Check the output is valid
assert os.path.exists(myResult)
L = readSafeLayer(myResult)
kwds = L.get_keywords()
myMessage = 'Extra keyword was not found in %s: %s' % (myResult, kwds)
assert kwds['zoot'] == 'animal', myMessage
[docs] 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 = clipLayer(myRasterLayer,
myBoundingBox,
myResolution,
theExtraKeywords=myExtraKeywords)
mySafeLayer = readSafeLayer(myResult)
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 nanallclose(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'
assert nanallclose(myManualData, myScaledData,
rtol=1.0e-8, atol=1.0e-8), myMessage
# Check that an exception is raised for bad arguments
try:
mySafeLayer.get_data(scaling='bad')
except GetDataError:
pass
else:
myMessage = 'String argument should have raised exception'
raise Exception(myMessage)
try:
mySafeLayer.get_data(scaling='(1, 3)')
except GetDataError:
pass
else:
myMessage = 'Tuple argument should have raised exception'
raise Exception(myMessage)
# Check None option without keyword datatype == 'density'
mySafeLayer.keywords['datatype'] = 'undefined'
myUnscaledData = mySafeLayer.get_data(scaling=None)
myMessage = 'Data should not have changed'
assert nanallclose(myNativeData, myUnscaledData,
rtol=1.0e-12, atol=1.0e-12), myMessage
# Try with None and density keyword
mySafeLayer.keywords['datatype'] = 'density'
myUnscaledData = mySafeLayer.get_data(scaling=None)
myMessage = 'Resampled raster was not rescaled correctly'
assert nanallclose(myScaledData, myUnscaledData,
rtol=1.0e-12, atol=1.0e-12), myMessage
mySafeLayer.keywords['datatype'] = 'counts'
myUnscaledData = mySafeLayer.get_data(scaling=None)
myMessage = 'Data should not have changed'
assert nanallclose(myNativeData, myUnscaledData,
rtol=1.0e-12, atol=1.0e-12), myMessage
[docs] 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:
clipLayer(myRasterLayer,
myBoundingBox,
myResolution,
theExtraKeywords=myExtraKeywords)
except InvalidProjectionError:
pass
else:
myMessage = 'Should have raised InvalidProjectionError'
raise Exception(myMessage)
[docs] def test_extentToKml(self):
"""Test if extent to KML is working."""
myExtent = [100.03, -1.14, 100.81, -0.73]
myKmlFilename = extentToKml(myExtent)
assert os.path.exists(myKmlFilename)
myFile = open(myKmlFilename)
myKml = myFile.read()
myFile.close()
myMessage = 'Generated KML was not as expected: %s' % myKml
assert '<?xml version' in myKml, myMessage
for myValue in myExtent:
assert str(myValue) in myKml, myMessage
[docs] def test_vectorProjections(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
myVectorLayer = QgsVectorLayer(VECTOR_PATH2,
'OSM Buildings',
'ogr')
myMessage = 'Failed to load osm buildings'
assert myVectorLayer is not None, myMessage
assert myVectorLayer.isValid()
setCanvasCrs(GEOCRS, True)
setJakartaGeoExtent()
myClipRect = [106.52, -6.38, 107.14, -6.07]
# Clip the vector to the bbox
myResult = clipLayer(myVectorLayer, myClipRect)
assert(os.path.exists(myResult))
[docs] def test_explodeMultiPolygonGeometry(self):
"""Test exploding POLY multipart to single part geometries works"""
myGeometry = QgsGeometry.fromWkt('MULTIPOLYGON(((-0.966314 0.445890,'
'-0.281133 0.555729,-0.092839 0.218369,-0.908780 0.035305,-0.966314'
' 0.445890)),((-0.906164 0.003923,-0.077148 0.197447,0.043151'
' -0.074533,-0.882628 -0.296824,-0.906164 0.003923)))')
myCollection = explodeMultiPartGeometry(myGeometry)
myMessage = 'Expected 2 parts from multipart polygon geometry'
assert len(myCollection) == 2, myMessage
[docs] def test_explodeMultiLineGeometry(self):
"""Test exploding LINES multipart to single part geometries works"""
myGeometry = QgsGeometry.fromWkt('MULTILINESTRING((-0.974159 0.526961,'
' -0.291594 0.644645), (-0.411893 0.728331,'
' -0.160834 0.568804))')
myCollection = explodeMultiPartGeometry(myGeometry)
myMessage = 'Expected 2 parts from multipart line geometry'
assert len(myCollection) == 2, myMessage
[docs] def test_explodeMultiPointGeometry(self):
"""Test exploding POINT multipart to single part geometries works"""
myGeometry = QgsGeometry.fromWkt('MULTIPOINT((-0.966314 0.445890),'
'(-0.281133 0.555729))')
myCollection = explodeMultiPartGeometry(myGeometry)
myMessage = 'Expected 2 parts from multipart point geometry'
assert len(myCollection) == 2, myMessage
[docs] def test_clipGeometry(self):
"""Test that we can clip a geometry using another geometry."""
myGeometry = QgsGeometry.fromPolyline([
QgsPoint(10, 10),
QgsPoint(20, 20),
QgsPoint(30, 30),
QgsPoint(40, 40)])
myClipPolygon = QgsGeometry.fromPolygon([[
QgsPoint(20, 20),
QgsPoint(20, 30),
QgsPoint(30, 30),
QgsPoint(30, 20),
QgsPoint(20, 20)]])
myResult = clipGeometry(myClipPolygon, myGeometry)
if qgisVersion() > 10800:
myExpectedWkt = 'LINESTRING(20.0 20.0, 30.0 30.0)'
else:
myExpectedWkt = ('LINESTRING(20.000000 20.000000, '
'30.000000 30.000000)')
# There should only be one feature that intersects this clip
# poly so this assertion should work.
self.assertEqual(myExpectedWkt, str(myResult.exportToWkt()))
# Now poly on poly clip test
myClipPolygon = QgsGeometry.fromWkt(
'POLYGON((106.8218 -6.1842,106.8232 -6.1842,'
'106.82304963 -6.18317148,106.82179736 -6.18314774,'
'106.8218 -6.1842))')
myGeometry = QgsGeometry.fromWkt(
'POLYGON((106.8216869 -6.1852067,106.8213233 -6.1843051,'
'106.8212891 -6.1835559,106.8222566 -6.1835184,'
'106.8227557 -6.1835076,106.8228588 -6.1851204,'
'106.8216869 -6.1852067))')
myResult = clipGeometry(myClipPolygon, myGeometry)
if qgisVersion() > 10800:
myExpectedWkt = (
'POLYGON((106.82179833 -6.18353616,106.8222566 -6.1835184,'
'106.8227557 -6.1835076,106.82279996 -6.1842,'
'106.8218 -6.1842,106.82179833 -6.18353616))')
else:
myExpectedWkt = (
'POLYGON((106.821798 -6.183536,106.822257 -6.183518,'
'106.822756 -6.183508,106.822800 -6.184200,'
'106.821800 -6.184200,106.821798 -6.183536))')
# There should only be one feature that intersects this clip
# poly so this assertion should work.
self.assertEqual(myExpectedWkt, str(myResult.exportToWkt()))
# Now point on poly test clip
myGeometry = QgsGeometry.fromWkt('POINT(106.82241 -6.18369)')
myResult = clipGeometry(myClipPolygon, myGeometry)
if qgisVersion() > 10800:
myExpectedWkt = 'POINT(106.82241 -6.18369)'
else:
myExpectedWkt = 'POINT(106.822410 -6.183690)'
# There should only be one feature that intersects this clip
# poly so this assertion should work.
self.assertEqual(myExpectedWkt, str(myResult.exportToWkt()))
[docs] def test_clipVectorHard(self):
"""Vector layers can be hard clipped.
Hard clipping will remove any dangling, non intersecting elements.
"""
myVectorLayer = QgsVectorLayer(VECTOR_PATH3,
'OSM Buildings',
'ogr')
myMessage = 'Failed to load osm buildings'
assert myVectorLayer is not None, myMessage
assert myVectorLayer.isValid()
setCanvasCrs(GEOCRS, True)
setJakartaGeoExtent()
myClipRect = [106.8218, -6.1842, 106.8232, -6.1830]
# Clip the vector to the bbox
myResult = clipLayer(myVectorLayer, myClipRect, theHardClipFlag=True)
# Check the output is valid
assert(os.path.exists(myResult))
if __name__ == '__main__':
suite = unittest.makeSuite(ClipperTest, 'test')
runner = unittest.TextTestRunner(verbosity=2)
runner.run(suite)
