# ------------------- The OpenQuake Model Building Toolkit --------------------
# Copyright (C) 2022 GEM Foundation
# _______ _______ __ __ _______ _______ ___ _
# | || | | |_| || _ || || | | |
# | _ || _ | ____ | || |_| ||_ _|| |_| |
# | | | || | | ||____|| || | | | | _|
# | |_| || |_| | | || _ | | | | |_
# | || | | ||_|| || |_| | | | | _ |
# |_______||____||_| |_| |_||_______| |___| |___| |_|
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Affero General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more
# details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# -----------------------------------------------------------------------------
# vim: tabstop=4 shiftwidth=4 softtabstop=4
# coding: utf-8
"""
Module :mod:`openquake.wkf.catalogue` contains tools for working with catalogue
files in the hmtk format
"""
import os
import numpy as np
import pandas as pd
import geopandas as gpd
from typing import Type
from shapely.geometry import Point
from openquake.hmtk.seismicity.catalogue import Catalogue
from openquake.wkf.utils import create_folder, get_list
from openquake.hmtk.parsers.catalogue.gcmt_ndk_parser import ParseNDKtoGCMT
[docs]
def to_df(cat: Type[Catalogue]) -> pd.DataFrame:
"""
Converts an :class:`openquake.hmtk.seismicity.catalogue.Catalogue` instance
into a dataframe
:param cat:
The catalogue instance
:returns:
The dataframe with the catalogue
"""
df = pd.DataFrame()
for key in cat.data:
if key not in ['comment', 'flag'] and len(cat.data[key]):
df.loc[:, key] = cat.data[key]
return df
[docs]
def from_df(df, end_year=None) -> Type[Catalogue]:
"""
Converts a dataframe into a
:class:`openquake.hmtk.seismicity.catalogue.Catalogue` instance
:param df:
The dataframe with the catalogue
:returns:
The catalogue instance
"""
cat = Catalogue()
for column in df:
if (column in Catalogue.FLOAT_ATTRIBUTE_LIST or
column in Catalogue.INT_ATTRIBUTE_LIST):
cat.data[column] = df[column].to_numpy()
else:
cat.data[column] = df[column]
cat.end_year = np.max(df.year) if end_year is None else end_year
return cat
[docs]
def create_subcatalogues(fname_polygons: str, fname_cat: str, folder_out: str,
source_ids: str = []):
"""
Given a catalogue and a gis-file with polygons (e.g. shapefile or
.geojson), this code creates for each polygon a subcatalogue with the
earthquakes with epicenters in the polygon.
:param fname_polygons:
The name of the gis file containing the polygons.
:param fname_cat:
The name of the file with the catalogue (hmtk formatted)
:param folder_out:
The name of the output folder where to create the output .csv files
containing the subcatalogues
:param source_ids:
[optional] The list of source ids to be considered. If omitted all the
polygons will be considered.
"""
# Create output folder
create_folder(folder_out)
if len(source_ids) > 0:
source_ids = get_list(source_ids)
# Create geodataframe with the catalogue
df = pd.read_csv(fname_cat)
gdf = gpd.GeoDataFrame(df, crs='epsg:4326', geometry=[Point(xy) for xy
in zip(df.longitude, df.latitude)])
# Read polygons
polygons_gdf = gpd.read_file(fname_polygons)
# check there are no duplicate ids
assert(polygons_gdf['id'].is_unique)
# explode for any multipolygons
polygons_gdf = polygons_gdf.explode()
unique_srcs = np.unique(polygons_gdf.id)
# Iterate over sources to catch any multipolygons
# These can be sneaky! Might not be clear in QGIS
out_fnames = []
for src in unique_srcs:
polys = polygons_gdf[polygons_gdf.id == src]
# Check if we care about this source
if len(source_ids) > 0 and src not in source_ids:
continue
all_within = gpd.GeoDataFrame(columns=gdf.columns, geometry='geometry')
for idx, poly in polys.iterrows():
df = pd.DataFrame({'Name': [poly.id], 'Polygon': [poly.geometry]})
gdf_poly = gpd.GeoDataFrame(df, geometry='Polygon', crs='epsg:4326')
within = gpd.sjoin(gdf, gdf_poly, predicate='within')
all_within = pd.concat([all_within, within])
# Create output file
if isinstance(poly.id, int):
fname = f'subcatalogue_zone_{poly.id:d}.csv'
else:
fname = f'subcatalogue_zone_{poly.id}.csv'
out_fname = os.path.join(folder_out, fname)
out_fnames.append(out_fname)
all_within.to_csv(out_fname, index=False)
return out_fnames
[docs]
def get_dataframe(fname: str) -> pd.DataFrame:
"""
Creates a dataframe with the information included in a .ndk formatted
file. For a description of the .ndk format see: https://www.globalcmt.org/
:param fname:
Name of the .ndk file
:returns:
A dataframe with the information in the .ndk file
"""
parser = ParseNDKtoGCMT(fname)
cat_gcmt = parser.read_file()
df = pd.DataFrame({k: cat_gcmt.data[k] for k in cat_gcmt.data.keys()})
return df
[docs]
def create_gcmt_files(fname_polygons: str, gcmt_filename: str, folder_out: str,
depth_max: float = 600.0, depth_min: float = 0.0):
# Create output folder
create_folder(folder_out)
# Create geodataframe with the catalogue
print(os.path.abspath(gcmt_filename))
tmp = get_dataframe(gcmt_filename)
# Filter depths
df = tmp[(tmp.depth > depth_min) & (tmp.depth <= depth_max)]
if len(df) < 0:
return []
# Create geodataframe
gdf = gpd.GeoDataFrame(df, crs='epsg:4326', geometry=[Point(xy) for xy
in zip(df.longitude, df.latitude)])
# Read polygons
polygons_gdf = gpd.read_file(fname_polygons)
# Iterate over sources
fnames_list = []
for idx, poly in polygons_gdf.iterrows():
df = pd.DataFrame({'Name': [poly.id], 'Polygon': [poly.geometry]})
gdf_poly = gpd.GeoDataFrame(df, geometry='Polygon', crs='epsg:4326')
within = gpd.sjoin(gdf, gdf_poly, predicate='within')
if len(df) < 1:
continue
# Create output file
if isinstance(poly.id, int):
fname = 'subcatalogue_zone_{:d}.csv'.format(poly.id)
else:
fname = 'subcatalogue_zone_{:s}.csv'.format(poly.id)
out_fname = os.path.join(folder_out, fname)
within.to_csv(out_fname, index=False)
fnames_list.append(out_fname)
return fnames_list