## LSDMap_Subplots.py
##=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## These functions are tools to deal with creating nice subplots from multiple
## rasters
##=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## FJC
## 22/12/2016
##=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
from __future__ import absolute_import, division, print_function, unicode_literals
from glob import glob
import os.path
import numpy as np
import matplotlib.pyplot as pp
import string
import matplotlib.image as mpimg
import matplotlib.cm as cmx
from matplotlib import rcParams
from . import LSDMap_GDALIO as LSDMap_IO
from . import LSDMap_BasicPlotting as LSDMap_BP
from . import labels as lsdlabels
#==============================================================================
# Convert cm to inch for figure sizing
#------------------------------------------------------------------------------
def cm2inch(value):
"""
Convert cm to inch for figure sizing
"""
return value/2.54
#==============================================================================
# Function to create nice field sites figure from all the hillshades in a folder
# N_HSFiles = number of field sites
# Also reads in a series of map files to show locations of each site. At the moment
# you have to manually indicate on these maps where the field site is - would
# be nice to automate this when I have some time to mess around with it.
# FJC 22/12/16
#------------------------------------------------------------------------------
def field_sites(DataDirectory, N_HSFiles, NRows, NCols, n_target_ticks):
# Set up fonts
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Liberation Sans']
rcParams['font.size'] = 6
# Read in the files for each site
HSFiles = sorted(glob(DataDirectory+'*_HS.bil'), key=str)
MapFiles = sorted(glob(DataDirectory+'*_map.eps'), key=str)
print(MapFiles)
n_files = len(HSFiles)+len(MapFiles)
print("Number of files = ", n_files)
# Now make the subplots
fig, ax = pp.subplots(NRows,NCols, figsize=(cm2inch(12),cm2inch(15)), frameon=False)
ax = ax.ravel()
#get a list to label the subfigures
alphabet = list(string.ascii_lowercase)
file_counter = 0
for i in range (n_files):
if i < N_HSFiles:
# first deal with the hillshade files. If there are more hillshade files simply
# change the value of N_HSFiles
print("The hillshade file name is: ", HSFiles[i])
hillshade_raster = LSDMap_IO.ReadRasterArrayBlocks(HSFiles[i])
hillshade_raster = np.ma.masked_where(hillshade_raster == -9999, hillshade_raster)
# now get the extent
extent_raster = LSDMap_IO.GetRasterExtent(HSFiles[i])
# get DEM info
CellSize,XMin,XMax,YMin,YMax = LSDMap_IO.GetUTMMaxMin(HSFiles[i])
print(YMin, YMax)
#plot the rasters
ax[i].imshow(hillshade_raster, extent = extent_raster, cmap=cmx.gray)
ax[i].text(0.05,0.97, alphabet[i], horizontalalignment='left', verticalalignment='top', bbox=dict(facecolor='white', edgecolor='k', pad=3), fontsize = 8, transform=ax[i].transAxes)
#change ticks
# xlocs = ax[i].xaxis.get_ticklocs()
# ylocs = ax[i].yaxis.get_ticklocs()
# new_xlocs,new_ylocs,new_x_labels,new_y_labels = LSDMap_BP.GetTicksForUTM(HSFiles[i],xlocs.max(),xlocs.min(),ylocs.max(),ylocs.min(),n_target_ticks)
#
# # change the location of the ticks depending on subplot placement
# if i < 2:
# ax[i].xaxis.tick_top()
# if i % 2 != 0:
# ax[i].yaxis.tick_right()
ax[i].set_xticklabels([])
ax[i].set_yticklabels([])
#ax[i].tick_params(axis='x', pad=7)
#ax[i].tick_params(axis='y', pad=7)
if i >= N_HSFiles:
print("The map file name is: ", MapFiles[file_counter])
# add in the location maps for the sites (e.g. USA, UK)
img = mpimg.imread(MapFiles[file_counter])
ax[i].imshow(img)
ax[i].text(0.02,1.00, alphabet[i], horizontalalignment='left', verticalalignment='top', fontsize = 8, transform=ax[i].transAxes)
file_counter=file_counter+1
#remove border
ax[i].axis('off')
# Save figure
# Add a big subplot to get a common x and y label for the subplots
fig.add_subplot(111, frameon=False)
# hide tick and tick label of the big axes
pp.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
#pp.xlabel('Easting (m)', fontsize=8, labelpad=-122)
#pp.ylabel('Northing (m)', fontsize=8, labelpad=10, position=(0.0,0.67))
pp.tight_layout(pad=0.1, w_pad = 0.1, h_pad = 0.2)
OutputFigureName = "field_sites"
OutputFigureFormat = 'pdf'
pp.savefig(DataDirectory+OutputFigureName + '.' + OutputFigureFormat, format=OutputFigureFormat, dpi=300)
#pp.show()
#==============================================================================
# Function to create comparison plots for floodplain mapping between the published
# flood maps and the geometric method
# FJC 22/12/16
#------------------------------------------------------------------------------
def multiple_flood_maps(DataDirectory):
"""
Make nice subplots of floodplain rasters for different field sites
"""
import seaborn as sns
# Set up fonts
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Liberation Sans']
rcParams['font.size'] = 12
# Read in the files for each site
FPFiles = sorted(glob(DataDirectory+'*_FP*.bil'), key=str)
n_files = len(FPFiles)
print("Number of files = ", n_files)
# Now make the subplots
fig, ax = pp.subplots(2,3, figsize=(cm2inch(15),cm2inch(11)))
ax = ax.ravel()
#use seaborn to get a nice color palette
cmap_oranges = sns.light_palette("#ff8f66", input="hex", as_cmap=True, reverse=True)
cmap_ice = sns.light_palette("#00ffff", input="hex", as_cmap=True, reverse=True)
alphabet = list(string.ascii_lowercase)
for i in range (n_files):
print("The floodplain file name is: ", FPFiles[i])
# get the name of the field site
fname = FPFiles[i].split('/')
print(fname)
split_fname = fname[-1].split('_')
print(split_fname)
HSFile = DataDirectory+split_fname[1]+'_'+split_fname[2]+"_HS.bil"
print(HSFile)
hillshade_raster = LSDMap_IO.ReadRasterArrayBlocks(HSFile)
FP_raster = LSDMap_IO.ReadRasterArrayBlocks(FPFiles[i])
FP_raster = np.ma.masked_where(FP_raster <= 0, FP_raster)
# now get the extent
extent_raster = LSDMap_IO.GetRasterExtent(HSFile)
# get DEM info
CellSize,XMin,XMax,YMin,YMax = LSDMap_IO.GetUTMMaxMin(HSFile)
print(YMin, YMax)
#plot the rasters
ax[i].imshow(hillshade_raster, extent = extent_raster, cmap=cmx.gray)
if i < 3:
ax[i].imshow(FP_raster, extent = extent_raster, cmap=cmap_oranges, alpha=0.8)
else:
ax[i].imshow(FP_raster, extent = extent_raster, cmap=cmap_ice, alpha=0.6)
ax[i].text(0.03,0.97, alphabet[i], bbox=dict(facecolor='white', edgecolor='k', pad=5), horizontalalignment='left', verticalalignment='top', transform=ax[i].transAxes)
#scalebars.add_scalebar(ax[i], matchx=False, sizex=500.0, loc=3, borderpad =1, lw=3, matchy=False, hidex=False, hidey=False)
ax[i].set_xticklabels([])
ax[i].set_yticklabels([])
ax[i].set_xticks([])
ax[i].set_yticks([])
ax[i].tick_params(axis='x', pad=3)
ax[i].tick_params(axis='y', pad=3)
# Save figure
pp.tight_layout(pad=0.5, h_pad=0, w_pad=0.1)
pp.subplots_adjust(wspace=0.05,hspace=0)
OutputFigureName = "Comparison_published_maps"
OutputFigureFormat = 'pdf'
pp.savefig(DataDirectory+OutputFigureName + '.' + OutputFigureFormat, format=OutputFigureFormat, transparent=True, dpi=300)
#==============================================================================
# Make subplots showing the difference between the mapped and predicted
# floodplain initiation points. Uses Fiona (yaaay) to read in the shapefiles
# FJC 05/01/17
#------------------------------------------------------------------------------
def flood_maps_with_shapefile(DataDirectory):
from fiona import collection
from descartes import PolygonPatch
# Set up fonts
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Liberation Sans']
rcParams['font.size'] = 8
# Read in the files for each site
HSFiles = sorted(glob(DataDirectory+'*_HS*.bil'), key=str)
FPFiles = sorted(glob(DataDirectory+'*_FIPs_FP*.shp'), key=str)
PointFiles = sorted(glob(DataDirectory+'*_FIPs_MP*.shp'), key=str)
n_files = len(FPFiles)
print("Number of files = ", n_files)
# Now make the subplots
fig, ax = pp.subplots(1,2, figsize=(cm2inch(12),cm2inch(7)))
ax = ax.ravel()
#get a list with the figure letterings
figure_letter = ["a", "b"]
titles = ["Mid Bailey Run, OH", "Coweeta, NC"]
for i in range (n_files):
print("The hillshade file name is: ", HSFiles[i])
print("The floodplain file name is: ", FPFiles[i])
print("The shapefile name is: ", PointFiles[i])
# get the name of the field site
split_fname = FPFiles[i].split('_FP')
print(split_fname[0])
hillshade_raster = LSDMap_IO.ReadRasterArrayBlocks(HSFiles[i])
# now get the extent
extent_raster = LSDMap_IO.GetRasterExtent(HSFiles[i])
# get DEM info
CellSize,XMin,XMax,YMin,YMax = LSDMap_IO.GetUTMMaxMin(HSFiles[i])
print(YMin, YMax)
# plot the raster
ax[i].imshow(hillshade_raster, extent = extent_raster, cmap=cmx.gray)
# plot the floodplain shapefile using fiona and descartes
with collection(FPFiles[i], 'r') as input:
for f in input:
ax[i].add_patch(PolygonPatch(f['geometry'], fc='blue', ec='blue', lw=0.1, alpha=0.8))
#plot the mapped points
with collection(PointFiles[i],'r') as input:
for point in input:
x = point['geometry']['coordinates'][0]
y = point['geometry']['coordinates'][1]
ax[i].scatter(x,y, c="red", s=15, zorder=100)
#ax[i].text(0.03,0.97, figure_letter[i], bbox=dict(facecolor='white', edgecolor='k', pad=3), horizontalalignment='left', verticalalignment='top', fontsize = 8, transform=ax[i].transAxes)
#change ticks
xlocs = ax[i].xaxis.get_ticklocs()
ylocs = ax[i].yaxis.get_ticklocs()
n_target_tics = 10
new_xlocs,new_ylocs,new_x_labels,new_y_labels = LSDMap_BP.GetTicksForUTM(HSFiles[i],xlocs.max(),xlocs.min(),ylocs.max(),ylocs.min(),n_target_tics)
ax[i].set_xticklabels(new_x_labels, rotation=30)
ax[i].set_yticklabels(new_y_labels)
#set axis limits
ax[i].set_xlim(XMin,XMax)
ax[i].set_ylim(YMin,YMax)
if i == 0:
ax[i].set_xlabel('Easting (m)', position=(1,0))
ax[i].set_ylabel('Northing (m)')
if i > 0:
ax[i].yaxis.tick_right()
ax[i].tick_params(axis='x', pad=2)
ax[i].tick_params(axis='y', pad=2)
ax[i].set_title(titles[i], fontsize=9)
# Save figure
pp.tight_layout(pad=0.5)
OutputFigureName = "Comparison_with_mapped_FIPs"
OutputFigureFormat = 'pdf'
pp.savefig(DataDirectory+OutputFigureName + '.' + OutputFigureFormat, format=OutputFigureFormat, dpi=300)
def findmaxval_multirasters(FileList):
"""
Loops through a list or array of rasters (np arrays)
and finds the maximum single value in the set of arrays.
"""
overall_max_val = 0
for i in range (len(FileList)):
raster_as_array = LSDMap_IO.ReadRasterArrayBlocks(FileList[i])
this_max_val = np.nanmax(raster_as_array)
if this_max_val > overall_max_val:
overall_max_val = this_max_val
print(overall_max_val)
return overall_max_val
def findminval_multirasters(FileList):
"""
Loops through a list or array of rasters (np arrays)
and finds the minimum single value in the set of arrays.
"""
overall_min_val = 0
for i in range (len(FileList)):
raster_as_array = LSDMap_IO.ReadRasterArrayBlocks(FileList[i])
this_min_val = np.nanmin(raster_as_array)
if this_min_val > overall_min_val:
overall_min_val = this_min_val
print(overall_min_val)
return overall_min_val
def MultiDrapeFloodMaps(DataDir, ElevationRaster, DrapeRasterWild, cmap,
drape_min_threshold=None, drape_max=None, cbar_label=None):
"""Creates a figure with multiple drape maps over a hillshade.
Plots flood extents from water depth rasters
draped over the catchment elevation raster
in a series of subplots
Takes a wildcard for the drapes
Expexts a fixed elevation raster, but this could be
modified in future.
Parameters:
DataDir (str): Path to the directory containing the data files
ElevationRaster (str): Name of the elevation raster used to create the hillshade
DrapeRasterWild (str): Wildcard string used to find all the drape files in the
directory.
cmap: Can be the string name of a colourmap, or a Colourmap object
drape_min (float, optional): Minimum value for the drape raster, i.e. values
below this threshold will be masked and not plotted.
drape_max (float, optional): Maximum value for the drape raster, i.e. values
above this value will be masked and not plotted.
cbar_label (str, optional): Label for the colourbar on the figure. This
is the colourbar for the drape colourmap.
Notes:
Consider, if plotting multiple datasets, how you
are going to deal with min a max values in the colur range.
imshow will automatically set vmin and vmax and stretch the colour bar
over this - which can be visually misleading. Ideally, you
want to have the same colour map used for *all* subplots, and
this is not default behaviour.
Note: If `drape_max` is not set, the function searches for the maximum value
in the range of rasters found by expanding the `DrapeRasterWild` argument
and searching for the maximum value out of all rasters found.
Raises:
Exception: If the maximum value in the drape maps could not be found.
"""
f, ax_arr = pp.subplots(2, 2, figsize=(10, 5), sharex=True, sharey=True)
ax_arr = ax_arr.ravel()
FPFiles = sorted(glob(DataDir+DrapeRasterWild), key=str)
n_files = len(FPFiles)
print("Number of files = ", n_files)
elev_raster_file = DataDir + ElevationRaster
hillshade = LSDMap_BP.Hillshade(elev_raster_file)
#hillshade_array = LSDP.ReadRasterArrayBlocks(elev_raster_file)
# now get the extent
extent_raster = LSDMap_IO.GetRasterExtent(elev_raster_file)
x_min = extent_raster[0]
x_max = extent_raster[1]
y_min = extent_raster[2]
y_max = extent_raster[3]
# now get the tick marks
n_target_tics = 5
xlocs,ylocs,new_x_labels,new_y_labels = LSDMap_BP.GetTicksForUTM(elev_raster_file,x_max,x_min,y_max,y_min,n_target_tics)
print("xmax: " + str(x_max))
print("xmin: " + str(x_min))
print("ymax: " + str(y_max))
print("ymin: " + str(y_min))
"""
Find the maximum water depth in all rasters.
You need this to normalize the colourscale accross
all plots when teh imshow is done later.
"""
try:
print("Calculating max drape raster value by scanning rasters...")
max_water_depth = findmaxval_multirasters(FPFiles)
drape_max = max_water_depth
except:
print("Something went wrong trying to obtain the max value in \
your drape raster file list.")
finally:
print("The drape(s) max value is set to: ", drape_max)
#im = mpimg.AxesImage()
for i in range(n_files):
print("The floodplain file name is: ", FPFiles[i])
FP_raster = LSDMap_IO.ReadRasterArrayBlocks(FPFiles[i])
#FP_raster = np.ma.masked_where(FP_raster <= 0, FP_raster)
filename = os.path.basename(FPFiles[i])
title = lsdlabels.make_line_label(filename)
print(title)
low_values_index = FP_raster < drape_min_threshold
FP_raster[low_values_index] = np.nan
im = ax_arr[i].imshow(hillshade, "gray", extent=extent_raster, interpolation="nearest")
"""
Now we can set vmax to be the maximum water depth we calcualted earlier, making our separate
subplots all have the same colourscale
"""
im = ax_arr[i].imshow(FP_raster, cmap, extent=extent_raster,
alpha=1.0, interpolation="nearest",
vmin=drape_min_threshold,
vmax=drape_max)
ax_arr[i].set_title(title)
pp.setp( ax_arr[i].xaxis.get_majorticklabels(), rotation=70 )
f.subplots_adjust(right=0.85)
cax = f.add_axes([0.9, 0.1, 0.03, 0.8])
cbar = f.colorbar(im, cax=cax)
cbar.set_label(cbar_label)
#cbar.set_ticks(np.linspace(0, 8, 8))
#cbar = LSDP.colours.colorbar_index(f, cax, 8, cmap,
# drape_min_threshold, drape_max)
#tick_locator = ticker.MaxNLocator(nbins=8)
#cbar.locator = tick_locator
#cbar.update_ticks()
f.text(0.5, 0.04, 'Easting (m)', ha='center', fontsize=17)
f.text(0.04, 0.5, 'Northing (m)', va='center', rotation='vertical', fontsize=17)
def MultiDrapeErodeDiffMaps(DataDir, ElevationRaster, DrapeRasterWild, cmap,
drape_min_threshold=None, cbar_label=None,
drape_max_threshold=None,
middle_mask_range=None):
"""Plots multiple drape maps of erosion/deposition (a DEM of difference)
over a hillshade raster of the basin.
Takes a wildcard for the drapes
Expexts a single elevation raster for the background hillshade,
but this could be modified in future.
Parameters:
DataDir (str): Path to the directory containing the data files
ElevationRaster (str): Name of the elevation raster used to create the hillshade
DrapeRasterWild (str): Wildcard string used to find all the drape files in the
directory.
cmap: Can be the string name of a colourmap, or a Colourmap object
drape_min_threshold (float, optional): Minimum value for the drape raster, i.e. values
below this threshold will be masked and not plotted.
drape_max_threshold (float, optional): Maximum value for the drape raster, i.e. values
above this value will be masked and not plotted.
cbar_label (str, optional): Label for the colourbar on the figure. This
is the colourbar for the drape colourmap.
middle_mask_range (tuple, optional): A tuple or list of two values, used
to mask an inner range of values in the drape raster.
e.g. if you pass `(-0.1, 0.1)` then all the values
in the range -0.1 to 0.1 will be masked and not plotted
on the final map. Use for masking very small values
either side of zero.
Notes:
Consider, if plotting multiple datasets, how you
are going to deal with min a max values in the colur range.
imshow will automatically set vmin and vmax and stretch the colour bar
over this - which can be visually misleading. Ideally, you
want to have the same colour map used for *all* subplots, and
this is not default behaviour.
Note: If `drape_max_threshold` is not set, the function searches for the maximum value
in the range of rasters found by expanding the `DrapeRasterWild` argument
and searching for the maximum value out of all rasters found.
Raises:
Exception: If the maximum value in the drape maps could not be found.
Author: DAV & FJC
"""
#import lsdmatplotlibextensions as mplext
f, ax_arr = pp.subplots(2, 2, figsize=(10, 5), sharex=True, sharey=True)
ax_arr = ax_arr.ravel()
FPFiles = sorted(glob(DataDir+DrapeRasterWild), key=str)
n_files = len(FPFiles)
print("Number of files = ", n_files)
elev_raster_file = DataDir + ElevationRaster
hillshade = LSDMap_BP.Hillshade(elev_raster_file)
#hillshade_array = LSDP.ReadRasterArrayBlocks(elev_raster_file)
# now get the extent
extent_raster = LSDMap_IO.GetRasterExtent(elev_raster_file)
x_min = extent_raster[0]
x_max = extent_raster[1]
y_min = extent_raster[2]
y_max = extent_raster[3]
# now get the tick marks
n_target_tics = 5
xlocs, ylocs, new_x_labels, new_y_labels = LSDMap_BP.GetTicksForUTM(
elev_raster_file,
x_max,
x_min,
y_max,
y_min,
n_target_tics)
print("xmax: " + str(x_max))
print("xmin: " + str(x_min))
print("ymax: " + str(y_max))
print("ymin: " + str(y_min))
"""
Find the maximum water depth in all rasters.
You need this to normalize the colourscale accross
all plots when teh imshow is done later.
"""
if drape_max_threshold is None:
try:
print("Calculating max drape raster value by scanning rasters...")
max_water_depth = findmaxval_multirasters(FPFiles)
drape_max_threshold = max_water_depth
except ValueError:
print("Something went wrong trying to obtain the max value in \
your drape raster file list.")
finally:
print("The drape(s) max value is set to: ", drape_max_threshold)
if drape_min_threshold is None:
try:
print("Calculating min drape raster value by scanning rasters...")
min_water_depth = findminval_multirasters(FPFiles)
drape_min_threshold = min_water_depth
except ValueError:
print("Something went wrong trying to obtain the min value in \
your drape raster file list.")
finally:
print("The drape(s) min value is set to: ", drape_min_threshold)
for i in range(n_files):
print("The floodplain file name is: ", FPFiles[i])
FP_raster = LSDMap_IO.ReadRasterArrayBlocks(FPFiles[i])
#FP_raster = np.ma.masked_where(FP_raster <= 0, FP_raster)
filename = os.path.basename(FPFiles[i])
title = lsdlabels.make_line_label(filename)
print(title)
# Mask the extreme high values
hi_values_index = FP_raster > drape_max_threshold
FP_raster[hi_values_index] = np.nan
# Mask the extreme low values
lo_values_index = FP_raster < drape_min_threshold
FP_raster[lo_values_index] = np.nan
# Mask the middle values that are really close to zero (i.e. if you
# have negative and positive values in the raster, such as in a DEM
# of difference with both erosion and deposition.)
if middle_mask_range is not None:
masked_mid_values_index = (np.logical_and(FP_raster > middle_mask_range[0],
FP_raster < middle_mask_range[1]))
FP_raster[masked_mid_values_index] = np.nan
im = ax_arr[i].imshow(hillshade, "gray", extent=extent_raster, interpolation="nearest")
"""
Now we can set vmax to be the maximum water depth we calcualted earlier, making our separate
subplots all have the same colourscale
"""
im = ax_arr[i].imshow(FP_raster, cmap, extent=extent_raster,
alpha=1.0, interpolation="nearest",
vmin=drape_min_threshold,
vmax=drape_max_threshold)
ax_arr[i].set_title(title)
pp.setp( ax_arr[i].xaxis.get_majorticklabels(), rotation=70 )
f.subplots_adjust(right=0.85)
cax = f.add_axes([0.9, 0.1, 0.03, 0.8])
cbar = f.colorbar(im, cax=cax)
cbar.set_label(cbar_label)
#cbar.set_ticks(np.linspace(0, 8, 8))
#cbar = mplext.colours.colorbar_index(f, cax, 8, cmap,
# drape_min_threshold, drape_max)
#tick_locator = ticker.MaxNLocator(nbins=8)
#cbar.locator = tick_locator
#cbar.update_ticks()
f.text(0.5, 0.04, 'Easting (m)', ha='center', fontsize=17)
f.text(0.04, 0.5, 'Northing (m)', va='center', rotation='vertical', fontsize=17)
