Skip to content

Commit

Permalink
Added 'show' map feature
Browse files Browse the repository at this point in the history
  • Loading branch information
sambowers committed Apr 6, 2018
1 parent 6136107 commit e886400
Show file tree
Hide file tree
Showing 2 changed files with 100 additions and 27 deletions.
10 changes: 5 additions & 5 deletions biota/IO.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ def outputGeoTiff(data, filename, geo_t, proj, output_dir = os.getcwd(), dtype =



def _buildMap(fig, ax, data, lat, lon, title ='', cbartitle = '', vmin = 10., vmax = 40., cmap = 'YlGn'):
def buildMap(fig, ax, data, lat, lon, title ='', cbartitle = '', vmin = None, vmax = None, cmap = None):
"""
Builds a standardised map for overviewFigure().
"""
Expand Down Expand Up @@ -149,20 +149,20 @@ def overviewFigure(tile_change, output = False, show = True):

# Plot a map of AGB at t1
ax1 = fig.add_subplot(2, 2, 1)
_buildMap(fig, ax1, AGB_t1, tile_change.lat, tile_change.lon, title = 'AGB %s'%str(tile_change.year_t1), cbartitle = 'tC/ha')
buildMap(fig, ax1, AGB_t1, tile_change.lat, tile_change.lon, title = 'AGB %s'%str(tile_change.year_t1), cbartitle = 'tC/ha', vmin = 10., vmax = 40., cmap = 'YlGn')

# Plot a map of AGB at t2
ax2 = fig.add_subplot(2, 2, 2, sharex = ax1, sharey = ax1)
_buildMap(fig, ax2, AGB_t2, tile_change.lat, tile_change.lon, title = 'AGB %s'%str(tile_change.year_t2), cbartitle = 'tC/ha')
buildMap(fig, ax2, AGB_t2, tile_change.lat, tile_change.lon, title = 'AGB %s'%str(tile_change.year_t2), cbartitle = 'tC/ha', vmin = 10., vmax = 40., cmap = 'YlGn')

# Plot a map of absolute AGB change
ax3 = fig.add_subplot(2, 2, 3, sharex = ax1, sharey = ax1)
_buildMap(fig, ax3, AGB_change, tile_change.lat, tile_change.lon, title = 'AGB change (%s-%s)'%(str(tile_change.data_t1.year),str(tile_change.year_t2)),
buildMap(fig, ax3, AGB_change, tile_change.lat, tile_change.lon, title = 'AGB change (%s-%s)'%(str(tile_change.data_t1.year),str(tile_change.year_t2)),
cbartitle = 'tC/ha', vmin = -10., vmax = 10., cmap = 'RdBu')

# Plot a map of % AGB change
ax4 = fig.add_subplot(2, 2, 4, sharex = ax1, sharey = ax1)
_buildMap(fig, ax4, change_code, tile_change.lat, tile_change.lon, title = 'Change type (%s-%s)'%(str(tile_change.data_t1.year),str(tile_change.year_t2)),
buildMap(fig, ax4, change_code, tile_change.lat, tile_change.lon, title = 'Change type (%s-%s)'%(str(tile_change.data_t1.year),str(tile_change.year_t2)),
vmin = 1., vmax = 6., cmap = 'Spectral')

plt.tight_layout()
Expand Down
117 changes: 95 additions & 22 deletions biota/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -332,7 +332,7 @@ def __rebin(self):
# Update number of looks
self.nLooks = self.nLooks * (self.downsample_factor ** 2)

def getMask(self, masked_px_count = False):
def getMask(self, masked_px_count = False, output = False, show = False):
"""
Loads the mask into a numpy array.
"""
Expand All @@ -347,9 +347,13 @@ def getMask(self, masked_px_count = False):
if self.downsample_factor != 1 and masked_px_count:
mask = skimage.measure.block_reduce(mask, (self.downsample_factor, self.downsample_factor), np.sum)

if output: self.__outputGeoTiff(mask, 'Mask', dtype = gdal.GDT_Byte)

if show: self.__showArray(mask, title = 'Mask', cmap = 'coolwarm')

return mask

def updateMask(self, filename, buffer_size = 0., classes = []):
def updateMask(self, filename, buffer_size = 0., classes = [], output = False, show = False):
"""
Function to add further pixels to mask based on a shapefile or GeoTiff file, with optional buffers.
"""
Expand All @@ -373,14 +377,18 @@ def updateMask(self, filename, buffer_size = 0., classes = []):
# Add the new raster masks to the existing mask
self.mask = np.logical_or(self.mask, mask)

if output: self.__outputGeoTiff(self.mask, 'Mask', dtype = gdal.GDT_Byte)

if show: self.__showArray(self.mask, title = 'Mask', cmap = 'coolwarm')

def resetMask(self):
"""
Function to reset a mask to the default.
"""

self.mask = self.getMask()

def getDN(self, polarisation = 'HV'):
def getDN(self, polarisation = 'HV', output = False, show = False):
"""
Loads DN (raw) values into a numpy array.
"""
Expand All @@ -405,9 +413,13 @@ def getDN(self, polarisation = 'HV'):
DN = np.zeros_like(DN_sum)
DN[self.mask == False] = (DN_sum.astype(np.float)[self.mask == False] / ((self.downsample_factor ** 2) - mask_sum[self.mask == False])).astype(np.int)

if output: self.__outputGeoTiff(DN, 'DN', dtype = gdal.GDT_Int32)

if show: self.__showArray(DN, title = 'DN', cbartitle = 'Digital Number', cmap = 'Spectral_r')

return np.ma.array(DN, mask = self.mask)

def getDOY(self, output = False):
def getDOY(self, output = False, show = False):
"""
Loads date values into a numpy array.
"""
Expand Down Expand Up @@ -441,16 +453,19 @@ def getDOY(self, output = False):
DOY = skimage.measure.block_reduce(DOY, (self.downsample_factor, self.downsample_factor), np.max)

if output: self.__outputGeoTiff(DOY, 'DOY', dtype = gdal.GDT_Int32)


if show: self.__showArray(DOY, title = 'DOY', cbartitle = 'Day', vmin = 0, vmax = 366, cmap = 'Spectral')

return DOY

def getGamma0(self, polarisation = 'HV', units = 'natural', output = False):
def getGamma0(self, polarisation = 'HV', units = 'natural', output = False, show = False):
"""
Calibrates data to gamma0 (baskscatter) in decibels or natural units.
"""

assert units == 'natural' or units == 'decibels', "Units must be 'natural' or 'decibels'. You input %s."%units

assert polarisation == 'HH' or polarisation == 'HV', "Polarisation must be 'HH' or 'HV'. You input %s."%polarisation

# Calibrate DN to units of dB
gamma0 = 10 * np.ma.log10(self.getDN(polarisation = polarisation).astype(np.float) ** 2) - 83. # units = decibels

Expand All @@ -466,9 +481,18 @@ def getGamma0(self, polarisation = 'HV', units = 'natural', output = False):

if output: self.__outputGeoTiff(gamma0, 'Gamma0')

if show:
# Different display settings depending on options
if polarisation == 'HH' and units == 'natural': vmin, vmax = 0, 0.15
if polarisation == 'HV' and units == 'natural': vmin, vmax = 0, 0.06
if polarisation == 'HH' and units == 'decibels': vmin, vmax = -15, -5
if polarisation == 'HV' and units == 'decibels': vmin, vmax = -20, -10

self.__showArray(gamma0, title = 'Gamma0 %s'%polarisation, cbartitle = units, vmin = vmin, vmax = vmax, cmap = 'Greys_r')

return gamma0

def getAGB(self, output = False):
def getAGB(self, output = False, show = False):
"""
Calibrates data to aboveground biomass (AGB).
Placeholder equation to calibrate backscatter (gamma0) to AGB (tC/ha).
Expand Down Expand Up @@ -496,9 +520,11 @@ def getAGB(self, output = False):

if output: self.__outputGeoTiff(self.AGB, 'AGB')

if show: self.__showArray(self.AGB, title = 'AGB', cbartitle = 'tC/ha', vmin = 0, vmax = 40, cmap = 'YlGn')

return self.AGB

def getWoodyCover(self, output = False):
def getWoodyCover(self, output = False, show = False):
"""
Get woody cover, based on a threshold of AGB.
min_area in ha
Expand Down Expand Up @@ -529,13 +555,15 @@ def getWoodyCover(self, output = False):
if output:

# Convert data to integer type for output
WoodyCover_out = WoodyCover.astype(np.uint8)
WoodyCover_out = self.WoodyCover.astype(np.uint8)

self.__outputGeoTiff(WoodyCover_out, 'WoodyCover', dtype = gdal.GDT_Byte)

if show: self.__showArray(self.WoodyCover, title = 'Woody Cover', cbartitle = '', vmin = 0, vmax = 1, cmap = 'summer_r')

return self.WoodyCover

def getForestPatches(self, output = False):
def getForestPatches(self, output = False, show = False):
"""
Get numbered forest patches, based on woody cover threshold.
"""
Expand All @@ -551,15 +579,17 @@ def getForestPatches(self, output = False):
min_pixels = int(round(self.area_threshold / (self.yRes * self.xRes * 0.0001)))

# Get areas that meet that threshold
_, location_id = biota.indices.getContiguousAreas(WoodyCover, True, min_pixels = min_pixels)
_, ForestPatches = biota.indices.getContiguousAreas(WoodyCover, True, min_pixels = min_pixels)

# Tidy up masked pixels
location_id.data[location_id.mask] = self.nodata
ForestPatches.data[ForestPatches.mask] = self.nodata

# Save output to class
self.ForestPatches = location_id
self.ForestPatches = ForestPatches

if output: self.__outputGeoTiff(self.ForestPatches * 1, 'ForestPatches', dtype = gdal.GDT_Int32)

if output: self.__outputGeoTiff(location_id * 1, 'ForestPatches', dtype = gdal.GDT_Int32)
if show: self.__showArray(self.ForestPatches, title = 'Forest patches', cbartitle = 'Patch ID', cmap = 'Spectral')

return self.ForestPatches

Expand All @@ -578,7 +608,23 @@ def __outputGeoTiff(self, data, output_name, dtype = 6):

# Write to disk
biota.IO.outputGeoTiff(data, filename, self.geo_t, self.proj, output_dir = self.output_dir, dtype = dtype, nodata = nodata)


def __showArray(self, data, title = '', cbartitle = '', vmin = None, vmax = None, cmap = None):
'''
'''

# Set up new figure
fig = plt.figure(figsize = (7, 6))
ax = fig.add_subplot(1, 1, 1)

# Plot map
biota.IO.buildMap(fig, ax, data, self.lat, self.lon, title = title, cbartitle = cbartitle, vmin = vmin, vmax = vmax, cmap = cmap)

# Show in window
plt.show()

# Tidy up
plt.close()

class LoadChange(object):
"""
Expand Down Expand Up @@ -677,7 +723,7 @@ def __getNodata(self, dtype = 6):

return nodata

def getAGBChange(self, output = False):
def getAGBChange(self, output = False, show = False):
'''
'''

Expand All @@ -688,12 +734,14 @@ def getAGBChange(self, output = False):

self.AGB_change = AGB_change

if output: self.__outputGeoTiff(AGB_change, 'AGBChange')
if output: self.__outputGeoTiff(self.AGB_change, 'AGBChange')

if show: self.__showArray(self.AGB_change, title = 'AGB Change', cbartitle = 'tC/ha', vmin = -10, vmax = 10, cmap = 'RdBu')

return self.AGB_change


def getChangeType(self, output = False):
def getChangeType(self, output = False, show = False):
'''
Returns pixels that meet change detection thresholds for country.
Expand Down Expand Up @@ -771,9 +819,13 @@ def getChangeType(self, output = False):
self.ChangeType = change_type
self.ChangeCode = change_code

if output:

self.__outputGeoTiff(self.changeCode, 'ChangeType', dtype = gdal.GDT_Byte)
if output: self.__outputGeoTiff(self.ChangeCode, 'ChangeType', dtype = gdal.GDT_Byte)

if show:
# Hide minor gain, minor loss and nonforest in display output
change_code_display = np.ma.array(self.ChangeCode, mask = np.zeros_like(self.ChangeCode, dtype = np.bool))
change_code_display.mask[np.logical_or(np.logical_or(change_code_display == 3, change_code_display == 4), change_code_display == 0)] = True
self.__showArray(change_code_display, title = 'Change type', cbartitle = 'Class', vmin = 1, vmax = 6, cmap = 'Spectral')

return self.ChangeType

Expand Down Expand Up @@ -802,6 +854,8 @@ def getAreaSum(self, proportion = False, output = False):
for change in totals:
totals[change] = totals[change] / change_hectares.sum()

if output: print 'TODO'

return totals


Expand All @@ -824,6 +878,8 @@ def getAGBSum(self, proportion = False, output = False):
for change_type in totals:
totals[change_type] = totals[change_type] / (self.data_t1.getAGB() * self.xRes * self.yRes * 0.0001).sum()

if output: print 'TODO'

return totals


Expand All @@ -839,5 +895,22 @@ def __outputGeoTiff(self, data, output_name, dtype = 6):

# Write to disk
biota.IO.outputGeoTiff(data, filename, self.geo_t, self.proj, output_dir = self.output_dir, dtype = dtype, nodata = nodata)

def __showArray(self, data, title = '', cbartitle = '', vmin = None, vmax = None, cmap = None):
'''
'''

# Set up new figure
fig = plt.figure(figsize = (7, 6))
ax = fig.add_subplot(1, 1, 1)

# Plot map
biota.IO.buildMap(fig, ax, data, self.lat, self.lon, title = title, cbartitle = cbartitle, vmin = vmin, vmax = vmax, cmap = cmap)

# Show in window
plt.show()

# Tidy up
plt.close()


0 comments on commit e886400

Please sign in to comment.