Source code for csep.core.regions

# Python imports
import itertools
import os
from itertools import compress
from xml.etree import ElementTree as ET

# Third-party imports
import numpy
import numpy as np
import mercantile
from shapely import geometry
from shapely.ops import unary_union

# PyCSEP imports
from csep.utils.calc import bin1d_vec, cleaner_range, first_nonnan, last_nonnan
from csep.utils.scaling_relationships import WellsAndCoppersmith

from csep.models import Polygon

def california_relm_collection_region(dh_scale=1, magnitudes=None, name="relm-california-collection", use_midpoint=True):
    """ Return collection region for California RELM testing region

    Args:
        dh_scale (int): factor of two multiple to change the grid size
        mangitudes (array-like): array representing the lower bin edges of the magnitude bins
        name (str): human readable identifer
        use_midpoints (bool): if true, treat values in file as midpoints. default = true.

    Returns:
        :class:`csep.core.spatial.CartesianGrid2D`

    Raises:
        ValueError: dh_scale must be a factor of two

    """
    if dh_scale % 2 != 0 and dh_scale != 1:
        raise ValueError("dh_scale must be a factor of two or dh_scale must equal unity.")

    # we can hard-code the dh because we hard-code the filename
    dh = 0.1
    root_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
    filepath = os.path.join(root_dir, 'artifacts', 'Regions', 'RELMCollectionArea.dat')
    points = numpy.loadtxt(filepath)
    if use_midpoint:
        origins = numpy.array(points) - dh / 2
    else:
        origins = numpy.array(points)

    if dh_scale > 1:
        origins = increase_grid_resolution(origins, dh, dh_scale)
        dh = dh / dh_scale

    # turn points into polygons and make region object
    bboxes = compute_vertices(origins, dh)
    relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name)

    if magnitudes is not None:
        relm_region.magnitudes = magnitudes

    return relm_region

