Skip to content

Commit

Permalink
ISSUE-103 Address Landscape.io issues in pycoal/environment.py
Browse files Browse the repository at this point in the history
  • Loading branch information
lewismc committed Jun 11, 2017
1 parent 940bddd commit 644a1ea
Show file tree
Hide file tree
Showing 2 changed files with 67 additions and 66 deletions.
91 changes: 46 additions & 45 deletions pycoal/environment.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,132 +22,133 @@

class EnvironmentalCorrelation:


def __init__(self):
"""
Construct a new ``EnvironmentalCorrelation`` object.
"""
pass

def intersectProximity(self, miningFilename, vectorFilename, proximity, correlatedFilename):
def intersect_proximity(self, mining_filename, vector_filename, proximity, correlated_filename):
"""
Generate an environmental correlation image containing pixels from the
mining classified image detected within a given distance of features
within a vector layer.
Args:
miningImage (str): filename of the mining classified image
vectorLayer (str): filename of vector layer
mining_filename (str): filename of the mining classified image
vector_filename (str): filename of vector layer
proximity (float): distance in meters
correlatedImage (str): filename of the correlated image
correlated_filename (str): filename of the correlated image
"""

# get path and file names
outputDirectory = dirname(abspath(correlatedFilename))
miningName = splitext(basename(abspath(miningFilename)))[0]
vectorName = splitext(basename(abspath(vectorFilename)))[0]
output_directory = dirname(abspath(correlated_filename))
mining_name = splitext(basename(abspath(mining_filename)))[0]
vector_name = splitext(basename(abspath(vector_filename)))[0]

# rasterize the vector features to the same dimensions as the mining image
featureHeaderName = outputDirectory + '/' + miningName + '_' + vectorName + '.hdr'
self.createEmptyCopy(miningFilename, featureHeaderName)
featureImageName = featureHeaderName[:-4] + '.img'
self.rasterize(vectorFilename, featureImageName)
feature_header_name = output_directory + '/' + mining_name + '_' + vector_name + '.hdr'
self.create_empty_copy(mining_filename, feature_header_name)
feature_image_name = feature_header_name[:-4] + '.img'
self.rasterize(vector_filename, feature_image_name)

# generate a proximity map from the features
proximityHeaderName = outputDirectory + '/' + miningName + '_' + vectorName + '_proximity.hdr'
proximityImageName = proximityHeaderName[:-4] + '.img'
self.proximity(featureImageName, proximityImageName)
proximity_header_name = output_directory + '/' + mining_name + '_' + vector_name + '_proximity.hdr'
proximity_image_name = proximity_header_name[:-4] + '.img'
self.proximity(feature_image_name, proximity_image_name)

# load mining and proximity images and initialize environmental correlation array
miningImage = spectral.open_image(miningFilename)
proximityImage = spectral.open_image(proximityHeaderName)
correlatedImage = numpy.zeros(shape=miningImage.shape, dtype=numpy.uint16)
mining_image = spectral.open_image(mining_filename)
proximity_image = spectral.open_image(proximity_header_name)
correlated_image = numpy.zeros(shape=mining_image.shape, dtype=numpy.uint16)

# get average pixel size
if miningImage.metadata.get('map info')[10][-6:].lower() == 'meters':
xPixelSize = float(miningImage.metadata.get('map info')[5])
yPixelSize = float(miningImage.metadata.get('map info')[6])
pixelSize = (xPixelSize + yPixelSize) / 2
if mining_image.metadata.get('map info')[10][-6:].lower() == 'meters':
x_pixel_size = float(mining_image.metadata.get('map info')[5])
y_pixel_size = float(mining_image.metadata.get('map info')[6])
pixel_size = (x_pixel_size + y_pixel_size) / 2
else:
raise ValueError('Mining image units not in meters.')

# intersect features within proximity
for x in range(miningImage.shape[0]):
for y in range(miningImage.shape[1]):
if miningImage[x,y,0]==1 and proximityImage[x,y,0]*pixelSize<=proximity:
correlatedImage[x,y,0] = miningImage[x,y,0]
for x in range(mining_image.shape[0]):
for y in range(mining_image.shape[1]):
if mining_image[x,y,0]==1 and proximity_image[x,y,0]*pixel_size<=proximity:
correlated_image[x,y,0] = mining_image[x,y,0]

# save the environmental correlation image
spectral.io.envi.save_classification(
correlatedFilename,
correlatedImage,
class_names=miningImage.metadata.get('class names'),
correlated_filename,
correlated_image,
class_names=mining_image.metadata.get('class names'),
metadata={
'data ignore value': 0,
'description': 'COAL '+pycoal.version+' environmental correlation image.',
'map info': miningImage.metadata.get('map info')
'map info': mining_image.metadata.get('map info')
})

def createEmptyCopy(self, sourceFilename, destinationFilename):
def create_empty_copy(self, source_filename, destination_filename):
"""
Create an empty copy of a COAL classified image with the same size.
Args:
sourceFilename (str): filename of the source image
destinationFilename (str): filename of the destination image
source_filename (str): filename of the source image
destination_filename (str): filename of the destination image
"""

# open the source image
source = spectral.open_image(sourceFilename)
source = spectral.open_image(source_filename)

