Source code for obspy.io.shapefile.core

# -*- coding: utf-8 -*-
import datetime
import re
import warnings

from obspy import Catalog, UTCDateTime
from obspy.core.event import Origin, Magnitude
from obspy.core.inventory import Inventory
from obspy.core.util.misc import to_int_or_zero

try:
    import shapefile
except ImportError as e:
    HAS_PYSHP = False
    PYSHP_VERSION = None
    PYSHP_VERSION_AT_LEAST_1_2_11 = False
    IMPORTERROR_MSG = str(e) + (
        ". ObsPy's write support for shapefiles requires the 'pyshp' module "
        "to be installed in addition to the general ObsPy dependencies.")
else:
    HAS_PYSHP = True
    try:
        PYSHP_VERSION = list(map(to_int_or_zero,
                                 shapefile.__version__.split('.')))
    except AttributeError:
        PYSHP_VERSION = None
        PYSHP_VERSION_AT_LEAST_1_2_11 = False
    else:
        PYSHP_VERSION_AT_LEAST_1_2_11 = PYSHP_VERSION >= [1, 2, 11]
PYSHP_VERSION_WARNING = (
    'pyshp versions < 1.2.11 are buggy, e.g. in writing numerical values to '
    'the dbf table, so e.g. timestamp float values might lack proper '
    'precision. You should update to a newer pyshp version.')