[docs]def california_relm_region(dh_scale=1, magnitudes=None, name="relm-california", use_midpoint=True): """ Returns class representing California testing region. This region can be used to create gridded datasets for earthquake forecasts. The XML file appears to use the midpoint, and the .dat file uses the origin in the "lower left" corner. Args: dh_scale: can resample this grid by factors of 2 Returns: :class:`csep.core.spatial.CartesianGrid2D` Raises: ValueError: dh_scale must be a factor of two """ if dh_scale % 2 != 0 and dh_scale != 1: raise ValueError("dh_scale must be a factor of two or dh_scale must equal unity.") # use default file path from python package root_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) filepath = os.path.join(root_dir, 'artifacts', 'Regions', 'csep-forecast-template-M5.xml') csep_template = os.path.expanduser(filepath) points, dh = parse_csep_template(csep_template) if use_midpoint: origins = numpy.array(points) - dh / 2 else: origins = numpy.array(points) if dh_scale > 1: origins = increase_grid_resolution(origins, dh, dh_scale) dh = dh / dh_scale # turn points into polygons and make region object bboxes = compute_vertices(origins, dh) relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name) if magnitudes is not None: relm_region.magnitudes = magnitudes return relm_region
[docs]def italy_csep_region(dh_scale=1, magnitudes=None, name="csep-italy", use_midpoint=True): """ Returns class representing Italian testing region. This region can be used to create gridded datasets for earthquake forecasts. The region is defined by the file 'forecast.italy.M5.xml' and contains a spatially gridded region with 0.1° x 0.1° cells. Args: dh_scale: can resample this grid by factors of 2 magnitudes (array-like): bin edges for magnitudes. if provided, will be bound to the output region class. this argument provides a short-cut for creating space-magnitude regions. name (str): human readable identify given to the region use_midpoint (bool): if true, treat values in file as midpoints. default = true. Returns: :class:`csep.core.spatial.CartesianGrid2D` Raises: ValueError: dh_scale must be a factor of two """ if dh_scale % 2 != 0 and dh_scale != 1: raise ValueError("dh_scale must be a factor of two or dh_scale must equal unity.") # use default file path from python package root_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) filepath = os.path.join(root_dir, 'artifacts', 'Regions', 'forecast.italy.M5.xml') csep_template = os.path.expanduser(filepath) points, dh = parse_csep_template(csep_template) if use_midpoint: origins = numpy.array(points) - dh / 2 else: origins = numpy.array(points) if dh_scale > 1: origins = increase_grid_resolution(origins, dh, dh_scale) dh = dh / dh_scale # turn points into polygons and make region object bboxes = compute_vertices(origins, dh) italy_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name) if magnitudes is not None: italy_region.magnitudes = magnitudes return italy_region
def italy_csep_collection_region(dh_scale=1, magnitudes=None, name="csep-italy-collection", use_midpoint=True): """ Return collection region for Italy CSEP collection region Args: dh_scale (int): factor of two multiple to change the grid size mangitudes (array-like): array representing the lower bin edges of the magnitude bins name (str): human readable identifer use_midpoint (bool): if true, treat values in file as midpoints. default = true. Returns: :class:`csep.core.spatial.CartesianGrid2D` Raises: ValueError: dh_scale must be a factor of two """ if dh_scale % 2 != 0 and dh_scale != 1: raise ValueError("dh_scale must be a factor of two or dh_scale must equal unity.") # we can hard-code the dh because we hard-code the filename dh = 0.1 root_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) filepath = os.path.join(root_dir, 'artifacts', 'Regions', 'italy.collection.nodes.dat') points = numpy.loadtxt(filepath) if use_midpoint: origins = numpy.array(points) - dh / 2 else: origins = numpy.array(points) if dh_scale > 1: origins = increase_grid_resolution(origins, dh, dh_scale) dh = dh / dh_scale # turn points into polygons and make region object bboxes = compute_vertices(origins, dh) relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name) if magnitudes is not None: relm_region.magnitudes = magnitudes return relm_region def nz_csep_region(dh_scale=1, magnitudes=None, name="csep-nz", use_midpoint=True): """ Return collection region for the New Zealand CSEP testing region Args: dh_scale (int): factor of two multiple to change the grid size mangitudes (array-like): array representing the lower bin edges of the magnitude bins name (str): human readable identifer use_midpoints (bool): if true, treat values in file as midpoints. default = true. Returns: :class:`csep.core.spatial.CartesianGrid2D` Raises: ValueError: dh_scale must be a factor of two """ if dh_scale % 2 != 0 and dh_scale != 1: raise ValueError("dh_scale must be a factor of two or dh_scale must equal unity.") # we can hard-code the dh because we hard-code the filename dh = 0.1 root_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) filepath = os.path.join(root_dir, 'artifacts', 'Regions', 'nz.testing.nodes.dat') points = numpy.loadtxt(filepath) if use_midpoint: origins = numpy.array(points) - dh / 2 else: origins = numpy.array(points) if dh_scale > 1: origins = increase_grid_resolution(origins, dh, dh_scale) dh = dh / dh_scale # turn points into polygons and make region object bboxes = compute_vertices(origins, dh) nz_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name) if magnitudes is not None: nz_region.magnitudes = magnitudes return nz_region def nz_csep_collection_region(dh_scale=1, magnitudes=None, name="csep-nz-collection", use_midpoint=True): """ Return collection region for the New Zealand CSEP collection region Args: dh_scale (int): factor of two multiple to change the grid size mangitudes (array-like): array representing the lower bin edges of the magnitude bins name (str): human readable identifer use_midpoints (bool): if true, treat values in file as midpoints. default = true. Returns: :class:`csep.core.spatial.CartesianGrid2D` Raises: ValueError: dh_scale must be a factor of two """ if dh_scale % 2 != 0 and dh_scale != 1: raise ValueError("dh_scale must be a factor of two or dh_scale must equal unity.") # we can hard-code the dh because we hard-code the filename dh = 0.1 root_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) filepath = os.path.join(root_dir, 'artifacts', 'Regions', 'nz.collection.nodes.dat') points = numpy.loadtxt(filepath) if use_midpoint: origins = numpy.array(points) - dh / 2 else: origins = numpy.array(points) if dh_scale > 1: origins = increase_grid_resolution(origins, dh, dh_scale) dh = dh / dh_scale # turn points into polygons and make region object bboxes = compute_vertices(origins, dh) nz_collection_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name) if magnitudes is not None: nz_collection_region.magnitudes = magnitudes return nz_collection_region
[docs]def global_region(dh=0.1, name="global", magnitudes=None): """ Creates a global region used for evaluating gridded forecasts on the global scale. The gridded region corresponds to the Args: dh: Returns: csep.utils.CartesianGrid2D: """ # generate latitudes lons = cleaner_range(-180.0, 179.9, dh) lats = cleaner_range(-90, 89.9, dh) coords = itertools.product(lons,lats) region = CartesianGrid2D([Polygon(bbox) for bbox in compute_vertices(coords, dh)], dh, name=name) if magnitudes is not None: region.magnitudes = magnitudes return region
[docs]def magnitude_bins(start_magnitude, end_magnitude, dmw): """ Returns array holding magnitude bin edges. The output from this function is monotonically increasing and equally spaced bin edges that can represent magnitude bins. Args: start_magnitude (float) end_magnitude (float) dmw (float): magnitude spacing Returns: bin_edges (numpy.ndarray) """ # convert to integers to prevent accumulating floating point errors const = 10000 start = numpy.floor(const * start_magnitude) end = numpy.floor(const * end_magnitude) d = const * dmw return numpy.arange(start, end + d / 2, d) / const
[docs]def create_space_magnitude_region(region, magnitudes): """Simple wrapper to create space-magnitude region """ if not (isinstance(region, CartesianGrid2D) or isinstance(region, QuadtreeGrid2D)) : raise TypeError("region must be CartesianGrid2D") # bind to region class if magnitudes is None: raise ValueError("magnitudes should not be None if creating space-magnitude region.") region.magnitudes = magnitudes region.num_mag_bins = len(region.magnitudes) return region
[docs]def parse_csep_template(xml_filename): """ Reads CSEP XML template file and returns the lat/lon values for the forecast. Returns: list of tuples where tuple is (lon, lat) """ tree = ET.parse(xml_filename) root = tree.getroot() points = [] for cell in root.iter('{http://www.scec.org/xml-ns/csep/forecast/0.1}cell'): points.append((float(cell.attrib['lon']), float(cell.attrib['lat']))) # get cell spacing data = root.find('{http://www.scec.org/xml-ns/csep/forecast/0.1}forecastData') dh_elem = data.find('{http://www.scec.org/xml-ns/csep/forecast/0.1}defaultCellDimension') dh_lat = float(dh_elem.attrib['latRange']) dh_lon = float(dh_elem.attrib['lonRange']) if not numpy.isclose(dh_lat, dh_lon): raise ValueError("dh_lat must equal dh_lon. grid needs to be regular.") return points, dh_lat
[docs]def increase_grid_resolution(points, dh, factor): """ Takes a set of origin points and returns a new set with higher grid resolution. assumes the origin point is in the lower left corner. the new dh is dh / factor. This implementation requires that the decimation factor be a multiple of 2. Args: points: list of (lon,lat) tuples dh: old grid spacing factor: amount to reduce Returns: points: list of (lon,lat) tuples with spacing dh / scale """ # short-circuit recursion if factor == 1: return points # handle edge cases assert factor % 2 == 0 assert factor >= 1 # first start out new_points = set() new_dh = dh / 2 for point in points: bbox = compute_vertex(point, new_dh) for pnt in bbox: new_points.add(pnt) # call function again with new_points, new_dh, new_factor new_factor = factor / 2 return increase_grid_resolution(list(new_points), new_dh, new_factor)
[docs]def masked_region(region, polygon): """ Build a new region based off the coordinates in the polygon. Args: region: CartesianGrid2D object polygon: Polygon object Returns: new_region: CartesianGrid2D object """ # contains is true if spatial cell in region is inside the polygon contains = polygon.contains(region.midpoints()) # compress only returns elements that are true, effectively removing elements outside of the polygons new_polygons = list(compress(region.polygons, contains)) # create new region with the spatial cells inside the polygon return CartesianGrid2D(new_polygons, region.dh)
[docs]def generate_aftershock_region(mainshock_mw, mainshock_lon, mainshock_lat, num_radii=3, region=california_relm_region, **kwargs): """ Creates a spatial region around a given epicenter The method uses the Wells and Coppersmith scaling relationship to determine the average fault length and creates a circular region centered at (mainshock_lon, mainshock_lat) with radius = num_radii. Args: mainshock_mw (float): magnitude of mainshock mainshock_lon (float): epicentral longitude mainshock_lat (float): epicentral latitude num_radii (float/int): number of radii of circular region region (callable): returns :class:`csep.utils.spatial.CartesianGrid2D` **kwargs (dict): passed to region callable Returns: :class:`csep.utils.spatial.CartesianGrid2D` """ rupture_length = WellsAndCoppersmith.mag_length_strike_slip(mainshock_mw) * 1000 aftershock_polygon = Polygon.from_great_circle_radius((mainshock_lon, mainshock_lat), num_radii * rupture_length, num_points=100) aftershock_region = masked_region(region(**kwargs), aftershock_polygon) return aftershock_region
def grid_spacing(vertices): """ Figures out the length and Args: vertices: Vertices describe a single node in grid. Returns: dh: grid spacing Raises: ValueError """ # get first two vertices a = vertices[0] b = vertices[1] # compute both differences, because unless point is the same one is bound to be the dh d1 = numpy.abs(b[0] - a[0]) d2 = numpy.abs(b[1] - a[1]) if not numpy.allclose(d1, d2): raise ValueError("grid spacing must be regular for cartesian grid.") dh = numpy.max([d1, d2]) # this would happen if the same point is repeated twice if dh == 0: raise ValueError("Problem computing grid spacing cannot be zero.") return dh def compute_vertex(origin_point, dh, tol=numpy.finfo(float).eps): """ Computes the bounding box of a rectangular polygon given its origin points and spacing dh. Args: origin_points: list of tuples, where tuple is (x, y) dh: spacing tol: used to eliminate overlapping polygons in the case of a rectangular mesh, defaults to the machine tolerance. Returns: list of polygon edges """ bbox = ((origin_point[0], origin_point[1]), (origin_point[0], origin_point[1] + dh - tol), (origin_point[0] + dh - tol, origin_point[1] + dh - tol), (origin_point[0] + dh - tol, origin_point[1])) return bbox def compute_vertices(origin_points, dh, tol=numpy.finfo(float).eps): """ Wrapper function to compute vertices for multiple points. Default tolerance is set to machine precision of floating point number. Args: origin_points: 2d ndarray Notes: (x,y) should be accessible like: #>>> x_coords = origin_points[:,0] #>>> y_coords = origin_points[:,1] """ return list(map(lambda x: compute_vertex(x, dh, tol=tol), origin_points)) def _bin_catalog_spatio_magnitude_counts(lons, lats, mags, n_poly, mask, idx_map, binx, biny, mag_bins, tol=0.00001): """ Returns a list of event counts as ndarray with shape (n_poly, n_cat) where each value represents the event counts within the polygon. Using [:, :, 1] index of the mask, we store the mapping between the index of n_poly and that polygon in the mask. Additionally, the polygons are ordered such that the index of n_poly in the result corresponds to the index of the polygons. Eventually, we can make a structure that could contain both of these, but the trade-offs will need to be compared against performance. """ # index in cartesian grid for events in data. note, this has a different index than the # vector of polygons. this mapping is stored in [:,:,1] index of mask # index in 2d grid idx = bin1d_vec(lons, binx) idy = bin1d_vec(lats, biny) mag_idxs = bin1d_vec(mags, mag_bins, tol=tol, right_continuous=True) # start with zero event counts in each bin event_counts = numpy.zeros((n_poly, len(mag_bins))) # does not seem that we can vectorize this part skipped = [] for i in range(idx.shape[0]): if not mask[idy[i], idx[i]] and idy[i] != -1 and idx[i] != -1 and mag_idxs[i] != -1: # getting spatial bin from mask hash_idx = int(idx_map[idy[i], idx[i]]) mag_idx = mag_idxs[i] # update event counts in that polygon event_counts[(hash_idx, mag_idx)] += 1 else: skipped.append((lons[i], lats[i], mags[i])) return event_counts, skipped def _bin_catalog_spatial_counts(lons, lats, n_poly, mask, idx_map, binx, biny): """ Returns a list of event counts as ndarray with shape (n_poly) where each value represents the event counts within the polygon. Using [:, :, 1] index of the mask, we store the mapping between the index of n_poly and that polygon in the mask. Additionally, the polygons are ordered such that the index of n_poly in the result corresponds to the index of the polygons. We can make a structure that could contain both of these, but the trade-offs will need to be compared against performance. """ ai, bi = binx, biny # will return negative idx = bin1d_vec(lons, ai) idy = bin1d_vec(lats, bi) # bin1d returns -1 if outside the region # todo: think about how to change this behavior for less confusions, bc -1 is an actual value that can be chosen bad = (idx == -1) | (idy == -1) | (mask[idy,idx] == 1) # this can be memory optimized by keeping short list and storing index, only for case where n/2 events event_counts = numpy.zeros(n_poly) # selecting the indexes into polygons correspoding to lons and lats within the grid hash_idx = idx_map[idy[~bad],idx[~bad]].astype(int) # aggregate in counts numpy.add.at(event_counts, hash_idx, 1) return event_counts def _bin_catalog_probability(lons, lats, n_poly, mask, idx_map, binx, biny): """ Returns a list of event counts as ndarray with shape (n_poly) where each value represents the event counts within the polygon. Using [:, :, 1] index of the mask, we store the mapping between the index of n_poly and that polygon in the mask. Additionally, the polygons are ordered such that the index of n_poly in the result corresponds to the index of the polygons. We can make a structure that could contain both of these, but the trade-offs will need to be compared against performance. """ ai, bi = binx, biny # returns -1 if outside of the bbox idx = bin1d_vec(lons, ai) idy = bin1d_vec(lats, bi) bad = (idx == -1) | (idy == -1) | (mask[idy, idx] == 1) event_counts = numpy.zeros(n_poly) # [:,:,1] is a mapping from the polygon array to cartesian grid hash_idx = idx_map[idy[~bad],idx[~bad]].astype(int) # dont accumulate just set to one for probability event_counts[hash_idx] = 1 return event_counts
[docs]class CartesianGrid2D: """Represents a 2D cartesian gridded region. The class provides functions to query onto an index 2D Cartesian grid and maintains a mapping between space coordinates defined by polygons and the index into the polygon array. Custom regions can be easily created by using the from_polygon classmethod. This function will accept an arbitrary closed polygon and return a CartesianGrid class with only points inside the polygon to be valid. """
[docs] def __init__(self, polygons, dh, name='cartesian2d', mask=None): self.polygons = polygons self.poly_mask = mask self.dh = dh self.name = name a, xs, ys = self._build_bitmask_vec() # in mask, True = bad value and False = good value self.bbox_mask = a[:,:,0] # contains the mapping from polygon_index to the mask self.idx_map = a[:,:,1] # index values of polygons array into the 2d cartesian grid, based on the midpoint. self.xs = xs self.ys = ys # Bounds [origin, top_right] orgs = self.origins() self.bounds = numpy.column_stack((orgs, orgs + dh))
def __eq__(self, other): return self.to_dict() == other.to_dict() @property def num_nodes(self): """ Number of polygons in region """ return len(self.polygons) def get_index_of(self, lons, lats): """ Returns the index of lons, lats in self.polygons Args: lons: ndarray-like lats: ndarray-like Returns: idx: ndarray-like """ idx = bin1d_vec(numpy.array(lons), self.xs) idy = bin1d_vec(numpy.array(lats), self.ys) if numpy.any(idx == -1) or numpy.any(idy == -1): raise ValueError("at least one lon and lat pair contain values that are outside of the valid region.") if numpy.any(self.bbox_mask[idy, idx] == 1): raise ValueError("at least one lon and lat pair contain values that are outside of the valid region.") return self.idx_map[idy,idx].astype(numpy.int) def get_location_of(self, indices): """ Returns the polygon associated with the index idx. Args: idx: index of polygon in region Returns: Polygon """ indices = list(indices) polys = [self.polygons[idx] for idx in indices] return polys def get_masked(self, lons, lats): """Returns bool array lons and lats are not included in the spatial region. .. note:: The ordering of lons and lats should correspond to the ordering of the lons and lats in the data. Args: lons: array-like lats: array-like Returns: idx: array-like """ idx = bin1d_vec(lons, self.xs) idy = bin1d_vec(lats, self.ys) # handles the case where values are outside of the region bad_idx = numpy.where((idx == -1) | (idy == -1)) mask = self.bbox_mask[idy, idx].astype(bool) # manually set values outside region mask[bad_idx] = True return mask def get_cartesian(self, data): """Returns 2d ndrray representation of the data set, corresponding to the bounding box. Args: data: """ assert len(data) == len(self.polygons) results = numpy.zeros(self.bbox_mask.shape[:2]) ny = len(self.ys) nx = len(self.xs) for i in range(ny): for j in range(nx): if self.bbox_mask[i, j] == 0: idx = int(self.idx_map[i, j]) results[i, j] = data[idx] else: results[i, j] = numpy.nan return results def get_bbox(self): """ Returns rectangular bounding box around region. """ return (self.xs.min(), self.xs.max()+self.dh, self.ys.min(), self.ys.max()+self.dh) def midpoints(self): """ Returns midpoints of rectangular polygons in region """ return numpy.array([poly.centroid() for poly in self.polygons]) def origins(self): """ Returns origins of rectangular polygons in region """ return numpy.array([poly.origin for poly in self.polygons]) def to_dict(self): adict = { 'name': str(self.name), 'dh': float(self.dh), 'polygons': [{'lat': float(poly.origin[1]), 'lon': float(poly.origin[0])} for poly in self.polygons], 'class_id': self.__class__.__name__ } return adict @classmethod def from_dict(cls, adict): """ Creates a region object from a dictionary """ origins = adict.get('polygons', None) dh = adict.get('dh', None) magnitudes = adict.get('magnitudes', None) name = adict.get('name', 'CartesianGrid2D') if origins is None: raise AttributeError("cannot create region object without origins") if dh is None: raise AttributeError("cannot create region without dh") if origins is not None: try: origins = numpy.array([[adict['lon'], adict['lat']] for adict in origins]) except: raise TypeError('origins must be numpy array like.') if magnitudes is not None: try: magnitudes = numpy.array(magnitudes) except: raise TypeError('magnitudes must be numpy array like.') out = cls.from_origins(origins, dh=dh, magnitudes=magnitudes, name=name) return out @classmethod def from_origins(cls, origins, dh=None, magnitudes=None, name=None): """Creates instance of class from 2d numpy.array of lon/lat origins. Note: Grid spacing should be constant in the entire region. This condition is not explicitly checked for for performance reasons. Args: origins (numpy.ndarray like): [:,0] = lons and [:,1] = lats magnitudes (numpy.array like): optional, if provided will bind magnitude information to the class. Returns: cls """ # ensure we can access the lons and lats try: lons = origins[:,0] lats = origins[:,1] except (TypeError): raise TypeError("origins must be of type numpy.array or be numpy array like.") # dh must be regular, no explicit checking. if dh is None: dh2 = numpy.abs(lons[1]-lons[0]) dh1 = numpy.abs(lats[1]-lats[0]) dh = numpy.max([dh1, dh2]) region = CartesianGrid2D([Polygon(bbox) for bbox in compute_vertices(origins, dh)], dh, name=name) if magnitudes is not None: region.magnitudes = magnitudes return region def _build_bitmask_vec(self): """ same as build mask but using vectorized calls to bin1d """ # build bounding box of set of polygons based on origins nd_origins = numpy.array([poly.origin for poly in self.polygons]) bbox = [(numpy.min(nd_origins[:, 0]), numpy.min(nd_origins[:, 1])), (numpy.max(nd_origins[:, 0]), numpy.max(nd_origins[:, 1]))] # get midpoints for hashing midpoints = numpy.array([poly.centroid() for poly in self.polygons]) # set up grid over bounding box xs = cleaner_range(bbox[0][0], bbox[1][0], self.dh) ys = cleaner_range(bbox[0][1], bbox[1][1], self.dh) # set up mask array, 1 is index 0 is mask a = numpy.ones([len(ys), len(xs), 2]) # set all indices to nan a[:, :, 1] = numpy.nan # bin1d returns the index of polygon within the cartesian grid idx = bin1d_vec(midpoints[:, 0], xs) idy = bin1d_vec(midpoints[:, 1], ys) for i in range(len(self.polygons)): a[idy[i], idx[i], 1] = int(i) # build mask in dim=0; here masked values are 1. see note below. if idx[i] >= 0 and idy[i] >= 0: if self.poly_mask is not None: # note: csep1 gridded forecast file format convention states that a "1" indicates a valid cell, which is the opposite # of the masking criterion if self.poly_mask[i] == 1: a[idy[i], idx[i], 0] = 0 else: a[idy[i], idx[i], 0] = 0 return a, xs, ys def tight_bbox(self, precision=4): # creates tight bounding box around the region poly = np.array([i.points for i in self.polygons]) sorted_idx = np.sort(np.unique(poly, return_index=True, axis=0)[1], kind='stable') unique_poly = poly[sorted_idx] # merges all the cell polygons into one polygons = [geometry.Polygon(np.round(i, precision)) for i in unique_poly] joined_poly = unary_union(polygons) bounds = np.array([i for i in joined_poly.boundary.xy]).T return bounds def get_cell_area(self): """ Compute the area of each polygon in sq. kilometers. Returns: out (numpy.array): numpy array containing cell area in km^2 """ area = numpy.zeros(self.num_nodes) for idx, origin in enumerate(self.origins()): top_right = origin + self.dh area[idx] = geographical_area_from_bounds(origin[0], origin[1], top_right[0], top_right[1]) return area
def geographical_area_from_bounds(lon1, lat1, lon2, lat2): """ Computes area of spatial cell identified by origin coordinate and top right cooridnate. The functions computes area only for square/rectangle bounding box by based on spherical earth assumption. Args: lon1,lat1 : Origin coordinates lon2,lat2: Top right coordinates Returns: Area of cell in Km2 """ if lon1 == lon2 or lat1 == lat2: return 0 else: earth_radius_km = 6371. R2 = earth_radius_km ** 2 rad_per_deg = numpy.pi / 180.0e0 strip_area_steradian = 2 * numpy.pi * (1.0e0 - numpy.cos((90.0e0 - lat1) * rad_per_deg)) \ - 2 * numpy.pi * (1.0e0 - numpy.cos((90.0e0 - lat2) * rad_per_deg)) area_km2 = strip_area_steradian * R2 / (360.0 / (lon2 - lon1)) return area_km2 def quadtree_grid_bounds(quadk): """ Computes the bottom-left and top-right coordinates corresponding to every quadkey Args: qk : Array of Strings Quadkeys. Returns: grid_coords : Array of floats [lon1,lat1,lon2,lat2] """ origin_lat = [] origin_lon = [] top_right_lon = [] top_right_lat = [] for i in range(len(quadk)): origin_lon.append(mercantile.bounds(mercantile.quadkey_to_tile(quadk[i])).west) origin_lat.append(mercantile.bounds(mercantile.quadkey_to_tile(quadk[i])).south) top_right_lon.append(mercantile.bounds(mercantile.quadkey_to_tile(quadk[i])).east) top_right_lat.append(mercantile.bounds(mercantile.quadkey_to_tile(quadk[i])).north) grid_origin = numpy.column_stack((numpy.array(origin_lon), numpy.array(origin_lat))) grid_top_right = numpy.column_stack((numpy.array(top_right_lon), numpy.array(top_right_lat))) grid_bounds = numpy.column_stack((grid_origin, grid_top_right)) return grid_bounds def compute_vertex_bounds(bound_point, tol=numpy.finfo(float).eps): """ Wrapper function to compute vertices using bounding points for multiple points. Default tolerance is set to machine precision of floating point number. Args: bounding points: nx4 ndarray [lon_origin, lat_origin, lon_top_right, lat_origin] Notes: (x,y) should be accessible like: #>>> origin coords = origin_points[:,0:1] #>>> Top right coords = origin_points[:,2:3] """ bbox = ((bound_point[0], bound_point[1]), (bound_point[0], bound_point[3] - tol), (bound_point[2] - tol, bound_point[3] - tol), (bound_point[2] - tol, bound_point[1])) return bbox def compute_vertices_bounds(bounds, tol=numpy.finfo(float).eps): """ Wrapper function to compute vertices using bounding points for multiple points. Default tolerance is set to machine precision of floating point number. Args: bounding points: nx4 ndarray [lon_origin, lat_origin, lon_top_right, lat_origin] Notes: (x,y) should be accessible like: #>>> origin coords = origin_points[:,0:1] #>>> Top right coords = origin_points[:,2:3] """ return list(map(lambda x: compute_vertex_bounds(x, tol=tol), bounds)) def _create_tile(quadk, threshold, zoom, lon, lat, qk, num): """ **Alert: This Function uses GLOBAL variable (qk) and (num). Provides multi-resolution quadtree spatial grid based on seismic density. It takes in a starting quadtree Tile (Quadkey), then keeps on increasing the zoom-level of every Tile (or dividing cell) recursively, unless every cell meets the cell dividion criteria. The primary criterion of dividing a parent cell into 4 child cells is a threshold on seismic denisity. The cells are divided unless evevry cell cas number of earthquakes less than "threshold". The cell division of any also stops if it reaches maximum zoom-level (zoom) Args: quadk : String 0, 1, 2, 3 or any desired starting level of Quad key. threshold : int Max number of earthquakes/cell allowed zoom: int Maximum zoom level allowed for a quadkey lon : float longitudes of earthquakes in catalog lat : float latitude of earthquakes in catalog Returns: """ boundary = mercantile.bounds(mercantile.quadkey_to_tile(quadk)) eqs = numpy.logical_and(numpy.logical_and(lon >= boundary.west, lat >= boundary.south), numpy.logical_and(lon < boundary.east, lat < boundary.north)) num_eqs = numpy.size(lat[eqs]) # global qk # global num # Setting the Min Threshold of Area 1 sq. km. Instead of Depth if num_eqs > threshold and len(quadk) < zoom: # #qk_area_km(quadk)>4: # print('inside If, Current Quad key ', quadk) # print('Length of Quadkey ', len(quadk)) # # print('Num of Eqs ', num_eqs) _create_tile(quadk + '0', threshold, zoom, lon, lat, qk, num) _create_tile(quadk + '1', threshold, zoom, lon, lat, qk, num) _create_tile(quadk + '2', threshold, zoom, lon, lat, qk, num) _create_tile(quadk + '3', threshold, zoom, lon, lat, qk, num) else: # print('inside ELSE, Current Quad key ', quadk) # print('Num of Eqs ', num_eqs) # qk = numpy.append(qk, quadk) qk.append(quadk) # num = numpy.append(num, num_eqs) num.append(num_eqs) def _create_tile_fix_len(quadk, zoom, qk): """ ***Alert: This Function uses GLOBAL variable (qk). Provides single-resolution quadtree grid. It takes in a starting quadkey (or Quadrant of Globe), then keeps on keeps on dividing it into 4 children unless the maximum zoom-level is achieved Parameters ---------- quadk : String 0, 1, 2, 3 or any desired starting level of Quad key. zoom : TYPE Length of Quad Key OR Depth of grid. Returns ------- None. """ if len(quadk) < zoom: # print('inside If, Current Quad key ', quadk) # print('Len of QK: ', len(quadk)) _create_tile_fix_len(quadk + '0', zoom, qk) _create_tile_fix_len(quadk + '1', zoom, qk) _create_tile_fix_len(quadk + '2', zoom, qk) _create_tile_fix_len(quadk + '3', zoom, qk) else: # print('inside ELSE, Current Quad key ', quadk) # print('Num of Eqs ', num_eqs) # qk = numpy.append(qk, quadk) qk.append(quadk) class QuadtreeGrid2D: """ Respresents a 2D quadtree gridded region. The class provides functionality to generate multi-resolution or single-resolution quadtree grid. It also enables users to load already available quadtree grird. It also provides functions to query onto an index 2D grid ad maintains mapping between space coordinates and defined polygons and the index into the polygon array. Note: It is replica of CartesianGrid2D class but with quadtree approach, with implementation of all the relevant functions required to CSEP1 tests """ def __init__(self, polygons, quadkeys, bounds, name='QuadtreeGrid2d', mask=None): """ Args: polygons: Represents the object of class "polygons" defined through a collection of vertices. This polygon is 2d and vertices are obtained as corner points of quadtree tile. quadkeys: Unique identifier of each quadtree tile. Quadkeys of every tile defines a grid cell. This is the first thing computed while acquiring quadtree grid. Rest can be computed from this. bounds: number of cells x [lon1, lat1, lon2, lat2], corresponding to origin coordinates and top right coordinates fo each grid cell name: Name of grid mask: Masked cells. NotImplemented yet. Always keep it none """ self.polygons = polygons self.quadkeys = quadkeys self.bounds = bounds self.cell_area = [] self.poly_mask = mask self.name = name # a, xs, ys = self._get_idx_map_xs_ys() # self.xs = xs # self.ys = ys # self.idx_map = a @property def num_nodes(self): """ Number of polygons in region """ return len(self.polygons) def get_cell_area(self): """ Calls function geographical_area_from_bounds and computes area of each grid cell. It also modified class variable "self.cell_area" It iterates over all the cells of grid and passes bounding coordinates of every cell to function geographical_area_from_bounds """ cell_area = numpy.array([geographical_area_from_bounds(bb[0],bb[1],bb[2],bb[3]) for bb in self.bounds]) self.cell_area = cell_area return self.cell_area def get_index_of(self, lons, lats): """ Returns the index of lons, lats in self.polygons Args: lons: ndarray-like lats: ndarray-like Returns: idx: ndarray-like """ # If its array or many coords if isinstance(lons, (list, numpy.ndarray)): idx = [] for i in range(len(lons)): idx = numpy.append(idx, self._find_location(lons[i], lats[i])) idx = idx.astype(int) return idx # It its just one Lon/Lon if isinstance(lons, (int, float)): idx = self._find_location(lons, lats) return idx return None def _find_location(self, lon, lat): """ Takes in single Lon and Lat and finds its Polygon Index. Returns: index number of polyons """ loc = numpy.logical_and(numpy.logical_and(lon >= self.bounds[:, 0], lat >= self.bounds[:, 1]), numpy.logical_and(lon < self.bounds[:, 2], lat < self.bounds[:, 3])) if len(numpy.where(loc == True)[0]) > 0: return numpy.where(loc == True)[0][0] else: return numpy.where(loc == True)[0] def get_location_of(self, indices): """ Returns the polygon associated with the index idx. Args: idx: index of polygon in region Returns: Polygon """ indices = list(indices) polys = [self.polygons[idx] for idx in indices] return polys def _get_spatial_counts(self, catalog, mag_bins=None): """ Gets the number of earthquakes in each cell for available catalog. Uses QuadtreeGrid2D.get_index_of function to map every earthquake location to its corresponding cell Args: catalog: CSEP Catalog mag_bins: Magnitude discritization used in earthquake forecast mdoel Note: mag_bins are only required to filter catalog for minimum magnitude Return: spatial counts: Number of earthquakes in each cell """ if mag_bins is None or mag_bins == []: mag_bins = catalog.magnitudes if min(catalog.get_magnitudes()) < min(mag_bins): print("-----Warning-----") print("Catalog contains magnitudes below the min magnitude range") print("Filtering catalog with Magnitude: ", min(mag_bins)) catalog.filter('magnitude >= ' + str(min(mag_bins))) if min(catalog.get_latitudes()) < self.get_bbox()[2] or max(catalog.get_latitudes()) > self.get_bbox()[3]: print("----Warning---") print("Catalog exceeds grid bounds, so catalog filtering") catalog.filter('latitude < ' + str(self.get_bbox()[3])) catalog.filter('latitude > ' + str(self.get_bbox()[2])) lon = catalog.get_longitudes() lat = catalog.get_latitudes() out = numpy.zeros(len(self.quadkeys)) idx = self.get_index_of(lon, lat) numpy.add.at(out, idx, 1) return out def _get_spatial_magnitude_counts(self, catalog, mag_bins=None): """ Gets the number of earthquakes in for each spatio-magnitude bin for available catalog Uses QuadtreeGrid2D.get_index_of function to map every earthquake location to its corresponding cell Uses bin1d_vec function to map earthquake magnitude to its respecrtive bin. Args: catalog: CSEPCatalog mag_bins: Magnitude discritization used in earthquake forecast model Return: Spatial-magnitude counts """ if mag_bins is None or mag_bins == []: mag_bins = catalog.magnitudes if min(catalog.get_magnitudes()) < min(mag_bins): print("-----Warning-----") print("Catalog contains magnitudes below the min magnitude range") print("Filtering catalog with Magnitude: ", min(mag_bins)) catalog.filter('magnitude >= ' + str(min(mag_bins))) if min(catalog.get_latitudes()) < self.get_bbox()[2] or max(catalog.get_latitudes()) > self.get_bbox()[3]: print("----Warning---") print("Catalog exceeds grid bounds filtering events outside of the region boundary") catalog.filter('latitude < ' + str(self.get_bbox()[3])) catalog.filter('latitude > ' + str(self.get_bbox()[2])) lon = catalog.get_longitudes() lat = catalog.get_latitudes() mag = catalog.get_magnitudes() out = numpy.zeros([len(self.quadkeys), len(mag_bins)]) idx_loc = self.get_index_of(lon, lat) idx_mag = bin1d_vec(mag, mag_bins, tol=0.00001, right_continuous=True) numpy.add.at(out, (idx_loc, idx_mag), 1) return out def get_bbox(self): """ Returns rectangular bounding box around region. """ # return (self.xs.min(), self.xs.max(), self.ys.min(), self.ys.max()) return (min(self.bounds[:, 0]), max(self.bounds[:, 2]), min(self.bounds[:, 1]), max(self.bounds[:, 3])) def midpoints(self): """ Returns midpoints of rectangular polygons in region """ return numpy.array([poly.centroid() for poly in self.polygons]) def origins(self): """ Returns origins of rectangular polygons in region """ return numpy.array([poly.origin for poly in self.polygons]) def to_dict(self): adict = { 'name': str(self.name), 'polygons': [{'lat': float(poly.origin[1]), 'lon': float(poly.origin[0])} for poly in self.polygons] } return adict def save_quadtree(self, filename): """ Saves the quadtree grid (quadkeys) in a text file Args: filename (str): filename to store file """ numpy.savetxt(filename, self.quadkeys, delimiter=',', fmt='%s') @classmethod def from_catalog(cls, catalog, threshold, zoom=11, magnitudes=None, name=None): """ Creates instance of class from 2d numpy.array of lon/lat of Catalog. Provides multi-resolution quadtree spatial grid based on seismic density. It starts from whole globe as 4 cells (Quadkeys:'0','1','2','3'), then keeps on increasing the zoom-level of every Tile recursively, unless every cell meets the division criteria. The primary criterion of dividing a parent cell into 4 child cells is a threshold on seismic density. The cells are divided unless every cell has number of earthquakes less than "threshold". The division of a cell also stops if it reaches maximum zoom-level (zoom) Args: catalog (CSEPCatalog): catalog used to create quadtree threshold (int): Max earthquakes allowed per cells zoom (int): Max zoom allowed for a cell magnitudes (array-like): left end values of magnitude discretization Returns: instance of QuadtreeGrid2D """ lon = catalog.get_longitudes() lat = catalog.get_latitudes() qk = [] num = [] _create_tile('0', threshold, zoom, lon, lat, qk, num) _create_tile('1', threshold, zoom, lon, lat, qk, num) _create_tile('2', threshold, zoom, lon, lat, qk, num) _create_tile('3', threshold, zoom, lon, lat, qk, num) qk = numpy.array(qk) bounds = quadtree_grid_bounds(qk) region = QuadtreeGrid2D( [Polygon(bbox) for bbox in compute_vertices_bounds(bounds)], qk, bounds, name=name) if magnitudes is not None: region.magnitudes = magnitudes return region @classmethod def from_single_resolution(cls, zoom, magnitudes=None, name=None): """ Creates instance of class at single-resolution using provided zoom-level. Provides single-resolution quadtree grid. It starts from whole globe as 4 cells (Quadkeys:'0','1','2','3'), then keeps on keeps on dividing every cell into 4 children unless the maximum zoom-level is achieved Args: zoom: Max zoom allowed for a cell magnitude: magnitude discretization Returns: instance of QuadtreeGrid2D """ qk = [] _create_tile_fix_len('0', zoom, qk) _create_tile_fix_len('1', zoom, qk) _create_tile_fix_len('2', zoom, qk) _create_tile_fix_len('3', zoom, qk) qk = numpy.array(qk) bounds = quadtree_grid_bounds(qk) region = QuadtreeGrid2D([Polygon(bbox) for bbox in compute_vertices_bounds(bounds)], qk, bounds, name=name) if magnitudes is not None: region.magnitudes = magnitudes return region @classmethod def from_quadkeys(cls, quadk, magnitudes=None, name=None): """ Creates instance of class from available quadtree grid. Args: quadk (list): List of quad keys strings corresponding to an already available quadtree grid magnitudes (array-like): left end-points of magnitude discretization Returns: instance of QuadtreeGrid2D """ bounds = quadtree_grid_bounds(numpy.array(quadk)) region = QuadtreeGrid2D([Polygon(bbox) for bbox in compute_vertices_bounds(bounds)], quadk, bounds, name=name) if magnitudes is not None: region.magnitudes = magnitudes return region def _get_idx_map_xs_ys(self): print('inside _get_idx_map') nd_origins = numpy.array([poly.origin for poly in self.polygons]) xs = numpy.unique(nd_origins[:, 0]) ys = numpy.unique(nd_origins[:, 1]) ny = len(ys) nx = len(xs) #Get the index map a = numpy.zeros([ny, nx]) for i in range(nx): for j in range(ny): idx = self.get_index_of(xs[i], ys[j]) a[j, i] = idx return a, xs, ys def get_cartesian(self, data): """ Returns 2d ndrray representation of the data set, corresponding to the bounding box. Args: data (numpy.array): array of values corresponding to cells in the quadtree region Returns: results (numpy.array): 2d numpy array with rates on cartesian grid """ a, xs, ys = self._get_idx_map_xs_ys() self.xs = xs self.ys = ys self.idx_map = a assert len(data) == len(self.polygons) ny = len(self.ys) nx = len(self.xs) results = numpy.zeros([ny, nx]) for i in range(nx): for j in range(ny): idx = int(self.idx_map[j,i]) results[j, i] = data[idx] return results def tight_bbox(self): # creates tight bounding box around the region, probably a faster way to do this. ny, nx = self.idx_map.shape asc = [] desc = [] for j in range(ny): row = self.idx_map[j, :] argmin = first_nonnan(row) argmax = last_nonnan(row) # points are stored clockwise poly_min = self.polygons[int(row[argmin])].points asc.insert(0, poly_min[0]) asc.insert(0, poly_min[1]) poly_max = self.polygons[int(row[argmax])].points lat_0 = poly_max[2][1] lat_1 = poly_max[3][1] # last two points are 'right hand side of polygon' if lat_0 < lat_1: desc.append(poly_max[2]) desc.append(poly_max[3]) else: desc.append(poly_max[3]) desc.append(poly_max[2]) # close the loop poly = np.array(asc + desc) sorted_idx = np.sort(np.unique(poly, return_index=True, axis=0)[1], kind='stable') unique_poly = poly[sorted_idx] unique_poly = np.append(unique_poly, [unique_poly[0, :]], axis=0) return unique_poly def california_quadtree_region(magnitudes=None, name="california-quadtree"): """ Returns object of QuadtreeGrid2D representing quadtree grid for California RELM testing region. The grid is already generated at zoom-level = 12 and it is loaded through classmethod: QuadtreeGrid2D.from_quadkeys The grid cells at zoom level 12 are selected using the external boundary of RELM california region. This grid can be used to create gridded datasets for earthquake forecasts. Args: magnitudes: Magnitude discretization name: string Returns: :class:`csep.core.spatial.QuadtreeGrid2D """ # use default file path from python package root_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) filepath = os.path.join(root_dir, 'artifacts', 'Regions', 'california_qk_zoom=12.txt') qk = numpy.genfromtxt(filepath, delimiter=',', dtype='str') california_region = QuadtreeGrid2D.from_quadkeys(qk, magnitudes=magnitudes, name=name) return california_region