# create an empty array of the same dimensions
destination = numpy.zeros(shape=source.shape, dtype=numpy.uint16)

# save it with source metadata
spectral.io.envi.save_classification(
destinationFilename,
destination_filename,
destination,
class_names=['No data','Data'],
metadata={
'data ignore value': 0,
'map info': source.metadata.get('map info')
})

def rasterize(self, vectorFilename, featureFilename):
def rasterize(self, vector_filename, feature_filename):
"""
Burn features from a vector image onto a raster image.
Args:
vectorFilename (str): filename of the vector image
featureFilename (str): filename of the raster image
vector_filename (str): filename of the vector image
feature_filename (str): filename of the raster image
"""

# assume the layer has the same name as the image
layerName = splitext(basename(vectorFilename))[0]
layerName = splitext(basename(vector_filename))[0]

# convert vector features into nonzero pixels of the output file
returncode = call(['gdal_rasterize',
'-burn', '1',
'-l', layerName,
vectorFilename,
featureFilename])
vector_filename,
feature_filename])

# detect errors
if returncode != 0:
raise RuntimeError('Could not rasterize vector.')

def proximity(self, featureFilename, proximityFilename):
def proximity(self, feature_filename, proximity_filename):
"""
Generate a proximity map from the features.
Args:
featureFilename (str): filename of the feature image
proximityFilename (str): filename of the proximity image
feature_filename (str): filename of the feature image
proximity_filename (str): filename of the proximity image
"""

# generate an ENVI proximity map with georeferenced units
returncode = call(['gdal_proximity.py',
featureFilename,
proximityFilename,
feature_filename,
proximity_filename,
'-of', 'envi'])

# detect errors
Expand Down
42 changes: 21 additions & 21 deletions pycoal/tests/environment_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,34 +45,34 @@
# Clip and choose the input, clip, and output layers.

# test files for proximity intersection test
miningFilename = 'images/ang20150420t182050_corr_v1e_img_class_mining_cut.hdr'
vectorFilename = 'images/NHDFlowline_cut.shp'
mining_filename = 'images/ang20150420t182050_corr_v1e_img_class_mining_cut.hdr'
vector_filename = 'images/NHDFlowline_cut.shp'
proximity = 10.0
correlatedFilename = 'images/ang20150420t182050_corr_v1e_img_class_mining_cut_NHDFlowline_corr.hdr'
testFilename = 'images/ang20150420t182050_corr_v1e_img_class_mining_cut_NHDFlowline_corr_test.hdr'
correlated_filename = 'images/ang20150420t182050_corr_v1e_img_class_mining_cut_NHDFlowline_corr.hdr'
test_filename = 'images/ang20150420t182050_corr_v1e_img_class_mining_cut_NHDFlowline_corr_test.hdr'

# remove generated files
def _test_intersectProximity_teardown():
miningName = splitext(basename(abspath(miningFilename)))[0]
vectorName = splitext(basename(abspath(vectorFilename)))[0]
outputDirectory = 'images'
featureHeaderName = outputDirectory + '/' + miningName + '_' + vectorName + '.hdr'
featureImageName = featureHeaderName[:-4] + '.img'
proximityHeaderName = outputDirectory + '/' + miningName + '_' + vectorName + '_proximity.hdr'
proximityImageName = proximityHeaderName[:-4] + '.img'
testImageName = testFilename[:-4] + '.img'
_remove_files([featureHeaderName, featureImageName,
proximityHeaderName, proximityImageName,
testFilename, testImageName])
def _test_intersect_proximity_teardown():
mining_name = splitext(basename(abspath(mining_filename)))[0]
vector_name = splitext(basename(abspath(vector_filename)))[0]
output_directory = 'images'
feature_header_name = output_directory + '/' + mining_name + '_' + vector_name + '.hdr'
feature_image_name = feature_header_name[:-4] + '.img'
proximity_header_name = output_directory + '/' + mining_name + '_' + vector_name + '_proximity.hdr'
proximity_image_name = proximity_header_name[:-4] + '.img'
test_image_name = test_filename[:-4] + '.img'
_remove_files([feature_header_name, feature_image_name,
proximity_header_name, proximity_image_name,
test_filename, test_image_name])

# verify that proximity intersection produces expected results
@with_setup(None, _test_intersectProximity_teardown)
@with_setup(None, _test_intersect_proximity_teardown)
@skipIf(environ.get('CONTINUOUS_INTEGRATION'), 'Skip test because GDAL not installed on server.')
def test_intersectProximity():
def test_intersect_proximity():
ec = environment.EnvironmentalCorrelation()
ec.intersectProximity(miningFilename, vectorFilename, proximity, testFilename)
expected = spectral.open_image(correlatedFilename)
actual = spectral.open_image(testFilename)
ec.intersect_proximity(mining_filename, vector_filename, proximity, test_filename)
expected = spectral.open_image(correlated_filename)
actual = spectral.open_image(test_filename)
assert numpy.array_equal(expected.asarray(), actual.asarray())
assert actual.metadata.get('description') == 'COAL '+pycoal.version+' environmental correlation image.'
assert expected.metadata.get('class names') == actual.metadata.get('class names')
Expand Down

0 comments on commit 644a1ea

Please sign in to comment.