Source code for shakelib.plotting.contour
# usgs imports
from impactutils.colors.cpalette import ColorPalette
# third party imports
from shapely.geometry import MultiLineString
from shapely.geometry import mapping
from scipy.ndimage.filters import median_filter
from skimage import measure
import numpy as np
from openquake.hazardlib import imt
[docs]def contour(imtdict, imtype, filter_size, gmice):
"""
Generate contours of a specific IMT and return as a Shapely
MultiLineString object.
Args:
container (ShakeMapOutputContainer): ShakeMapOutputContainer
with ShakeMap output data.
imtype (str): String containing the name of an Intensity
Measure Type found in container.
filter_size (int): Integer filter (see
https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.ndimage.filters.median_filter.html)
Returns:
list: List of dictionaries containing two fields
- geometry: GeoJSON-like representation of one of the objects
in https://toblerity.org/fiona/manual.html#geometry-types
- properties: Dictionary of properties describing that
feature.
Raises:
NotImplementedError -- if the user attempts to contour a data file
with sets of points rather than grids.
""" # noqa
oqimt = imt.from_string(imtype)
intensity_colormap = ColorPalette.fromPreset("mmi")
grid = imtdict["mean"]
metadata = imtdict["mean_metadata"]
if imtype == "MMI":
sgrid = grid
units = "mmi"
elif imtype == "PGV":
sgrid = np.exp(grid)
units = "cms"
else:
sgrid = np.exp(grid) * 100.0
units = "pctg"
if filter_size > 0:
fgrid = median_filter(sgrid, size=int(filter_size))
else:
fgrid = sgrid
interval_type = "log"
if imtype == "MMI":
interval_type = "linear"
if np.all(np.isnan(fgrid)): # data is totally empty; no contours
intervals = np.array([])
else:
grid_min = np.nanmin(fgrid)
grid_max = np.nanmax(fgrid)
if grid_max - grid_min:
intervals = getContourLevels(grid_min, grid_max, itype=interval_type)
else: # data is totally flat; don't draw any contours
intervals = np.array([])
lonstart = metadata["xmin"]
latstart = metadata["ymin"]
lonend = metadata["xmax"]
if lonend < lonstart:
lonstart -= 360
lonspan = np.abs(lonend - lonstart)
latspan = np.abs(metadata["ymax"] - latstart)
nlon = metadata["nx"]
nlat = metadata["ny"]
line_strings = [] # dictionary of MultiLineStrings and props
for cval in intervals:
contours = measure.find_contours(fgrid, cval)
#
# Convert coords to geographic coordinates; the coordinates
# are returned in row, column order (i.e., (y, x))
#
new_contours = []
plot_contours = []
for ic, coords in enumerate(contours): # coords is a line segment
#
# This greatly reduces the number of points in the contours
# without changing their shape too much
#
coords = measure.approximate_polygon(coords, filter_size / 20)
mylons = np.around(coords[:, 1] * lonspan / nlon + lonstart, decimals=6)
mylats = np.around(
(nlat - coords[:, 0]) * latspan / nlat + latstart, decimals=6
)
contours[ic] = np.hstack(
(mylons[:].reshape((-1, 1)), mylats[:].reshape((-1, 1)))
)
plot_contours.append(contours[ic])
new_contours.append(contours[ic].tolist())
if len(new_contours):
mls = MultiLineString(new_contours)
props = {"value": cval, "units": units}
if imtype == "MMI":
pass
elif imtype == "PGV":
lcval = np.log(cval)
else:
lcval = np.log(cval / 100)
if gmice:
mmival = gmice.getMIfromGM(np.array([lcval]), oqimt)[0][0]
elif imtype == "MMI":
mmival = cval
else:
mmival = 1
color_array = np.array(intensity_colormap.getDataColor(mmival))
color_rgb = np.array(color_array[0:3] * 255, dtype=int).tolist()
props["color"] = "#%02x%02x%02x" % tuple(color_rgb)
if imtype == "MMI":
if (cval * 2) % 2 == 1:
props["weight"] = 4
else:
props["weight"] = 2
else:
props["weight"] = 4
line_strings.append({"geometry": mapping(mls), "properties": props})
return line_strings
[docs]def getContourLevels(dmin, dmax, itype="log"):
"""
Get contour levels given min/max values and desired IMT.
Use itype='log' for any IMT that is logarithmically distributed, such as
PGA, PGV, and Sa. Linear for MMI.
Args:
dmin (float): Minimum value of data to contour.
dmax (float): Maximum value of data to contour.
itype (str): Interval type; default is 'log', anythign else
indicates linear intervals.
Returns:
ndarray: Numpy array of contour levels.
"""
if itype == "log":
# Within-decade label values
dec_inc = np.array([1, 2, 5], dtype=float)
# Upper and lower decades
lower_dec = np.floor(np.log10(dmin))
upper_dec = np.ceil(np.log10(dmax))
# Don't make a crazy number of contours if there are very small
# values in the input
if lower_dec < upper_dec - 4:
lower_dec = upper_dec - 4
# Array of decades
decades = np.arange(lower_dec, upper_dec + 1)
# Construct levels
levels = np.concatenate([np.power(10, d) * dec_inc for d in decades])
levels = levels[(levels < dmax) & (levels > dmin)]
if np.size(levels) == 0:
levels = np.array([(dmin + dmax) / 2])
else:
# MMI contours are every 0.5 units
levels = np.arange(np.ceil(dmin * 2) / 2, np.floor(dmax * 2) / 2 + 0.5, 0.5)
return levels