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