2025/bin/gfs2fep

#!/usr/bin/env seiscomp-python
# -*- coding: utf-8 -*-

############################################################################
# Copyright (C) gempa GmbH                                                 #
# All rights reserved.                                                     #
# Contact: gempa GmbH (seiscomp-dev@gempa.de)                              #
#                                                                          #
# GNU Affero General Public License Usage                                  #
# This file may be used under the terms of the GNU Affero                  #
# Public License version 3.0 as published by the Free Software Foundation  #
# and appearing in the file LICENSE included in the packaging of this      #
# file. Please review the following information to ensure the GNU Affero   #
# Public License version 3.0 requirements will be met:                     #
# https://www.gnu.org/licenses/agpl-3.0.html.                              #
#                                                                          #
# Other Usage                                                              #
# Alternatively, this file may be used in accordance with the terms and    #
# conditions contained in a signed written agreement between you and       #
# gempa GmbH.                                                              #
############################################################################


import datetime
import getopt
import sys

from typing import TextIO

from seiscomp import geo


# -----------------------------------------------------------------------------
def printHelp():
    msg = """
gfs2fep - converts a SeisComP GeoFeatureSet file (GeoJSON or BNA) to FEP format

usage: {} [OPTIONS]

    -h, --help
        print this help message

    -i, --input
        input file (default: -)

    -o, --output
        output fep file (default: -)

    -a, --append
        append fep data to output file instead of overriding it

    -p, --precision (default: unrestricted)
        number of decimal places of coordintes"""
    print(msg.format(sys.argv[0]), file=sys.stderr)
    sys.exit(0)


# -----------------------------------------------------------------------------
def error(code, msg):
    print(f"error ({code}): {msg}", file=sys.stderr)
    sys.exit(code)


# -----------------------------------------------------------------------------
def run():
    if len(sys.argv) == 1:
        printHelp()

    inFile = "-"
    outFile = None
    append = False
    precision = None
    opts, _ = getopt.getopt(
        sys.argv[1:], "hi:o:ap:", ["help", "input=", "output=", "append", "precision"]
    )
    for o, a in opts:
        if o in ("-h", "--help"):
            printHelp()

        if o in ("-i", "--input"):
            inFile = a

        if o in ("-o", "--output"):
            outFile = a

        if o in ("-a", "--append"):
            append = True

        if o in ("-p", "--precision"):
            precision = max(int(a), 0)

    gfs = geo.GeoFeatureSet()
    if not gfs.readFile(inFile, None):
        error(1, f"Could not read from file '{inFile}'")

    # combine features sharing the same name
    featureDict = {}
    for f in gfs.features():
        if not f.closedPolygon():
            print(
                f"warning: feature not a closed polygon: {f.name()}",
                file=sys.stderr,
            )
        if f.name() in featureDict:
            featureDict[f.name()].append(f)
        else:
            featureDict[f.name()] = [f]

    # output is set to stdout or a file name if specified
    if outFile and outFile != "-":
        try:
            with open(outFile, "a" if append else "w", encoding="utf8") as fp:
                writeFEPFile(featureDict, inFile, fp, precision)
        except Exception as e:
            error(2, e)

    else:
        writeFEPFile(featureDict, inFile, sys.stdout, precision)
        sys.stdout.flush()


# -----------------------------------------------------------------------------
def writeFEPFile(featureDict: dict, inFile: str, fp: TextIO, precision: int = None):
    def _print(data: str):
        print(data, file=fp)

    if precision:

        def _printVertex(v):
            print(
                f"{v.longitude():.{precision}f} {v.latitude():.{precision}f}", file=fp
            )

    else:

        def _printVertex(v):
            print(f"{v.longitude()} {v.latitude()}", file=fp)

    _print(f"# created from file: {inFile}")
    _print(
        f"# created on {str(datetime.datetime.now())} by gfs2fep.py - (C) gempa GmbH"
    )
    _print("# LON LAT")

    # write fep
    for name, features in featureDict.items():
        # print("{}: {}".format(len(features), name))
        vCount = 0
        fStart = features[0].vertices()[0]
        v = fStart

        # iterate over features sharing name
        for f in features:
            # vertex array contains vertices of main land and sub features
            vertices = f.vertices()

            # sub feature array holds indices of starting points
            endIndices = list(f.subFeatures()) + [len(vertices)]

            # iterate of main land and sub features
            i = 0
            for iEnd in endIndices:
                vStart = vertices[i]
                while i < iEnd:
                    v = vertices[i]
                    _printVertex(v)
                    vCount += 1
                    i += 1

                # end sub feature on sub feature start
                v = vStart
                _printVertex(v)
                vCount += 1

                # go back to start of main land
                if v != vertices[0]:
                    v = vertices[0]
                    _printVertex(v)
                    vCount += 1

            # go back to start of first feature
            if v != fStart:
                v = fStart
                _printVertex(v)
                vCount += 1

        # end fep region
        _print(f"99.0 99.0 {vCount}")
        _print(f"L  {name}")


# -----------------------------------------------------------------------------
if __name__ == "__main__":
    run()