From fc37c4c67eb2e52486fc83d743ab4705dadd16b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20C=2E=20Riven=C3=A6s?= Date: Mon, 18 Dec 2023 18:11:57 +0100 Subject: [PATCH] ENH: improve i/o for IJXYZ format --- src/clib/xtg/libxtg.h | 30 -- src/clib/xtg/surf_import_ijxyz.c | 483 ------------------ src/clib/xtg/surf_import_ijxyz_tmpl.c | 141 ----- src/xtgeo/surface/_regsurf_export.py | 45 +- src/xtgeo/surface/_regsurf_ijxyz_parser.py | 267 ++++++++++ src/xtgeo/surface/_regsurf_import.py | 134 +---- src/xtgeo/surface/regular_surface.py | 10 +- tests/test_surface/test_file_io.py | 4 +- tests/test_surface/test_ijxyz_spec.py | 194 +++++++ .../test_regular_surface_bytesio.py | 27 + 10 files changed, 569 insertions(+), 766 deletions(-) delete mode 100644 src/clib/xtg/surf_import_ijxyz.c delete mode 100644 src/clib/xtg/surf_import_ijxyz_tmpl.c create mode 100644 src/xtgeo/surface/_regsurf_ijxyz_parser.py create mode 100644 tests/test_surface/test_ijxyz_spec.py diff --git a/src/clib/xtg/libxtg.h b/src/clib/xtg/libxtg.h index 8c5550fb2..6e6e623fa 100644 --- a/src/clib/xtg/libxtg.h +++ b/src/clib/xtg/libxtg.h @@ -264,36 +264,6 @@ surf_import_irap_ascii(FILE *fc, long n_swig_np_dbl_aout_v1, // nmap int option); -int -surf_import_ijxyz(FILE *fc, - int mode, - int *swig_int_out_p1, // *nx - int *swig_int_out_p2, // *ny - long *swig_lon_out_p1, // *ndef, - double *swig_dbl_out_p1, // *xori, - double *swig_dbl_out_p2, // *yori, - double *swig_dbl_out_p3, // *xinc, - double *swig_dbl_out_p4, // *yinc, - double *swig_dbl_out_p5, // *rot, - int *swig_np_int_aout_v1, // *ilines, - long n_swig_np_int_aout_v1, // nilines or ncol, - int *swig_np_int_aout_v2, // *xlines, - long n_swig_np_int_aout_v2, // nxlines or nrow, - double *swig_np_dbl_aout_v1, // *p_map_v - long n_swig_np_dbl_aout_v1, // nmap = nrow*ncol - int *swig_int_out_p3, // yflip - int option); - -int -surf_import_ijxyz_tmpl(FILE *fc, - int *swig_np_int_in_v1, // *ilines, - long n_swig_np_int_in_v1, // nilines or ncol, - int *swig_np_int_in_v2, // *xlines, - long n_swig_np_int_in_v2, // nxlines or nrow, - double *swig_np_dbl_aout_v1, // *p_map_v - long n_swig_np_dbl_aout_v1, // nmap = nrow*ncol - int option); - void surf_import_petromod_bin(FILE *fhandle, int mode, diff --git a/src/clib/xtg/surf_import_ijxyz.c b/src/clib/xtg/surf_import_ijxyz.c deleted file mode 100644 index 3cd76e893..000000000 --- a/src/clib/xtg/surf_import_ijxyz.c +++ /dev/null @@ -1,483 +0,0 @@ -/* -**************************************************************************************** -* -* Import maps on seismic points format (OW specific perhaps?) -* -*************************************************************************************** -*/ - -#include "libxtg.h" -#include "libxtg_.h" -#include "logger.h" -#include - -#define MAXIX 1000000 - -/* -**************************************************************************************** - * - * NAME: - * surf_import_ijxyz.c - * - * DESCRIPTION: - * Import a map on DSG I J XYZ format. This format is a coordinate value - * per cell. Undefined cells are simply not mentioned. The routine is - * a bit complicated, and numerous "options" may occur. Example: - * 1347 2020 532398.6557089972 6750702.972141342 11712.976 - * 1340 2049 532758.1090287864 6750842.3494704 7387.8145 - * 1340 2048 532748.6891376793 6750834.13270089 7485.6885 - * 1347 2029 532483.4347289609 6750776.923066932 7710.251 - * 1347 2028 532474.0148378538 6750768.706297422 7850.605 - * 1347 2027 532464.5949467467 6750760.489527912 7640.044 - * 1347 2026 532455.1750556396 6750752.272758402 7700.624 - * 1347 2025 532445.7551645326 6750744.055988892 8482.951 - * - * In some cases headersand footers are present. Assume footer when - * line is like this, and I J may be float, not INT - * @File_Version: 4 - @Coordinate_Type_is: 3 - @Export_Type_is: 1 - @Number_of_Projects 1 - @Project_Type_Name: , 3,TROLL_WEST, - @Project_Unit_is: meters , ST_ED50_UTMXXX , PROJECTED_COORDINATE_SYSTEM - #File_Version____________-> 4 - #Project_Name____________-> TROLL_WEST - #Horizon_remark_size_____-> 0 - - #End_of_Horizon_ASCII_Header_ - - - * - * See there is no system in fastest vs slowest ILINE vs XLINE! - * - * ARGUMENTS: - * filename i File name, character string - * mode i If 0, then scan only for dimensions - * mx i/o Map dimension X (I) - * my i/o Map dimension Y (J) - * xori o X origin coordinate - * yori o Y origin coordinate - * xinc o X increment - * yinc o Y increment - * rot o Rotation (degrees, from X axis, anti-clock) - * ilines o Vector of inline numbers - * xlines o Vector of xline numbers - * p_surf_v o 1D pointer to map/surface values pointer array - * yflip o YFLIP indicator (1 for normal, -1 for reverse) - * option i For later use - * debug i Debug level - * - * RETURNS: - * Function: 0: upon success. If problems <> 0: - * - * TODO/ISSUES/BUGS: - * Code could benefit from a refactorization - * - * LICENCE: - * See XTGeo license - *************************************************************************************** - */ - -/* -**************************************************************************************** -* Local routine to scan for NX, NY -*************************************************************************************** -*/ - -void -_scan_dimensions(FILE *fd, int *nx, int *ny) -{ - int inum, nrow, ncol, iline, xline, iok; - int ilinemin, ilinemax, xlinemin, xlinemax; - int ispacing, xspacing, ispace, mispace, mxspace; - int *itmp, *xtmp; - float rdum, filine, fxline; - char lbuffer[132] = ""; - - itmp = calloc(MAXIX, 4); - xtmp = calloc(MAXIX, 4); - - nrow = 0; - ncol = 0; - - for (inum = 0; inum < MAXIX; inum++) - itmp[inum] = 0; - for (inum = 0; inum < MAXIX; inum++) - xtmp[inum] = 0; - - ilinemin = 999999999; - ilinemax = -99999999; - xlinemin = 999999999; - xlinemax = -99999999; - iok = 1; - while (fgets(lbuffer, 132, (FILE *)fd)) { - if (strncmp(lbuffer, "\n", 1) == 0) - continue; - lbuffer[strcspn(lbuffer, "\n")] = 0; - if (strncmp(lbuffer, "#", 1) == 0) - continue; - if (strncmp(lbuffer, "@", 1) == 0) - continue; - if (strncmp(lbuffer, "E", 1) == 0) - continue; - - iok = sscanf(lbuffer, "%f %f %f %f %f", &filine, &fxline, &rdum, &rdum, &rdum); - - iline = (int)(filine + 0.01); - xline = (int)(fxline + 0.01); - - if (iok > 5) - logger_error(LI, FI, FU, "Wrong file format for map file?"); - - if (iline < ilinemin) - ilinemin = iline; - if (iline > ilinemax) - ilinemax = iline; - if (xline < xlinemin) - xlinemin = xline; - if (xline > xlinemax) - xlinemax = xline; - itmp[iline] = 1; - xtmp[xline] = 1; - } - - mispace = (ilinemax - ilinemin) / 4; - mxspace = (xlinemax - xlinemin) / 4; - - /* find minimum spacing in INLINES */ - ispacing = 0; - for (ispace = 1; ispace < mispace; ispace++) { - for (inum = ilinemin; inum < ilinemax - mispace; inum++) { - if (itmp[inum] == 1 && itmp[inum + ispace] == 1) { - ispacing = ispace; - break; - } - } - if (ispacing > 0) - break; - } - - /* find minimum spacing in XLINES */ - xspacing = 0; - for (ispace = 1; ispace < mxspace; ispace++) { - for (inum = xlinemin; inum < xlinemax - mxspace; inum++) { - if (xtmp[inum] == 1 && xtmp[inum + ispace] == 1) { - xspacing = ispace; - break; - } - } - if (xspacing > 0) - break; - } - - *nx = (ilinemax - ilinemin) / ispacing + 1; - *ny = (xlinemax - xlinemin) / xspacing + 1; - - free(itmp); - free(xtmp); -} - -/* -**************************************************************************************** -* Local routine to collect all values form file as 1D buffers -*************************************************************************************** -*/ - -long -_collect_values(FILE *fd, - int *ilinesb, - int *xlinesb, - double *xbuffer, - double *ybuffer, - double *zbuffer, - int *ilmin, - int *ilmax, - int *xlmin, - int *xlmax) -{ - int ilinemin, ilinemax, xlinemin, xlinemax, iok; - int iline, xline; - float filine, fxline; - long ixnum; - double xval, yval, zval; - char lbuffer[132]; - - ilinemin = 999999999; - ilinemax = -99999999; - xlinemin = 999999999; - xlinemax = -99999999; - - ixnum = 0; - iok = 1; - while (fgets(lbuffer, 132, (FILE *)fd)) { - if (strncmp(lbuffer, "\n", 1) == 0) - continue; - lbuffer[strcspn(lbuffer, "\n")] = 0; - if (strncmp(lbuffer, "#", 1) == 0) - continue; - if (strncmp(lbuffer, "@", 1) == 0) - continue; - if (strncmp(lbuffer, "E", 1) == 0) - continue; - - iok = - sscanf(lbuffer, "%f %f %lf %lf %lf", &filine, &fxline, &xval, &yval, &zval); - - iline = (int)(filine + 0.01); - xline = (int)(fxline + 0.01); - - ilinesb[ixnum] = iline; - xlinesb[ixnum] = xline; - xbuffer[ixnum] = xval; - ybuffer[ixnum] = yval; - zbuffer[ixnum] = zval; - - if (iline < ilinemin) - ilinemin = iline; - if (iline > ilinemax) - ilinemax = iline; - if (xline < xlinemin) - xlinemin = xline; - if (xline > xlinemax) - xlinemax = xline; - - ixnum++; - } - - *ilmin = ilinemin; - *ilmax = ilinemax; - *xlmin = xlinemin; - *xlmax = xlinemax; - - return ixnum; -} - -/* -**************************************************************************************** -* Compute map vectors -*/ - -int -_compute_map_vectors(int ilinemin, - int ilinemax, - int ncol, - int xlinemin, - int xlinemax, - int nrow, - long ixnummax, - int *ilinesb, - int *xlinesb, - double *xbuffer, - double *ybuffer, - double *zbuffer, - double *xcoord, - double *ycoord, - /* result vectors: */ - int *ilines, - int *xlines, - double *p_map_v) -{ - int i, itrue, jtrue; - long ic, ixnum; - int ilinestep, xlinestep; - - ilinestep = (ilinemax - ilinemin) / (ncol - 1); - xlinestep = (xlinemax - xlinemin) / (nrow - 1); - - /* set the *ilines and *xlines results vectors */ - for (i = 0; i < ncol; i++) - ilines[i] = ilinemin + i * ilinestep; - for (i = 0; i < nrow; i++) - xlines[i] = xlinemin + i * xlinestep; - - /* setting Z values; UNDEF initially */ - for (ic = 0; ic < ncol * nrow; ic++) - p_map_v[ic] = UNDEF; - - for (ixnum = 0; ixnum < ixnummax; ixnum++) { - itrue = (ilinesb[ixnum] / ilinestep) - ilinemin / ilinestep + 1; /* base = 1 */ - jtrue = (xlinesb[ixnum] / xlinestep) - xlinemin / xlinestep + 1; /* base = 1 */ - - /* get the C order index */ - ic = x_ijk2ic(itrue, jtrue, 1, ncol, nrow, 1, 0); - if (ic < 0) { - throw_exception("Index outside boundary in surf_import_ijxyz"); - return EXIT_FAILURE; - } - - p_map_v[ic] = zbuffer[ixnum]; - xcoord[ic] = xbuffer[ixnum]; - ycoord[ic] = ybuffer[ixnum]; - } - return EXIT_SUCCESS; -} - -/* -**************************************************************************************** -* Compute map properties in XTGeo format, need to deduce from "incomplete -* data". This version is rather basic, it looks only on one step. It may be -* more precise to look at long vectors, but that is more complicated to -* get correct. -*************************************************************************************** -*/ - -int -_compute_map_props(int ncol, - int nrow, - double *xcoord, - double *ycoord, - double *p_map_v, - double *xori, - double *yori, - double *xinc, - double *yinc, - double *rot, - int *yflip) -{ - - int okstatus = 0, icol, jrow; - long icn0, icn1, icn2, ic0 = 0, jc0 = 1; - double xc0 = 0.0, xc1 = 0.0, xc2 = 0.0, yc0 = 0.0, yc1 = 0.0, yc2 = 0.0; - double a1rad, a2rad, roty, sinusrad; - - okstatus = 0; - - for (icol = 1; icol < ncol; icol++) { - for (jrow = 1; jrow < nrow; jrow++) { - - icn0 = x_ijk2ic(icol, jrow, 1, ncol, nrow, 1, 0); - icn1 = x_ijk2ic(icol + 1, jrow, 1, ncol, nrow, 1, 0); - icn2 = x_ijk2ic(icol, jrow + 1, 1, ncol, nrow, 1, 0); - if (icn0 < 0 || icn1 < 0 || icn2 < 0) { - throw_exception("Loop through surface gave index outside boundary in " - "surf_import_ijxyz"); - return EXIT_FAILURE; - } - - if (p_map_v[icn0] < UNDEF_LIMIT && p_map_v[icn1] < UNDEF_LIMIT && - p_map_v[icn2] < UNDEF_LIMIT) { - - xc0 = xcoord[icn0]; - xc1 = xcoord[icn1]; - xc2 = xcoord[icn2]; - - yc0 = ycoord[icn0]; - yc1 = ycoord[icn1]; - yc2 = ycoord[icn2]; - - ic0 = icol; - jc0 = jrow; - - okstatus = 1; - break; - } - if (okstatus == 1) - break; - } - if (okstatus == 1) - break; - } - - if (okstatus == 0) { - logger_error(LI, FI, FU, "Could not find info to deduce map properties"); - return -9; - } - - x_vector_info2(xc0, xc1, yc0, yc1, xinc, &a1rad, rot, 1); - x_vector_info2(xc0, xc2, yc0, yc2, yinc, &a2rad, &roty, 1); - - /* compute yflip: sin (y-x) = sin(y)*cos(x) - sin(x)*cos(y) */ - *yflip = 1; - sinusrad = sin(a2rad) * cos(a1rad) - sin(a1rad) * cos(a2rad); - if (sinusrad < 0) - *yflip = -1; - - /* find origin */ - surf_xyori_from_ij(ic0, jc0, xc0, yc0, xori, *xinc, yori, *yinc, ncol, nrow, *yflip, - *rot, 0); - - return EXIT_SUCCESS; -} - -int -surf_import_ijxyz(FILE *fd, - int mode, - int *nx, - int *ny, - long *ndef, - double *xori, - double *yori, - double *xinc, - double *yinc, - double *rot, - int *ilines, - long ncol, - int *xlines, - long nrow, - double *p_map_v, - long nmap, - int *yflip, - int option) -{ - - /* locals*/ - int iok = 0; - int ilinemin, ilinemax, xlinemin, xlinemax; - - double *xbuffer, *ybuffer, *zbuffer, *xcoord, *ycoord; - int *ilinesb, *xlinesb; - int nncol, nnrow; - - /* read header */ - logger_info(LI, FI, FU, "Entering routine %s", __FUNCTION__); - - fseek(fd, 0, SEEK_SET); - - /* ================================================================================= - * scan mode; to determine dimensions */ - if (mode == 0) { - _scan_dimensions(fd, nx, ny); - - return EXIT_SUCCESS; - } - - /* ================================================================================= - * Read mode; now dimensions shall be known */ - - nncol = (int)ncol; - nnrow = (int)nrow; - - *nx = nncol; - *ny = nnrow; - - ilinesb = calloc(ncol * nrow + 10, sizeof(int)); - xlinesb = calloc(ncol * nrow + 10, sizeof(int)); - xbuffer = calloc(ncol * nrow + 10, sizeof(double)); - ybuffer = calloc(ncol * nrow + 10, sizeof(double)); - zbuffer = calloc(ncol * nrow + 10, sizeof(double)); - xcoord = calloc(ncol * nrow + 10, sizeof(double)); - ycoord = calloc(ncol * nrow + 10, sizeof(double)); - - *ndef = _collect_values(fd, ilinesb, xlinesb, xbuffer, ybuffer, zbuffer, &ilinemin, - &ilinemax, &xlinemin, &xlinemax); - - iok = _compute_map_vectors(ilinemin, ilinemax, nncol, xlinemin, xlinemax, nnrow, - *ndef, ilinesb, xlinesb, xbuffer, ybuffer, zbuffer, - xcoord, ycoord, ilines, xlines, p_map_v); - - iok = _compute_map_props(nncol, nnrow, xcoord, ycoord, p_map_v, xori, yori, xinc, - yinc, rot, yflip); - - if (iok != 0) - logger_error(LI, FI, FU, "Error, cannot compute map props"); - - free(ilinesb); - free(xlinesb); - free(xbuffer); - free(ybuffer); - free(zbuffer); - free(xcoord); - free(ycoord); - - return EXIT_SUCCESS; -} diff --git a/src/clib/xtg/surf_import_ijxyz_tmpl.c b/src/clib/xtg/surf_import_ijxyz_tmpl.c deleted file mode 100644 index 90f8c25e9..000000000 --- a/src/clib/xtg/surf_import_ijxyz_tmpl.c +++ /dev/null @@ -1,141 +0,0 @@ -/* -**************************************************************************************** -* -* Import map value on seismic points format (OW specific perhaps?) on the -* premise that the map topology (through INLINE and XLINE is known); hence -* map values are based on INLINE/XLINE match. -* -*************************************************************************************** -*/ - -#include "libxtg.h" -#include "libxtg_.h" -#include "logger.h" - -/* -**************************************************************************************** -* -* NAME: -* surf_import_ijxyz_tmpl.c -* -* AUTHOR(S): -* -* -* DESCRIPTION: -* Import a map on DSG I J XYZ format. This format is a coordinate value -* per cell. Undefined cells are simply not mentioned. The routine is -* a bit complicated, and numerous "options" may occur. Example: -* 1347 2020 532398.6557089972 6750702.972141342 11712.976 -* 1340 2049 532758.1090287864 6750842.3494704 7387.8145 -* 1340 2048 532748.6891376793 6750834.13270089 7485.6885 -* 1347 2029 532483.4347289609 6750776.923066932 7710.251 -* 1347 2028 532474.0148378538 6750768.706297422 7850.605 -* 1347 2027 532464.5949467467 6750760.489527912 7640.044 -* 1347 2026 532455.1750556396 6750752.272758402 7700.624 -* 1347 2025 532445.7551645326 6750744.055988892 8482.951 -* -* In some cases headersand footers are present. Assume footer when -* line is like this, and I J may be float, not INT -* @File_Version: 4 -* @Coordinate_Type_is: 3 -* @Export_Type_is: 1 -* @Number_of_Projects 1 -* @Project_Type_Name: , 3,TROLL_WEST, -* @Project_Unit_is: meters , ST_ED50_UTM31N_P23031_T1133 , PROJECTED_.. -* #File_Version____________-> 4 -* #Project_Name____________-> TROLL_WEST -* #Horizon_remark_size_____-> 0 -* -* #End_of_Horizon_ASCII_Header_ -* -* See there is no system in fastest vs slowest ILINE vs XLINE! -* -* ARGUMENTS: -* filename i File name, character string -* ilines i Vector of inline numbers, as template -* xlines i Vector of xline numbers, as template -* p_surf_v o 1D pointer to map/surface values pointer array -* option i For later use -* debug i Debug level -* -* RETURNS: -* Function: 0: upon success. If problems <> 0: -* -1: inline/xline in file is outside input ilines, xlines -* -* TODO/ISSUES/BUGS: -* Code could benefit from a refactorization -* -* LICENCE: -* See XTGeo license -*************************************************************************************** -*/ - -int -surf_import_ijxyz_tmpl(FILE *fd, - int *ilines, - long ncol, - int *xlines, - long nrow, - double *p_map_v, - long nmap, - int option) -{ - - /* locals*/ - - int iok, iline, xline, nncol, nnrow, found, ili, xli; - long ind, ic; - float rdum, zval, filine, fxline; - char lbuffer[132]; - - nncol = (int)ncol; - nnrow = (int)nrow; - - for (ind = 0; ind < nnrow * nncol; ind++) - p_map_v[ind] = UNDEF; - - while (fgets(lbuffer, 132, (FILE *)fd)) { - if (strncmp(lbuffer, "\n", 1) == 0) - continue; - lbuffer[strcspn(lbuffer, "\n")] = 0; - - if (strncmp(lbuffer, "#", 1) == 0) - continue; - if (strncmp(lbuffer, "@", 1) == 0) - continue; - if (strncmp(lbuffer, "E", 1) == 0) - continue; - - iok = sscanf(lbuffer, "%f %f %f %f %f", &filine, &fxline, &rdum, &rdum, &zval); - - iline = (int)(filine + 0.01); - xline = (int)(fxline + 0.01); - - /* some sanity tests first */ - if (iline < ilines[0] || iline > ilines[nncol - 1] || xline < xlines[0] || - xline > xlines[nnrow - 1]) { - return -1; // handle this code as error in client - } - - found = 0; - for (ili = 0; ili < nncol; ili++) { - for (xli = 0; xli < nnrow; xli++) { - if (iline == ilines[ili] && xline == xlines[xli]) { - ic = x_ijk2ic(ili + 1, xli + 1, 1, nncol, nnrow, 1, 0); - if (ic < 0) { - throw_exception("Loop through surface gave index outside " - "boundary in surf_import_ijxyz_tmpl"); - return EXIT_FAILURE; - } - p_map_v[ic] = zval; - found = 1; - break; - } - } - if (found == 1) - break; - } - } - - return EXIT_SUCCESS; -} diff --git a/src/xtgeo/surface/_regsurf_export.py b/src/xtgeo/surface/_regsurf_export.py index bec365bb8..451838303 100644 --- a/src/xtgeo/surface/_regsurf_export.py +++ b/src/xtgeo/surface/_regsurf_export.py @@ -1,7 +1,9 @@ """Export RegularSurface data.""" +from __future__ import annotations import json import struct +from typing import TYPE_CHECKING import h5py import hdf5plugin @@ -12,6 +14,10 @@ from xtgeo.common import null_logger from xtgeo.common.constants import UNDEF_MAP_IRAPA, UNDEF_MAP_IRAPB +if TYPE_CHECKING: + from xtgeo.common.sys import _XTGeoFile + from xtgeo.surface.regular_surface import RegularSurface + logger = null_logger(__name__) @@ -204,8 +210,43 @@ def _export_irap_binary_cxtgeo(self, mfile): mfile.cfclose() -def export_ijxyz_ascii(self, mfile): - """Export to DSG IJXYZ ascii format.""" +def export_ijxyz_ascii(self: RegularSurface, mfile: _XTGeoFile) -> None: + # cxtgeo is default since twice as fas as python, but use python if memstreams + if mfile.memstream: + _export_ijxyz_ascii_python(self, mfile) + else: + _export_ijxyz_ascii_cxtgeo(self, mfile) + + +def _export_ijxyz_ascii_python(self: RegularSurface, mfile: _XTGeoFile) -> None: + """Export to DSG IJXYZ ascii format, using python.""" + # the python version is ~twice as slow as the cxtgeo version + + dfr = self.get_dataframe(ij=True, order="F") # order F since this was in cxtgeo + ix = dfr.IX - 1 # since dataframe indexing starts at 1 + jy = dfr.JY - 1 + inlines = self.ilines[ix] + xlines = self.xlines[jy] + dfr["IL"] = inlines + dfr["XL"] = xlines + dfr = dfr[["IL", "XL", "X_UTME", "Y_UTMN", "VALUES"]] + + fmt = "%f" + + if mfile.memstream: + buf = dfr.to_csv(sep="\t", float_format=fmt, index=False, header=False) + buf = buf.encode("latin1") + mfile.file.write(buf) + else: + dfr.to_csv(mfile.name, sep="\t", float_format=fmt, index=False, header=False) + + +def _export_ijxyz_ascii_cxtgeo(self: RegularSurface, mfile: _XTGeoFile) -> None: + """Export to DSG IJXYZ ascii format, using cxtgeo. + + Keep this for while since it was the original solution, and faster, + but consider remove in e.g. 4.0.""" + vals = self.get_values1d(fill_value=xtgeo.UNDEF) ier = _cxtgeo.surf_export_ijxyz( mfile.get_cfhandle(), diff --git a/src/xtgeo/surface/_regsurf_ijxyz_parser.py b/src/xtgeo/surface/_regsurf_ijxyz_parser.py new file mode 100644 index 000000000..44984aab0 --- /dev/null +++ b/src/xtgeo/surface/_regsurf_ijxyz_parser.py @@ -0,0 +1,267 @@ +"""IJXYZ files (OW/DSG) parsing. + +This format lacks explicit information such as origin, xinc, yinc, rotation etc; i.e. +all these setting must be computed from the input. + +Missing inline and xline values are assumed to be undefined. Examples: + + +======================================================================================== +815 1210 777574.220000 6736507.760000 1010.000000 +816 1210 777561.894910 6736521.890000 1010.000000 +817 1210 777549.569820 6736536.020000 1010.000000 +818 1210 777537.244731 6736550.150000 1010.000000 +819 1210 777524.919641 6736564.280000 1010.000000 +820 1210 777512.594551 6736578.410000 1010.000000 +821 1210 777500.269461 6736592.540000 1010.000000 +822 1210 777487.944371 6736606.670000 1010.000000 +823 1210 777475.619281 6736620.800000 1010.000000 +... +846 1211 777201.561266 6736954.005898 1010.000000 +847 1211 777189.236176 6736968.135898 1010.000000 +848 1211 777176.911086 6736982.265898 1010.000000 +849 1211 777164.585996 6736996.395898 1010.000000 +850 1211 777152.260906 6737010.525898 1010.000000 +851 1211 777139.935817 6737024.655898 1010.000000 +... + +Here inline is fastest and increasing and every inline and xline is applied + +======================================================================================== + +24 3 774906.625000 6736030.500000 1741.107300 +26 3 774946.017310 6736023.554073 1730.530884 +28 3 774985.409620 6736016.608146 1719.947876 +30 3 775024.801930 6736009.662219 1709.229370 +32 3 775064.194240 6736002.716292 1698.470337 +34 3 775103.586551 6735995.770364 1687.320435 +... +272 3 779791.271455 6735169.205039 1404.342041 +274 3 779830.663765 6735162.259112 1413.154175 +24 6 774915.307409 6736079.740388 1727.455078 +26 6 774954.699719 6736072.794461 1719.031128 +28 6 774994.092029 6736065.848533 1710.606201 +30 6 775033.484339 6736058.902606 1702.161255 +32 6 775072.876649 6736051.956679 1693.707153 +34 6 775112.268959 6736045.010752 1685.164673 + +Here every second inline and every third xline is applied + +======================================================================================== + + 1312 2015 655905.938642 6627044.343360 1671.420800 + 1312 2010 655858.839186 6627003.259513 1669.091400 + 1302 2028 656151.649337 6627009.862811 1680.437500 + 1312 2011 655868.259077 6627011.476282 1669.444000 + 1302 2029 656161.069228 6627018.079581 1681.000000 + 1302 2026 656132.809555 6626993.429272 1679.312500 + 1302 2027 656142.229446 6627001.646042 1679.875000 + 1302 2024 656113.969773 6626976.995733 1678.187500 + 1302 2025 656123.389664 6626985.212503 1678.750000 + +This seems to lack obvious patterns (mostly due to many undefied cells?). Hence we +need to detect minimum spacing and min/max of both ilines and xlines, unless a template +is applied. +""" +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Literal, Optional + +import numpy as np +from numpy.ma import MaskedArray + +import xtgeo +from xtgeo.common.calc import find_flip, xyori_from_ij +from xtgeo.common.log import null_logger + +logger = null_logger(__name__) + + +@dataclass +class SurfaceIJXYZ: + """Temporary class for parsing IJXYZ settings. + + Attributes: + xcoords_in: Input X coordinates. + ycoords_in: Input Y coordinates. + values_in: Input values. + ilines_in: Input inline values. + xlines_in: Input crossline values. + template: Geometrical template. + """ + + xcoords_in: MaskedArray + ycoords_in: MaskedArray + values_in: MaskedArray + ilines_in: np.ndarray + xlines_in: np.ndarray + template: Optional[xtgeo.RegularSurface | xtgeo.Cube] = None + + ncol: int = field(default=1, init=False) + nrow: int = field(default=1, init=False) + xori: float = field(default=0.0, init=False) + yori: float = field(default=0.0, init=False) + xinc: float = field(default=1.0, init=False) + yinc: float = field(default=1.0, init=False) + yflip: Literal[-1, 1] = field(default=1, init=False) + xcoords: MaskedArray = field(default_factory=MaskedArray, init=False) + ycoords: MaskedArray = field(default_factory=MaskedArray, init=False) + values: MaskedArray = field(default_factory=MaskedArray, init=False) + ilines: np.ndarray = field(default_factory=lambda: np.zeros(1), init=False) + xlines: np.ndarray = field(default_factory=lambda: np.zeros(1), init=False) + + def __post_init__(self) -> None: + self._parse_arrays() + if self.template: + self._map_on_template() + else: + self._calculate_geometrics() + + def _parse_arrays(self) -> None: + """Parse input and arrange basic data and arrays.""" + + logger.debug("Parse ijxyz arrays") + # both ilines and xlines may jump e.g. every second; detect minimum jump: + inline_mindiff = int(np.min(np.diff(np.unique(np.sort(self.ilines_in))))) + xline_mindiff = int(np.min(np.diff(np.unique(np.sort(self.xlines_in))))) + + i_index = np.divide(self.ilines_in, inline_mindiff) + i_index = (i_index - i_index.min()).astype("int32") + + j_index = np.divide(self.xlines_in, xline_mindiff) + j_index = (j_index - j_index.min()).astype("int32") + + self.ncol = int(i_index.max() + 1) + self.nrow = int(j_index.max() + 1) + + self.ilines = np.array( + range( + self.ilines_in.min(), + self.ilines_in.max() + inline_mindiff, + inline_mindiff, + ) + ).astype("int32") + + self.xlines = np.array( + range( + self.xlines_in.min(), + self.xlines_in.max() + xline_mindiff, + xline_mindiff, + ) + ).astype("int32") + + shape = (self.ncol, self.nrow) + + xvalues = np.full(shape, np.nan) + xvalues[i_index, j_index] = self.xcoords_in + self.xcoords = np.ma.masked_invalid(xvalues) + + yvalues = np.full(shape, np.nan) + yvalues[i_index, j_index] = self.ycoords_in + self.ycoords = np.ma.masked_invalid(yvalues) + + zvalues = np.full(shape, np.nan) + zvalues[i_index, j_index] = self.values_in + self.values = np.ma.masked_invalid(zvalues) + + def _map_on_template(self) -> None: + """An existing RegularSurface or Cube forms the geometrical template.""" + logger.debug("Parse ijxyz with template") + if not isinstance(self.template, (xtgeo.RegularSurface, xtgeo.Cube)): + raise ValueError( + "The provided template is not of type RegularSurface or Cube" + ) + + ilines_orig = self.ilines.copy() + xlines_orig = self.xlines.copy() + + il_step = int(np.diff(np.sort(ilines_orig)).min()) + xl_step = int(np.diff(np.sort(xlines_orig)).min()) + + # now ovveride with settings from template + for attr in ( + "ncol", + "nrow", + "xori", + "yori", + "xinc", + "yinc", + "rotation", + "yflip", + "ilines", + "xlines", + ): + setattr(self, attr, getattr(self.template, attr)) + + i_start = abs(int((self.ilines.min() - ilines_orig.min()) / il_step)) + j_start = abs(int((self.xlines.min() - xlines_orig.min()) / xl_step)) + + shape = (self.ncol, self.nrow) + + zvalues = np.full(shape, np.nan) + zvalues[ + i_start : i_start + self.values.shape[0], + j_start : j_start + self.values.shape[1], + ] = np.ma.filled(self.values, fill_value=np.nan) + + self.values = np.ma.masked_invalid(zvalues) + + def _calculate_geometrics(self) -> None: + """Compute the geometrics such as rotation, flip origin and increments.""" + + logger.debug("Compute ijxyz geometrics from arrays") + # along_east and along_north here refers to the case when rotation is 0.0 + dx_along_east = np.diff(self.xcoords, axis=0) + dy_along_east = np.diff(self.ycoords, axis=0) + dx_along_north = np.diff(self.xcoords, axis=1) + dy_along_north = np.diff(self.ycoords, axis=1) + + self.yflip = find_flip( + (dx_along_east.mean(), dy_along_east.mean(), 0.0), + (dx_along_north.mean(), dy_along_north.mean(), 0.0), + ) + + self.rotation = float( + np.degrees(np.arctan2(dy_along_east, dx_along_east)).mean() + ) + + distances_xinc = np.sqrt(dx_along_east**2 + dy_along_east**2) + distances_yinc = np.sqrt(dx_along_north**2 + dy_along_north**2) + self.xinc = float(distances_xinc.mean()) + self.yinc = float(distances_yinc.mean()) + + cols, rows = np.ma.where(~self.xcoords.mask) # get active indices + + first_active_col = int(cols[0]) + first_active_row = int(rows[0]) + + self.xori, self.yori = xyori_from_ij( + first_active_col, + first_active_row, + self.xcoords[first_active_col, first_active_row], + self.ycoords[first_active_col, first_active_row], + self.xinc, + self.yinc, + self.ncol, + self.nrow, + self.yflip, + float(self.rotation), + ) + + +def parse_ijxyz(mfile, template=None) -> tuple: + """Read IJXYZ data file and parse the content.""" + + if mfile.memstream: + mfile.file.seek(0) + + data = np.loadtxt(mfile.file, comments=["@", "#", "EOB"]) + + inline = data[:, 0].astype("int32") + xline = data[:, 1].astype("int32") + x_arr = data[:, 2].astype("float64") + y_arr = data[:, 3].astype("float64") + z_arr = data[:, 4].astype("float64") + + return SurfaceIJXYZ(x_arr, y_arr, z_arr, inline, xline, template=template) diff --git a/src/xtgeo/surface/_regsurf_import.py b/src/xtgeo/surface/_regsurf_import.py index 0a7dac9b1..da48b4e0c 100644 --- a/src/xtgeo/surface/_regsurf_import.py +++ b/src/xtgeo/surface/_regsurf_import.py @@ -1,7 +1,9 @@ """Import RegularSurface data.""" +from __future__ import annotations import json from struct import unpack +from typing import TYPE_CHECKING import h5py import numpy as np @@ -12,8 +14,13 @@ from xtgeo import _cxtgeo from xtgeo.common import XTGeoDialog, null_logger from xtgeo.common.constants import UNDEF_MAP_IRAPA, UNDEF_MAP_IRAPB +from xtgeo.surface._regsurf_ijxyz_parser import parse_ijxyz from xtgeo.surface._zmap_parser import parse_zmap +if TYPE_CHECKING: + from xtgeo.cube.cube1 import Cube + from xtgeo.surface.regular_surface import RegularSurface + xtg = XTGeoDialog() logger = null_logger(__name__) @@ -250,114 +257,27 @@ def _import_irap_ascii(mfile): return args -def import_ijxyz(mfile, template=None, **_): - """Import OW/DSG IJXYZ ascii format.""" - return _import_ijxyz_tmpl(mfile, template) if template else _import_ijxyz(mfile) - - -def _import_ijxyz(mfile): +def import_ijxyz( + mfile: xtgeo._XTGeoFile, + template: RegularSurface | Cube | None = None, + **_, +) -> dict: """Import OW/DSG IJXYZ ascii format.""" - # import of seismic column system on the form: - # 2588 1179 476782.2897888889 6564025.6954 1000.0 - # 2588 1180 476776.7181777778 6564014.5058 1000.0 - logger.debug("Read data from file... (scan for dimensions)") - - cfhandle = mfile.get_cfhandle() - - xlist = _cxtgeo.surf_import_ijxyz(cfhandle, 0, 1, 1, 1, 0) - - ( - ier, - ncol, - nrow, - _, - xori, - yori, - xinc, - yinc, - rot, - iln, - xln, - val, - yflip, - ) = xlist - - if ier != 0: - mfile.cfclose() - raise RuntimeError("Import from C is wrong...") - - # now real read mode - xlist = _cxtgeo.surf_import_ijxyz(cfhandle, 1, ncol, nrow, ncol * nrow, 0) - - ier, ncol, nrow, _, xori, yori, xinc, yinc, rot, iln, xln, val, yflip = xlist - - if ier != 0: - raise RuntimeError("Import from C is wrong...") - - logger.info(xlist) - - val = ma.masked_greater(val, xtgeo.UNDEF_LIMIT) - args = {} - args["xori"] = xori - args["xinc"] = xinc - args["yori"] = yori - args["yinc"] = yinc - args["ncol"] = ncol - args["nrow"] = nrow - args["rotation"] = rot - args["yflip"] = yflip - - args["values"] = val.reshape((args["ncol"], args["nrow"])) - - args["ilines"] = iln - args["xlines"] = xln - - mfile.cfclose() - return args - - -def _import_ijxyz_tmpl(mfile, template): - """Import OW/DSG IJXYZ ascii format, with a Cube or RegularSurface as template.""" - cfhandle = mfile.get_cfhandle() - - if isinstance(template, (xtgeo.cube.Cube, xtgeo.surface.RegularSurface)): - logger.info("OK template") - else: - raise ValueError(f"Template is of wrong type: {type(template)}") - - nxy = template.ncol * template.nrow - ier, val = _cxtgeo.surf_import_ijxyz_tmpl( - cfhandle, template.ilines, template.xlines, nxy, 0 - ) - - if ier == -1: - raise ValueError( - f"The file {mfile.name} and template map or cube has inconsistent " - "inline and/or xlines numbering. Try importing without template " - "and use e.g. resampling instead." - ) - - if ier != 0: - raise RuntimeError("Unknown error when trying to import the IJXYZ based file!") - - val = ma.masked_greater(val, xtgeo.UNDEF_LIMIT) - - args = {} - args["xori"] = template.xori - args["xinc"] = template.xinc - args["yori"] = template.yori - args["yinc"] = template.yinc - args["ncol"] = template.ncol - args["nrow"] = template.nrow - args["rotation"] = template.rotation - args["yflip"] = template.yflip - args["values"] = val.reshape((args["ncol"], args["nrow"])) - - args["ilines"] = template._ilines.copy() - args["xlines"] = template._xlines.copy() - - mfile.cfclose() - return args + ijxyz_data = parse_ijxyz(mfile, template) + + return { + "ncol": ijxyz_data.ncol, + "nrow": ijxyz_data.nrow, + "xori": ijxyz_data.xori, + "yori": ijxyz_data.yori, + "xinc": ijxyz_data.xinc, + "yinc": ijxyz_data.yinc, + "values": ijxyz_data.values, + "ilines": ijxyz_data.ilines, + "xlines": ijxyz_data.xlines, + "yflip": ijxyz_data.yflip, + "rotation": ijxyz_data.rotation, + } def import_petromod(mfile, **_): diff --git a/src/xtgeo/surface/regular_surface.py b/src/xtgeo/surface/regular_surface.py index 8ee089a3c..47f335a88 100644 --- a/src/xtgeo/surface/regular_surface.py +++ b/src/xtgeo/surface/regular_surface.py @@ -453,6 +453,14 @@ def _read_zmap_ascii(cls, mfile, values): args = _data_reader_factory("zmap_ascii")(mfile, values=values) return cls(**args) + @classmethod + def _read_ijxyz( + cls, mfile: xtgeo._XTGeoFile, template: xtgeo.RegularSurface | xtgeo.Cube | None + ): + mfile = xtgeosys._XTGeoFile(mfile) + args = _data_reader_factory("ijxyz")(mfile, template=template) + return cls(**args) + def __repr__(self): """Magic method __repr__.""" return ( @@ -938,7 +946,7 @@ def from_file( engine: Default is "cxtgeo" which use a C backend. Optionally a pure python "python" reader will be used, which in general is slower but may be safer when reading memory streams and/or threading. - Engine is relevant for Irap binary, Irap ascii and zmap. + Keyword engine is only relevant for Irap binary, Irap ascii and zmap. Returns: Object instance. diff --git a/tests/test_surface/test_file_io.py b/tests/test_surface/test_file_io.py index 6c0ba8044..16b6c8607 100644 --- a/tests/test_surface/test_file_io.py +++ b/tests/test_surface/test_file_io.py @@ -62,8 +62,8 @@ def surfaces(draw): "irap_binary", "irap_ascii", "zmap_ascii", - # "ijxyz", # This fails horribly "petromod", + "ijxyz", "zmap", "zmap_ascii", "xtgregsurf", @@ -84,6 +84,7 @@ def test_simple_io(input_val, expected_result, fformat, engine): "irap_ascii", "irap_binary", "zmap_ascii", + "ijxyz", ]: pytest.skip("Only one engine available") surf = RegularSurface(ncol=2, nrow=2, xinc=2.0, yinc=2.0, values=input_val) @@ -105,7 +106,6 @@ def test_simple_io(input_val, expected_result, fformat, engine): "irap_ascii", "zmap", "zmap_ascii", - # "ijxyz", # This fails horribly "petromod", "xtgregsurf", ], diff --git a/tests/test_surface/test_ijxyz_spec.py b/tests/test_surface/test_ijxyz_spec.py new file mode 100644 index 000000000..cd4972de0 --- /dev/null +++ b/tests/test_surface/test_ijxyz_spec.py @@ -0,0 +1,194 @@ +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pytest +import xtgeo +from xtgeo.common import XTGeoDialog +from xtgeo.surface._regsurf_ijxyz_parser import SurfaceIJXYZ + +xtg = XTGeoDialog() +logger = xtg.basiclogger(__name__) + +if not xtg.testsetup(): + raise SystemExit + +# ============================================================================= +# Do tests +# ============================================================================= + +TPATH = xtg.testpathobj + +SURF1 = TPATH / "surfaces/etc/h0.dat" + +SIMPLE_VALUES1 = [[2, 3, 4, 5], [6, 5, 4, 99], [44.0, 33.0, 22.0, 11.0]] +SIMPLE_MASK1 = [ + [True, False, False, False], + [True, False, True, False], + [False, False, True, False], +] + +SIMPLE_VALUES2 = range(300) +SIMPLE_MASK2 = [0, 1, 4, 5, 19, 95, 96, 97, 121] + +# real file(s) +IJXYZFILE1 = TPATH / "surfaces/etc/ijxyz1.dat" # OW IJXYZ format with comments +OTHERFILE1 = TPATH / "surfaces/reek/1/topreek_rota.gri" # for template test + + +@pytest.fixture(name="simple1") +def fixture_simple1(tmp_path: Path): + """Test reading small simple case.""" + surf = xtgeo.RegularSurface( + 3, + 4, + 55.0, + 133.0, + xori=12345.0, + yori=998877.0, + values=np.ma.MaskedArray(SIMPLE_VALUES1, mask=SIMPLE_MASK1), + rotation=33, + yflip=1, + ) + dfr = surf.get_dataframe(ij=True) + target = tmp_path / "ijxyz.dat" + dfr.to_csv(target, sep=" ", index=False, header=False) + + return (surf, np.loadtxt(target, comments=["@", "#", "EOB"])) + + +@pytest.fixture(name="simple2") +def fixture_simple2(tmp_path: Path): + """Test somewhat larger case""" + + values = np.array(SIMPLE_VALUES2) + values = np.ma.masked_where(np.isin(np.arange(values.size), SIMPLE_MASK2), values) + + surf = xtgeo.RegularSurface( + 30, + 10, + 50.0, + 100.0, + xori=10000.0, + yori=90000.0, + values=values, + rotation=0, + yflip=1, + ilines=np.array(range(233, 233 + 90, 3)), + xlines=np.array(range(44, 44 + 20, 2)), + ) + dfr = surf.get_dataframe(ij=True) + ix = dfr.IX - 1 + jy = dfr.JY - 1 + ilines = surf.ilines[ix] + xlines = surf.xlines[jy] + dfr["IL"] = ilines + dfr["XL"] = xlines + dfr = dfr[["IL", "XL", "X_UTME", "Y_UTMN", "VALUES"]] + target = tmp_path / "ijxyz2.dat" + dfr.to_csv(target, sep=" ", index=False, header=False) + + return (surf, np.loadtxt(target, comments=["@", "#", "EOB"])) + + +def test_simple1_ijxyz_read(simple1: tuple): + surf, data = simple1 # truth surf, and data arrays read from file + inline = data[:, 0].astype("int32") + xline = data[:, 1].astype("int32") + x_arr = data[:, 2].astype("float64") + y_arr = data[:, 3].astype("float64") + z_arr = data[:, 4].astype("float64") + + ijxyz = SurfaceIJXYZ(x_arr, y_arr, z_arr, inline, xline) + + assert ijxyz.yflip == surf.yflip + assert ijxyz.xori == pytest.approx(surf.xori) + assert ijxyz.ncol == surf.ncol + assert ijxyz.nrow == surf.nrow + assert ijxyz.xinc == pytest.approx(surf.xinc) + assert ijxyz.yinc == pytest.approx(surf.yinc) + np.testing.assert_array_almost_equal(ijxyz.values, surf.values) + np.testing.assert_array_almost_equal(ijxyz.values.mask, surf.values.mask) + + +def test_simple2_ijxyz_read(simple2: tuple): + surf, data = simple2 # truth surf, and data arrays read from file + inline = data[:, 0].astype("int32") + xline = data[:, 1].astype("int32") + x_arr = data[:, 2].astype("float64") + y_arr = data[:, 3].astype("float64") + z_arr = data[:, 4].astype("float64") + + ijxyz = SurfaceIJXYZ(x_arr, y_arr, z_arr, inline, xline) + + assert ijxyz.yflip == surf.yflip + assert ijxyz.xori == pytest.approx(surf.xori) + assert ijxyz.ncol == surf.ncol + assert ijxyz.nrow == surf.nrow + assert ijxyz.xinc == pytest.approx(surf.xinc) + assert ijxyz.yinc == pytest.approx(surf.yinc) + + np.testing.assert_array_almost_equal(ijxyz.values, surf.values) + np.testing.assert_array_almost_equal(ijxyz.values.mask, surf.values.mask) + + np.testing.assert_array_almost_equal(ijxyz.ilines, surf.ilines) + np.testing.assert_array_almost_equal(ijxyz.xlines, surf.xlines) + + +def test_ijxyzfile1_read(): + """Test a real file exported from OW""" + data = np.loadtxt(IJXYZFILE1, comments=["@", "#", "EOB"]) + inline = data[:, 0].astype("int32") + xline = data[:, 1].astype("int32") + x_arr = data[:, 2].astype("float64") + y_arr = data[:, 3].astype("float64") + z_arr = data[:, 4].astype("float64") + + ijxyz = SurfaceIJXYZ(x_arr, y_arr, z_arr, inline, xline) + + assert ijxyz.values.mean() == pytest.approx(5037.5840, abs=0.001) + assert ijxyz.ncol == 51 + assert ijxyz.yflip == -1 + + +def test_ijxyz_io_file_with_template(tmp_path): + """Test i/o with and without template""" + + surf0 = xtgeo.surface_from_file(OTHERFILE1) + + surf0.ilines += 22 + surf0.ilines *= 2 + + surf0.xlines += 33 + surf0.xlines *= 3 + + file1 = tmp_path / "other_x1.dat" + surf0.to_file(file1, fformat="ijxyz") + + surf1 = xtgeo.surface_from_file(file1, fformat="ijxyz") + + # assertions here are intentional; with undef cells, ijxy cannot reproduce + # ncol, nrow and xori, yori for original file due to lack of information + assert surf1.ncol == surf0.ncol - 29 + assert surf1.nrow == surf0.nrow - 9 + assert surf1.rotation == pytest.approx(surf0.rotation) + assert surf1.xinc == pytest.approx(surf0.xinc) + assert surf1.yinc == pytest.approx(surf0.yinc) + assert surf1.xori != pytest.approx(surf0.xori) + assert surf1.yori != pytest.approx(surf0.yori) + + # try ijxyz with template instead; here assertions are equal + surf2 = xtgeo.surface_from_file(file1, fformat="ijxyz", template=surf0) + assert surf2.ncol == surf0.ncol + assert surf2.nrow == surf0.nrow + assert surf2.rotation == pytest.approx(surf0.rotation) + assert surf2.xinc == pytest.approx(surf0.xinc) + assert surf2.yinc == pytest.approx(surf0.yinc) + assert surf2.xori == pytest.approx(surf0.xori) + assert surf2.yori == pytest.approx(surf0.yori) + np.testing.assert_array_equal(surf0.ilines, surf2.ilines) + np.testing.assert_array_equal(surf0.xlines, surf2.xlines) + + np.testing.assert_array_almost_equal(surf0.values, surf2.values) + assert surf1.values.mean() == pytest.approx(surf2.values.mean()) diff --git a/tests/test_surface/test_regular_surface_bytesio.py b/tests/test_surface/test_regular_surface_bytesio.py index 4dcfd1d11..9918a8e62 100644 --- a/tests/test_surface/test_regular_surface_bytesio.py +++ b/tests/test_surface/test_regular_surface_bytesio.py @@ -215,3 +215,30 @@ def test_export_import_hdf5_bytesio(tmp_path): xsurf4 = xtgeo.surface_from_file(stream, fformat="hdf") assert xsurf4.values.mean() == xsurf.values.mean() + + +def test_export_import_ijxyz_bytesio(tmp_path): + """Test ijxyz format via files and memory streams.""" + # just the input, and save as hdf5 + xsurf = xtgeo.surface_from_file(TESTSET2, fformat="irap_ascii") + assert xsurf.ncol == 554 + assert xsurf.nrow == 451 + + ijxyzfile = tmp_path / "ijxyz_tset2.fgr" + xsurf.to_file(ijxyzfile, fformat="ijxyz") + + xsurf2 = xtgeo.surface_from_file(ijxyzfile, fformat="ijxyz") + + # if surface is undefined in areas at border, the ijxyz format will not be + # able to know the original ncol, nrow, so a reduction may be expected + assert xsurf2.ncol == 525 + assert xsurf2.nrow == 442 + assert xsurf2.values.mean() == pytest.approx(xsurf.values.mean()) + + stream = io.BytesIO() + + xsurf.to_file(stream, fformat="ijxyz") + + xsurf3 = xtgeo.surface_from_file(stream, fformat="ijxyz") + assert xsurf3.ncol == xsurf2.ncol + assert xsurf3.values.mean() == pytest.approx(xsurf2.values.mean())