Questions about reading USG generated concentration file that has multiple components

[DataGeoranger]

Hello, my USG concentration files have two components (1 - Sulfate, 2 - Selenium). When I used flopy binary module to read the file, it only created the object for the second concentration. How do I get the object created for the first component?

import flopy.utils.binaryfile as bf
conc_file = 'test.con'
concobj = bf.HeadUFile(os.path.join(project_folder, conc_file), text = 'conc', verbose=True)

[cnicol-gwlogic]

I don't think this is possible with flopy in its current form. I suggest you use mfusg's "multifile" option to save one con file for each species.

[langevin-usgs]

Any chance this works?

conc_file = 'CONC CONC01.CON'
concobj = bf.HeadUFile(os.path.join(project_folder, conc_file), text='conc01', precision='single', verbose=True)

[cnicol-gwlogic]

Thanks @langevin-usgs - gave it a shot but another issue with that one:
Traceback (most recent call last): File "C:\Users\Chris\Downloads\test_con\test.py", line 5, in <module> concobj = bf.HeadUFile( File "C:\Python39\lib\site-packages\flopy\utils\binaryfile.py", line 857, in __init__ super().__init__(filename, precision, verbose, kwargs) File "C:\Python39\lib\site-packages\flopy\utils\binaryfile.py", line 439, in __init__ super().__init__(filename, precision, verbose, kwargs) File "C:\Python39\lib\site-packages\flopy\utils\datafile.py", line 208, in __init__ self._build_index() File "C:\Python39\lib\site-packages\flopy\utils\binaryfile.py", line 487, in _build_index self.file.seek(databytes, 1) OSError: [Errno 22] Invalid argument
I tried 3.3.7 and latest dev branch. I'm guessing there's something odd going on on the usg-transport side.

[langevin-usgs]

Sorry for the distraction, and thanks for giving that a shot. If the binary headfile reader is dying on the line if t.upper() not in asciiset: then I'm guessing the string that is written to text in the binary file has some uninterpretable characters, which should probably be spaces. That can happen on the Fortran side if the text string is not constructed properly. Seems like your guess about it being on the usg-transport side makes sense. We could probably hack flopy to handle that, but it would be cleaner to fix on the usg side. Thanks for looking into this.

[cnicol-gwlogic]

Beauty thanks for that pointer, I'll dig into it.

[DataGeoranger]

Thank you both @cnicol-gwlogic @langevin-usgs. I emailed Sorab about multifile option. will follow up

[ougx]

looks like the files are not flat binary files. If not written using access='stream' or ifort's binary format, it inserts record length before writing the actual data. maybe rerun the model with standard executables and try again.

[DataGeoranger]

Thank you, Shixiong!

[hjia1005]

@cnicol-gwlogic Chris, you helped me with this last year by referring me to a link posted maybe two years ago.
@DataGeoranger Below is the function to extract concentrations for multi-species from the binary file. It was modified from the code provided in that link but I can't find the source now.

def get_conc(ucnobj, text, nlay, sp):
    totim=dict_totim[sp]
    indices = np.where((ucnobj.recordarray["totim"] == totim) & (ucnobj.recordarray["text"] == text))[0]
    data = nlay * [None]
    data_total=[]
    for idx in indices:
        ipos = ucnobj.iposarray[idx]
        ilay = ucnobj.recordarray["ilay"][idx]
        nstrt = ucnobj.recordarray["ncol"][idx]
        nend = ucnobj.recordarray["nrow"][idx]
        npl = nend - nstrt + 1
        ucnobj.file.seek(ipos, 0)
        data[ilay-1] =bf.binaryread(ucnobj.file, ucnobj.realtype, shape=(npl,))
        data_total+=list(data[ilay-1])
    arr_data_total=np.array(data_total)
    return arr_data_total

My model has 5 layers and 65 stress periods. To meet my own requirements, I created a dictionary of sp and totim. I am sure it can be done in a smarter way.

lst_sp=list(np.arange(0,65,1))   
lst_totim=[31., 59., 90., 120., 151., 181., 212., 243., 273., 304., 334.,
           365., 396., 424., 455., 485., 516., 546., 577., 608., 638., 669., 699., 730.,
           761., 790., 821., 851., 882., 912., 943., 974., 1004., 1035., 1065., 1096.,
           1127., 1155., 1186., 1216.,1247., 1277., 1308., 1339., 1369., 1400., 1430.,
           1461., 1492., 1520., 1551., 1581., 1612., 1642., 1673., 1704., 1734., 1765.,
           1795., 1826., 2191., 4017., 7670., 25932., 62456.]
dict_totim=dict(zip(lst_sp,lst_totim))

ucnobj = bf.HeadUFile('path_to_your_concentration_file', precision='single', text='conc')
text_component1 = ucnobj.recordarray['text'][0]    # in my case, the text_component1 is "CONC01"
text_component2 = ucnobj.recordarray['text'][5]    # in my case, the text_component2 is "CONC02"

# Then you can extract concentration data for specified species and sp.
conc01_sp1 = get_conc(ucnobj, text=text_component1, nlay=5, sp=0)
conc01_sp65 = get_conc(ucnobj, text=text_component1, nlay=5, sp=64)
conc02_sp1 = get_conc(ucnobj, text=text_component2, nlay=5, sp=0)
conc02_sp65 = get_conc(ucnobj, text=text_component2, nlay=5,sp=64)

If you have more than two species, I guess you can use text_component3=ucnobj.recordarray['text'][10] for the third species.
The increment of 5 should be because of 5 layers in the model.

Worth mentioning that I tried to use ucnobj.get_times( ) and ucnobj.get_kstpkper( ) but both results looked suspicious, while results from hdobj.get_times( ) and hdobj.get_kstpkper( ) are correct.

I just tested the code and it still works. Hope it helps.

[DataGeoranger]

@hjia1005 Thanks!

[sorab]

Sorry for this confusion. This was a USG-Transport issue. We had changed to using UNFORMATTED form and STREAM access as an approach that was "compatible" with both Ifort and GNU compilers. Unfortunately, that created some issues with some of the files depending on how the other options for opening files were handled. I have changed it back to BINARY form and SEQUENTIAL access and the code should be on the GSI website in a day or two.

[DataGeoranger]

Thank you, Sorab! Look forward to the new version.