Get corners not working for segy_loader

[trhallam]

Get corners not working for segy_loader

[trhallam]

A couple of QC things you can try:

• Are there valid xy points the loaded xarray Dataset?
• Are the corners of your ds populated with values? If you loaded from a sparse segy data volume you may need to run the fill_cdp command first.

There is a good example of scraping and checking the headers here. https://segysak.readthedocs.io/en/latest/examples/example_segy_headers.html#Scraping-Headers

[ksurf1]

the corners are valid in il/xl but the corners_xy are all nan

... 'measurement_system': 'm', 'd3_domain': None, 'epsg': None, 'corner_points': ((1001, 3001), (1001, 4875), (2234, 4875), (2234, 3001)), 'corner_points_xy': ((nan, nan), (nan, nan), (nan, nan), (nan, nan)), 'source_file': 'path.segy', 'srd': None, 'datatype': None, 'percentiles': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 'coord_scalar': -100.0}

[trhallam]

X and Y corners are for the CDPX and CDPY values. IL and XL corners you can get directly from the dataset dimension vars.

If your X and Y values are not populating correctly it suggests to me that you might not have X and Y data in your SEGY or you have not selected the correct byte positions for the cdpx and cdpy variables.

Can you share the output of segy_header_scan(fn) to display the full output you might have to do this:

import pandas as pd
from IPython.display import display

with pd.option_context("display.max_rows", 89):
    display(segy_header_scan(fn))

[ksurf1]

scrape = scrape[scan[scan["std"] > 0].index]
scrape

	TRACE_SEQUENCE_FILE	TraceNumber	CDP	SourceX	SourceY	CDP_X	CDP_Y	INLINE_3D	CROSSLINE_3D
1001	1	3183	43518576	49289683	43518576	49289683	1001	3183
1001	2	3184	43519826	49289683	43519826	49289683	1001	3184
1001	3	3185	43521076	49289683	43521076	49289683	1001	3185
1001	4	3186	43522326	49289683	43522326	49289683	1001	3186
1001	5	3187	43523576	49289683	43523576	49289683	1001	3187
...	...	...	...	...	...	...	...	...
1025	9996	3548	43974826	49259683	43974826	49259683	1025	3548
1025	9997	3549	43976076	49259683	43976076	49259683	1025	3549
1025	9998	3550	43977326	49259683	43977326	49259683	1025	3550
1025	9999	3551	43978576	49259683	43978576	49259683	1025	3551
1025	10000	3552	43979826	49259683	43979826	49259683	1025	3552

[ksurf1]

A couple of QC things you can try:

• Are there valid xy points the loaded xarray Dataset?
• Are the corners of your ds populated with values? If you loaded from a sparse segy data volume you may need to run the fill_cdp command first.

There is a good example of scraping and checking the headers here. https://segysak.readthedocs.io/en/latest/examples/example_segy_headers.html#Scraping-Headers

ds.seis.fill_cdpna() did the trick

here are my corners after that:
((432910.76000091387, 492896.8299149062), (432910.75999001024, 477484.3295504351), (456335.7600132545, 492896.83054377494), (456335.76000235084, 477484.3301793039))

THANKS !

[trhallam]

You must have had a sparse volume with missing traces.

[ksurf1]

correct. That is the case for most data by the way. You should consider running this seis_fill_cdpna() by default !

[trhallam]

correct. That is the case for most data by the way. You should consider running this seis_fill_cdpna() by default !

Glad you got it sorted but I don't think this is appropriate, the default behaviour is just to load the data available in the segy file. fill_cdpna makes some assumptions to populate the empty traces with coordinates and that behaviour shouldn't be enforced on all users.

[ksurf1]

Ok, you may want to warn users when XYs are na in the xarray Dataset when it is converted to a ZGY.
Otherwise users will create ZGYs with nan corners without any warning.
Since fill_cdpna only alters XYs of null or missing traces, I don't see harm in regenerating the XYs since these traces are dead.
Anyway, thanks for your code, this is my goto to convert SEGY files to ZGY files.

Thanks again for your work. OpenVDS is next for me.

That might help others, I added this function to your _openzgy.py file (kinda quick and dirty):

def zgy_writer2(seisnc_dataset, filename, datatype="float", datarange=(-215, 215-1), dimension=None):

"""Write a seisnc dataset to ZGY file format.
Args:
    seisnc_dataset (xr.Dataset): This should be a seisnc dataset.
    filename (pathlib.Path/str): The filename to write to.
    datatype  (str): "float", "int16" or "int8"
    datarange (float,float): (-1500.0, 1500.0)
    dimension (str, optional): The dimension to write out, if
        None uses available dimenions twt first.
"""
assert seisnc_dataset.seis.is_3d()

# fill gaps in data
seisnc_dataset.seis.fill_cdpna()

seisnc_dataset.seis.calc_corner_points()
if seisnc_dataset.corner_points_xy:
    corners = tuple(seisnc_dataset.corner_points_xy[i] for i in [0, 3, 1, 2])
else:
    corners=[(0, 0), (0, 0), (0, 0), (0, 0)]

# print( corners )

dimension = _check_dimension(seisnc_dataset, dimension)
if dimension == "twt":
    zunitdim = UnitDimension(2001)
elif dimension == "depth":
    zunitdim = UnitDimension(2002)
else:
    zunitdim = UnitDimension(2000)
try:
    coord_scalar_mult = seisnc_dataset.attrs["coord_scalar_mult"]
except KeyError:
    coord_scalar_mult = 1.0

# dimensions
ni, nj, nk = (
    seisnc_dataset.dims[CoordKeyField.iline],
    seisnc_dataset.dims[CoordKeyField.xline],
    seisnc_dataset.dims[dimension],
)

# vertical
z0 = int(seisnc_dataset[dimension].values[0])
dz = int(seisnc_dataset.sample_rate)

# annotation
il0 = seisnc_dataset[CoordKeyField.iline].values[0]
xl0 = seisnc_dataset[CoordKeyField.xline].values[0]
dil = int(np.median(np.diff(seisnc_dataset[CoordKeyField.iline].values)))
dxl = int(np.median(np.diff(seisnc_dataset[CoordKeyField.xline].values)))

dtype_float = SampleDataType[datatype]

if datatype == "float":

    with ZgyWriter(
        str(filename),
        size=(ni, nj, nk),
        datatype=dtype_float,
        zunitdim=zunitdim,
        zstart=z0,
        zinc=dz,
        annotstart=(il0, xl0),
        annotinc=(dil, dxl),
        corners=corners,
    ) as writer:
        
        order = (CoordKeyField.iline, CoordKeyField.xline, dimension)

        writer.write(
            (0, 0, 0),
            seisnc_dataset[VariableKeyField.data]
            .transpose(*order)
            .values.astype(np.float32),
        )

else:
    
    with ZgyWriter(
        str(filename),
        size=(ni, nj, nk),
        datatype=dtype_float,
        datarange=datarange,
        zunitdim=zunitdim,
        zstart=z0,
        zinc=dz,
        annotstart=(il0, xl0),
        annotinc=(dil, dxl),
        corners=corners,
    ) as writer:

        order = (CoordKeyField.iline, CoordKeyField.xline, dimension)

        writer.write(
            (0, 0, 0),
            seisnc_dataset[VariableKeyField.data]
            .transpose(*order)
            .values.astype(np.float32),
        )