reproject_segy.py

Utility script to reproject/transform SEG-Y header coordinates between coordinate reference systems (CRS).

dms2dd(dms)

Convert DMS (Degrees, Minutes, Seconds) coordinate string to DD (Decimal Degrees).

Source code in pseudo_3D_interpolation\reproject_segy.py

def dms2dd(dms):
    """Convert DMS (Degrees, Minutes, Seconds) coordinate string to DD (Decimal Degrees)."""
    dms_str = str(dms)
    degrees, minutes, seconds = int(dms_str[:3]), int(dms_str[3:5]), int(dms_str[5:])
    if len(str(seconds)) > 2:
        seconds = float(str(seconds)[:2] + '.' + str(seconds)[2:])
    return float(degrees) + float(minutes) / 60 + float(seconds) / (60 * 60)

wrapper_reproject_segy(in_path, src_coords_bytes, dst_coords_bytes, args)

Reproject SEG-Y header coordinates for single SEG-Y file.

Parameters:

• in_path (str) –

  SEG-Y input path.

• src_coords_bytes (tuple) –

  Input byte position of coordinates in trace header.

• dst_coords_bytes (tuple) –

  Output byte position of coordinates in trace header.

• args (Namespace) –

  Input parameter.

Source code in pseudo_3D_interpolation\reproject_segy.py

def wrapper_reproject_segy(in_path, src_coords_bytes, dst_coords_bytes, args):
    """
    Reproject SEG-Y header coordinates for single SEG-Y file.

    Parameters
    ----------
    in_path : str
        SEG-Y input path.
    src_coords_bytes : tuple
        Input byte position of coordinates in trace header.
    dst_coords_bytes : tuple
        Output byte position of coordinates in trace header.
    args : argparse.Namespace
        Input parameter.

    """
    basepath, filename = os.path.split(in_path)
    basename, suffix = os.path.splitext(filename)
    xprint(f'Processing file < {filename} >', kind='info', verbosity=args.verbose)

    default_txt_suffix = 'reproj'

    if args.inplace is True:  # `inplace` parameter supersedes any `output_dir`
        xprint('Updating SEG-Y inplace', kind='warning', verbosity=args.verbose)
        path = in_path
    else:
        if args.output_dir is None:  # default behavior
            xprint('Creating copy of file in INPUT directory:\n', basepath, kind='info', verbosity=args.verbose)
            out_dir = basepath
        elif args.output_dir is not None and os.path.isdir(args.output_dir):
            xprint('Creating copy of file in OUTPUT directory:\n', args.output_dir, kind='info', verbosity=args.verbose)
            out_dir = args.output_dir
        else:
            raise FileNotFoundError(f'The output directory > {args.output_dir} < does not exist')

        if args.txt_suffix is not None:
            out_name = f'{basename}_{args.txt_suffix}'
        else:
            out_name = f'{basename}_{default_txt_suffix}'
        out_path = os.path.join(out_dir, f'{out_name}{suffix}')

        # sanity check
        if os.path.isfile(out_path):
            xprint('Output file already exists and will be removed!', kind='warning', verbosity=args.verbose)
            os.remove(out_path)

        copy2(in_path, out_path)
        path = out_path

    # read SEGY file
    with segyio.open(path, 'r+', strict=False, ignore_geometry=True) as file:
        # get scaled coordinates
        xcoords, ycoords, coordinate_units = scale_coordinates(file, src_coords_bytes)

        # get source & destination CRS
        crs_src = pyproj.crs.CRS(args.crs_src)
        if coordinate_units != 1 and crs_src.is_projected:
            xprint(
                'Forced source CRS to be geographic (WGS84 - EPSG:4326)!',
                kind='warning',
                verbosity=args.verbose,
            )
            crs_src = pyproj.crs.CRS('epsg:4326')
        crs_dst = pyproj.crs.CRS(args.crs_dst)

        # create CRS transformer
        transformer = pyproj.transformer.Transformer.from_crs(crs_src, crs_dst, always_xy=True)

        # convert coordinates
        xcoords_t, ycoords_t = transformer.transform(xcoords, ycoords, errcheck=True)

        # [OPTIONAL] smooth coordinates
        if args.smooth is not None:
            xcoords_t = smooth(xcoords_t, args.smooth)
            ycoords_t = smooth(ycoords_t, args.smooth)

        # update coordinates in trace header
        set_coordinates(
            file,
            xcoords_t,
            ycoords_t,
            crs_dst,
            dst_coords_bytes,
            coordinate_units=1,
            scaler=args.scalar_coords,
        )

    # update textual header
    text = get_textual_header(path)
    ## add info about new CRS
    info = f'EPSG:{crs_dst.to_epsg()}'
    text_updated = add_processing_info_header(text, info, prefix='CRS (PROJECTED)')
    # add info about byte position of reprojected coords
    info = f'REPROJECT (BYTES:{dst_coords_bytes[0]} {dst_coords_bytes[1]})'
    info += ' SMOOTHED' if args.smooth is not None else ''
    text_updated = add_processing_info_header(text_updated, info, prefix='_TODAY_')
    write_textual_header(path, text_updated)