[docs] def _write_shapefile(obj, filename, extra_fields=None, **kwargs): """ Write :class:`~obspy.core.inventory.inventory.Inventory` or :class:`~obspy.core.event.Catalog` object to a ESRI shapefile. :type obj: :class:`~obspy.core.event.Catalog` or :class:`~obspy.core.inventory.inventory.Inventory` :param obj: ObsPy object for shapefile output :type filename: str :param filename: Filename to write to. According to ESRI shapefile definition, multiple files with the following suffixes will be written: ".shp", ".shx", ".dbj", ".prj". If filename does not end with ".shp", it will be appended. Other files will be created with respective suffixes accordingly. :type extra_fields: list :param extra_fields: List of extra fields to write to the shapefile table. Each item in the list has to be specified as a tuple of: field name (i.e. name of database column, ``str``), field type (single character as used by ``pyshp``: ``'C'`` for string fields, ``'N'`` for integer/float fields - use precision ``None`` for integer fields, ``'L'`` for boolean fields), field width (``int``), field precision (``int``) and field values (``list`` of individual values, must have same length as given catalog object or as the sum of all station objects across all networks of a given inventory). """ if not HAS_PYSHP: raise ImportError(IMPORTERROR_MSG) if not PYSHP_VERSION_AT_LEAST_1_2_11: warnings.warn(PYSHP_VERSION_WARNING) if not filename.endswith(".shp"): filename += ".shp" if PYSHP_VERSION >= [2., 0, 0]: writer = shapefile.Writer(target=filename, shapeType=shapefile.POINT) else: writer = shapefile.Writer(shapeType=shapefile.POINT) writer.autoBalance = 1 # create the layer if isinstance(obj, Catalog): _add_catalog_layer(writer, obj, extra_fields=extra_fields) elif isinstance(obj, Inventory): _add_inventory_layer(writer, obj, extra_fields=extra_fields) else: msg = ("Object for shapefile output must be " "a Catalog or Inventory.") raise TypeError(msg) if PYSHP_VERSION >= [2.0, 0, 0]: writer.close() else: writer.save(filename) _save_projection_file(filename.rsplit('.', 1)[0] + '.prj')
[docs] def _add_catalog_layer(writer, catalog, extra_fields=None): """ :type writer: :class:`shapefile.Writer`. :param writer: pyshp Writer object :type catalog: :class:`~obspy.core.event.Catalog` :param catalog: Event data to add as a new layer. :type extra_fields: list :param extra_fields: List of extra fields to write to the shapefile table. For details see :func:`_write_shapefile()`. """ # [name, type, width, precision] # field name is 10 chars max # ESRI shapefile attributes are stored in dbf files, which can not # store datetimes, only dates, see: # http://www.gdal.org/drv_shapefile.html # use POSIX timestamp for exact origin time, set time of first pick # for events with no origin field_definitions = [ ["EventID", 'C', 100, None], ["OriginID", 'C', 100, None], ["MagID", 'C', 100, None], ["Date", 'D', None, None], ["OriginTime", 'N', 20, 6], ["FirstPick", 'N', 20, 6], ["Longitude", 'N', 16, 10], ["Latitude", 'N', 16, 10], ["Depth", 'N', 8, 3], ["MinHorUncM", 'N', 12, 3], ["MaxHorUncM", 'N', 12, 3], ["MaxHorAzi", 'N', 7, 3], ["OriUncDesc", 'C', 40, None], ["Magnitude", 'N', 8, 3], ] _create_layer(writer, field_definitions, extra_fields) if extra_fields: for name, type_, width, precision, values in extra_fields: if len(values) != len(catalog): msg = ("list of values for each item in 'extra_fields' must " "have same length as Catalog object") raise ValueError(msg) for i, event in enumerate(catalog): # try to use preferred origin/magnitude, fall back to first or use # empty one with `None` values in it origin = (event.preferred_origin() or event.origins and event.origins[0] or Origin(force_resource_id=False)) magnitude = (event.preferred_magnitude() or event.magnitudes and event.magnitudes[0] or Magnitude(force_resource_id=False)) t_origin = origin.time pick_times = [pick.time for pick in event.picks if pick.time is not None] t_pick = pick_times and min(pick_times) or None date = t_origin or t_pick feature = {} # setting fields with `None` results in values of `0.000` # need to really omit setting values if they are `None` if event.resource_id is not None: feature["EventID"] = str(event.resource_id) if origin.resource_id is not None: feature["OriginID"] = str(origin.resource_id) if t_origin is not None: # Use timestamp for exact timing feature["OriginTime"] = t_origin.timestamp if t_pick is not None: # Use timestamp for exact timing feature["FirstPick"] = t_pick.timestamp if date is not None: # ESRI shapefile attributes are stored in dbf files, which can # not store datetimes, only dates. We still need to use the # GDAL API with precision up to seconds (aiming at other output # drivers of GDAL; `100` stands for GMT) feature["Date"] = date.datetime if origin.latitude is not None: feature["Latitude"] = origin.latitude if origin.longitude is not None: feature["Longitude"] = origin.longitude if origin.depth is not None: feature["Depth"] = origin.depth / 1e3 if magnitude.mag is not None: feature["Magnitude"] = magnitude.mag if magnitude.resource_id is not None: feature["MagID"] = str(magnitude.resource_id) if origin.origin_uncertainty is not None: ou = origin.origin_uncertainty ou_description = ou.preferred_description if ou_description == 'uncertainty ellipse': feature["MinHorUncM"] = ou.min_horizontal_uncertainty feature["MaxHorUncM"] = ou.max_horizontal_uncertainty feature["MaxHorAzi"] = \ ou.azimuth_max_horizontal_uncertainty feature["OriUncDesc"] = ou_description elif ou_description == 'horizontal uncertainty': feature["MinHorUncM"] = ou.horizontal_uncertainty feature["MaxHorUncM"] = ou.horizontal_uncertainty feature["MaxHorAzi"] = 0.0 feature["OriUncDesc"] = ou_description else: msg = ('Encountered an event with origin uncertainty ' 'description of type "{}". This is not yet ' 'implemented for output as shapefile. No origin ' 'uncertainty will be added to shapefile for such ' 'events.').format(ou_description) warnings.warn(msg) if origin.latitude is not None and origin.longitude is not None: writer.point(origin.longitude, origin.latitude) if extra_fields: for name, _, _, _, values in extra_fields: feature[name] = values[i] _add_record(writer, feature)
[docs] def _add_inventory_layer(writer, inventory, extra_fields=None): """ :type writer: :class:`shapefile.Writer`. :param writer: pyshp Writer object :type inventory: :class:`~obspy.core.inventory.inventory.Inventory` :param inventory: Inventory data to add as a new layer. :type extra_fields: list :param extra_fields: List of extra fields to write to the shapefile table. For details see :func:`_write_shapefile()`. """ # [name, type, width, precision] # field name is 10 chars max # ESRI shapefile attributes are stored in dbf files, which can not # store datetimes, only dates, see: # http://www.gdal.org/drv_shapefile.html # use POSIX timestamp for exact origin time, set time of first pick # for events with no origin field_definitions = [ ["Network", 'C', 20, None], ["Station", 'C', 20, None], ["Longitude", 'N', 16, 10], ["Latitude", 'N', 16, 10], ["Elevation", 'N', 9, 3], ["StartDate", 'D', None, None], ["EndDate", 'D', None, None], ["Channels", 'C', 254, None], ] _create_layer(writer, field_definitions, extra_fields) station_count = sum(len(net) for net in inventory) if extra_fields: for name, type_, width, precision, values in extra_fields: if len(values) != station_count: msg = ("list of values for each item in 'extra_fields' must " "have same length as the count of all Stations " "combined across all Networks.") raise ValueError(msg) i = 0 for net in inventory: for sta in net: channel_list = ",".join(["%s.%s" % (cha.location_code, cha.code) for cha in sta]) feature = {} # setting fields with `None` results in values of `0.000` # need to really omit setting values if they are `None` if net.code is not None: feature["Network"] = net.code if sta.code is not None: feature["Station"] = sta.code if sta.latitude is not None: feature["Latitude"] = sta.latitude if sta.longitude is not None: feature["Longitude"] = sta.longitude if sta.elevation is not None: feature["Elevation"] = sta.elevation if sta.start_date is not None: # ESRI shapefile attributes are stored in dbf files, which # can not store datetimes, only dates. We still need to use # the GDAL API with precision up to seconds (aiming at # other output drivers of GDAL; `100` stands for GMT) feature["StartDate"] = sta.start_date.datetime if sta.end_date is not None: # ESRI shapefile attributes are stored in dbf files, which # can not store datetimes, only dates. We still need to use # the GDAL API with precision up to seconds (aiming at # other output drivers of GDAL; `100` stands for GMT) feature["EndDate"] = sta.end_date.datetime if channel_list: feature["Channels"] = channel_list if extra_fields: for name, _, _, _, values in extra_fields: feature[name] = values[i] if sta.latitude is not None and sta.longitude is not None: writer.point(sta.longitude, sta.latitude) _add_record(writer, feature) i += 1
wgs84_wkt = \ """ GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]] """ wgs84_wkt = re.sub(r'\s+', '', wgs84_wkt)
[docs] def _save_projection_file(filename): with open(filename, 'wt') as fh: fh.write(wgs84_wkt)
[docs] def _add_field(writer, name, type_, width, precision): # default field width is not set correctly for dates and booleans in # shapefile <=1.2.10, see # GeospatialPython/pyshp@ba61854aa7161fd7d4cff12b0fd08b6ec7581bb7 and # GeospatialPython/pyshp#71 so work around this if type_ == 'D': width = 8 precision = 0 elif type_ == 'L': width = 1 precision = 0 kwargs = dict(fieldType=type_, size=width, decimal=precision) # remove None's because shapefile.Writer.field() doesn't use None as # placeholder but the default values directly for key in list(kwargs.keys()): if kwargs[key] is None: kwargs.pop(key) writer.field(name, **kwargs)
[docs] def _create_layer(writer, field_definitions, extra_fields=None): # Add the fields we're interested in for name, type_, width, precision in field_definitions: _add_field(writer, name, type_, width, precision) field_names = [name for name, _, _, _ in field_definitions] # add custom fields if extra_fields is not None: for name, type_, width, precision, _ in extra_fields: if name in field_names: msg = "Conflict with existing field named '{}'.".format(name) raise ValueError(msg) _add_field(writer, name, type_, width, precision)
[docs] def _add_record(writer, feature): values = [] for key, type_, width, precision in writer.fields: value = feature.get(key) # various hacks for old pyshp < 1.2.11 if not PYSHP_VERSION_AT_LEAST_1_2_11: if type_ == 'C': # mimic pyshp 1.2.12 behavior of putting 'None' in string # fields for value of `None` if value is None: value = 'None' # older pyshp is not correctly writing dates as used nowadays # '%Y%m%d' (8 chars), work around this elif type_ == 'D': if isinstance(value, (UTCDateTime, datetime.date)): value = value.strftime('%Y%m%d') # work around issues with older pyshp, backport 1.2.12 behavior elif type_ == 'L': # logical: 1 byte - initialized to 0x20 (space) # otherwise T or F if value in [True, 1]: value = "T" elif value in [False, 0]: value = "F" else: value = ' ' # work around issues with older pyshp, backport 1.2.12 behavior elif type_ in ('N', 'F'): # numeric or float: number stored as a string, right justified, # and padded with blanks to the width of the field. if value in (None, ''): value = ' ' * width # QGIS NULL elif not precision: # caps the size if exceeds the field size value = format(value, "d")[:width].rjust(width) else: # caps the size if exceeds the field size value = format(value, ".%sf" % precision)[:width].rjust( width) # work around older pyshp not converting `None`s properly (e.g. for # float fields) elif value is None: value = '' values.append(value) writer.record(*values)
if __name__ == '__main__': import doctest doctest.testmod(exclude_empty=True)