main(argv=sys.argv)

Reproject trace header coordinats of SEG-Y file(s).

Source code in pseudo_3D_interpolation\reproject_segy.py

def main(argv=sys.argv):  # noqa
    """Reproject trace header coordinats of SEG-Y file(s)."""
    TIMESTAMP = datetime.datetime.now().isoformat(timespec='seconds').replace(':', '')
    SCRIPT = os.path.basename(__file__).split(".")[0]  # FIXME

    parser = define_input_args()
    args = parser.parse_args(argv[1:])  # exclude filename parameter at position 0
    xprint(args, kind='debug', verbosity=args.verbose)

    # input and output coordinate byte positions
    src_coords_bytes = TRACE_HEADER_COORDS.get(args.src_coords)
    dst_coords_bytes = TRACE_HEADER_COORDS.get(args.dst_coords)

    # check input file(s)
    in_path = args.input_path
    basepath, filename = os.path.split(in_path)
    basename, suffix = os.path.splitext(filename)
    if suffix == '':
        basepath = in_path
        basename, suffix = None, None  # noqa

    # (1) single input file
    if os.path.isfile(in_path) and (suffix != '.txt'):
        # perform coordinate transformation
        wrapper_reproject_segy(in_path, src_coords_bytes, dst_coords_bytes, args)
        sys.exit()

    # (2) input directory (multiple files)
    elif os.path.isdir(in_path):
        pattern = '*'
        pattern += f'{args.filename_suffix}' if args.filename_suffix is not None else pattern
        pattern += f'.{args.suffix}' if args.suffix is not None else '.sgy'
        file_list = glob.glob(os.path.join(in_path, pattern))

    # (3) file input is datalist (multiple files)
    elif os.path.isfile(in_path) and (suffix == '.txt'):
        with open(in_path, 'r') as datalist:
            file_list = datalist.readlines()
            file_list = [
                os.path.join(basepath, line.rstrip())
                if os.path.split(line.rstrip()) not in ['', '.']
                else line.rstrip()
                for line in file_list
            ]
    else:
        raise FileNotFoundError('Invalid input file')

    # compute files from options (2) or (3)
    if len(file_list) > 0:
        # redirect stdout to logfile
        logfile = os.path.join(basepath, f'{TIMESTAMP}_{SCRIPT}.log')
        with open(logfile, 'w', newline='\n') as f:
            with redirect_stdout(f):
                xprint(f'Processing total of < {len(file_list)} > files', kind='info', verbosity=args.verbose)
                for file_path in tqdm(
                    file_list,
                    desc='Reprojecting SEG-Y',
                    ncols=80,
                    total=len(file_list),
                    unit_scale=True,
                    unit=' files',
                ):
                    # perform coordinate transformation
                    wrapper_reproject_segy(file_path, src_coords_bytes, dst_coords_bytes, args)
        clean_log_file(logfile)
    else:
        sys.exit('No input files to process. Exit process.')

---

Last update: Monday, 03 July 2023 at 09:46:51