search

modflowpy / flopy

A Python package to create, run, and post-process MODFLOW-based models.
https://flopy.readthedocs.io
Other
517 stars 313 forks source link

Read timeseries in a FormattedHeadFile gives wrong results #116

Closed rabbl closed 8 years ago

rabbl commented 8 years ago

I have a 1,8GB fhd-results-file which I can read without problems with the get_data method. Reading the heads of a specific cell gives me the correct results. Reading the timeseries of the same cell gives me other values.

I'll post the test-script, the output of the script and the header and first row of the fhd-file.

Script:

import os
import flopy.modflow as mf
import flopy.utils as fu

name = 'Q_63'
layer = 4
row = 81
column = 87

workspace = os.path.join('data')
modelname = 'Transient_361703'

ml = mf.Modflow.load(modelname+'.nam', exe_name='mfnwt', model_ws=workspace, verbose=False)
fhd = fu.formattedfile.FormattedHeadFile(os.path.join(workspace, modelname+'.fhd'), precision='double', verbose=False)

# read the values with .get_data
print()
print("Read the values with .get_data")
print("T: time, L: Layer, R: Row, Head")
times = fhd.get_times()
for x in range(0, 19):
    heads = fhd.get_data(totim=times[x], mflay=layer-1)
    print("T:{}, L:{}, R:{}, C:{}, Head:{}".format(times[x], layer, row, column, heads[row-1][column-1]))

# read the timeserie
print()
print("Read the values with .get_ts")
print("T: time, L: Layer, R: Row, Head")
heads = fhd.get_ts((layer-1, row-1, column-1))
for x in range(0, 20):
    print("T:{}, L:{}, R:{}, C:{}, Head:{}".format(heads[x][0], layer, row, column, heads[x][1]))

Script-Output:

Read the values with .get_data
T: time, L: Layer, R: Row, Head
T:1.0, L:4, R:81, C:87, Head:-12.0397
T:2.0, L:4, R:81, C:87, Head:-12.2902
T:3.0, L:4, R:81, C:87, Head:-12.4616
T:4.0, L:4, R:81, C:87, Head:-12.5935
T:5.0, L:4, R:81, C:87, Head:-12.7015
T:6.0, L:4, R:81, C:87, Head:-12.7931
T:7.0, L:4, R:81, C:87, Head:-12.8725
T:8.0, L:4, R:81, C:87, Head:-12.9427
T:9.0, L:4, R:81, C:87, Head:-13.0054
T:10.0, L:4, R:81, C:87, Head:-13.0616
T:11.0, L:4, R:81, C:87, Head:-13.1115
T:12.0, L:4, R:81, C:87, Head:-13.1574
T:13.0, L:4, R:81, C:87, Head:-13.2001
T:14.0, L:4, R:81, C:87, Head:-13.2399
T:15.0, L:4, R:81, C:87, Head:-13.2771
T:16.0, L:4, R:81, C:87, Head:-13.3119
T:17.0, L:4, R:81, C:87, Head:-13.3446
T:18.0, L:4, R:81, C:87, Head:-13.3753
T:19.0, L:4, R:81, C:87, Head:-13.4042

Read the values with .get_ts
T: time, L: Layer, R: Row, Head
T:1.0, L:4, R:81, C:87, Head:-2e+20
T:2.0, L:4, R:81, C:87, Head:-2e+20
T:3.0, L:4, R:81, C:87, Head:-2e+20
T:4.0, L:4, R:81, C:87, Head:-2e+20
T:5.0, L:4, R:81, C:87, Head:-2e+20
T:6.0, L:4, R:81, C:87, Head:-2e+20
T:7.0, L:4, R:81, C:87, Head:-2e+20
T:8.0, L:4, R:81, C:87, Head:-2e+20
T:9.0, L:4, R:81, C:87, Head:-2e+20
T:10.0, L:4, R:81, C:87, Head:-2e+20
T:11.0, L:4, R:81, C:87, Head:-2e+20
T:12.0, L:4, R:81, C:87, Head:-2e+20
T:13.0, L:4, R:81, C:87, Head:-2e+20
T:14.0, L:4, R:81, C:87, Head:-2e+20
T:15.0, L:4, R:81, C:87, Head:-2e+20
T:16.0, L:4, R:81, C:87, Head:-2e+20
T:17.0, L:4, R:81, C:87, Head:-2e+20
T:18.0, L:4, R:81, C:87, Head:-2e+20
T:19.0, L:4, R:81, C:87, Head:-2e+20
T:20.0, L:4, R:81, C:87, Head:-2e+20

Process finished with exit code 0

FHD-File 

    1    1   1.000000E+00   1.000000E+00             HEAD   165   175     1 (10(1X1PE13.5))
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20  -2.00000E+20
...
...
...
langevin-usgs commented 8 years ago

Thanks, Ralf. I've asked the person who wrote the formatted file head reader to take a look.

langevin-usgs commented 8 years ago

Ralph, I believe your issue has been fixed on the develop branch. Please verify. Also, here is note from colleague about your problem:

This should solve Ralf's problem with get_ts. However, it looks like Ralf is using get_data incorrectly. Instead of:

heads = fhd.get_data(idx=totime, mflay=layer-1)

he might mean:

heads = fhd.get_data(totim=totime, mflay=layer-1)

Since layers have headers, each layer is indexed separately. Therefore passing in an index and a layer number is redundant.

rabbl commented 8 years ago

It's working. Great. Here the output of the same script:

Read the values with .get_data
T: time, L: Layer, R: Row, Head
T:1.0, L:4, R:81, C:87, Head:-12.0397
T:2.0, L:4, R:81, C:87, Head:-12.2902
T:3.0, L:4, R:81, C:87, Head:-12.4616
T:4.0, L:4, R:81, C:87, Head:-12.5935
T:5.0, L:4, R:81, C:87, Head:-12.7015
T:6.0, L:4, R:81, C:87, Head:-12.7931
T:7.0, L:4, R:81, C:87, Head:-12.8725
T:8.0, L:4, R:81, C:87, Head:-12.9427
T:9.0, L:4, R:81, C:87, Head:-13.0054
T:10.0, L:4, R:81, C:87, Head:-13.0616
T:11.0, L:4, R:81, C:87, Head:-13.1115
T:12.0, L:4, R:81, C:87, Head:-13.1574
T:13.0, L:4, R:81, C:87, Head:-13.2001
T:14.0, L:4, R:81, C:87, Head:-13.2399
T:15.0, L:4, R:81, C:87, Head:-13.2771
T:16.0, L:4, R:81, C:87, Head:-13.3119
T:17.0, L:4, R:81, C:87, Head:-13.3446
T:18.0, L:4, R:81, C:87, Head:-13.3753
T:19.0, L:4, R:81, C:87, Head:-13.4042

Read the values with .get_ts
T: time, L: Layer, R: Row, Head
T:1.0, L:4, R:81, C:87, Head:-12.0397
T:2.0, L:4, R:81, C:87, Head:-12.2902
T:3.0, L:4, R:81, C:87, Head:-12.4616
T:4.0, L:4, R:81, C:87, Head:-12.5935
T:5.0, L:4, R:81, C:87, Head:-12.7015
T:6.0, L:4, R:81, C:87, Head:-12.7931
T:7.0, L:4, R:81, C:87, Head:-12.8725
T:8.0, L:4, R:81, C:87, Head:-12.9427
T:9.0, L:4, R:81, C:87, Head:-13.0054
T:10.0, L:4, R:81, C:87, Head:-13.0616
T:11.0, L:4, R:81, C:87, Head:-13.1115
T:12.0, L:4, R:81, C:87, Head:-13.1574
T:13.0, L:4, R:81, C:87, Head:-13.2001
T:14.0, L:4, R:81, C:87, Head:-13.2399
T:15.0, L:4, R:81, C:87, Head:-13.2771
T:16.0, L:4, R:81, C:87, Head:-13.3119
T:17.0, L:4, R:81, C:87, Head:-13.3446
T:18.0, L:4, R:81, C:87, Head:-13.3753
T:19.0, L:4, R:81, C:87, Head:-13.4042

Process finished with exit code 0

Thanks for the advice to use totim. I had updated the script-file yet. When will you merge it to the master?

langevin-usgs commented 8 years ago

We will try to make a new flopy release in the next month or so at which time we will merge develop into master. Thanks for the bug report.

rabbl commented 8 years ago

In the actual master i have the same problem as reported again. Is this correct?

Read the values with .get_data
T: time, L: Layer, R: Row, Head
T:1.0, L:4, R:81, C:87, Head:-12.0397
T:2.0, L:4, R:81, C:87, Head:-12.2902
T:3.0, L:4, R:81, C:87, Head:-12.4616
T:4.0, L:4, R:81, C:87, Head:-12.5935
T:5.0, L:4, R:81, C:87, Head:-12.7015
T:6.0, L:4, R:81, C:87, Head:-12.7931
T:7.0, L:4, R:81, C:87, Head:-12.8725
T:8.0, L:4, R:81, C:87, Head:-12.9427
T:9.0, L:4, R:81, C:87, Head:-13.0054
T:10.0, L:4, R:81, C:87, Head:-13.0616
T:11.0, L:4, R:81, C:87, Head:-13.1115
T:12.0, L:4, R:81, C:87, Head:-13.1574
T:13.0, L:4, R:81, C:87, Head:-13.2001
T:14.0, L:4, R:81, C:87, Head:-13.2399
T:15.0, L:4, R:81, C:87, Head:-13.2771
T:16.0, L:4, R:81, C:87, Head:-13.3119
T:17.0, L:4, R:81, C:87, Head:-13.3446
T:18.0, L:4, R:81, C:87, Head:-13.3753
T:19.0, L:4, R:81, C:87, Head:-13.4042
()
Read the values with .get_ts
T: time, L: Layer, R: Row, Head
T:1.0, L:4, R:81, C:87, Head:-2e+20
T:2.0, L:4, R:81, C:87, Head:-2e+20
T:3.0, L:4, R:81, C:87, Head:-2e+20
T:4.0, L:4, R:81, C:87, Head:-2e+20
T:5.0, L:4, R:81, C:87, Head:-2e+20
T:6.0, L:4, R:81, C:87, Head:-2e+20
T:7.0, L:4, R:81, C:87, Head:-2e+20
T:8.0, L:4, R:81, C:87, Head:-2e+20
T:9.0, L:4, R:81, C:87, Head:-2e+20
T:10.0, L:4, R:81, C:87, Head:-2e+20
T:11.0, L:4, R:81, C:87, Head:-2e+20
T:12.0, L:4, R:81, C:87, Head:-2e+20
T:13.0, L:4, R:81, C:87, Head:-2e+20
T:14.0, L:4, R:81, C:87, Head:-2e+20
T:15.0, L:4, R:81, C:87, Head:-2e+20
T:16.0, L:4, R:81, C:87, Head:-2e+20
T:17.0, L:4, R:81, C:87, Head:-2e+20
T:18.0, L:4, R:81, C:87, Head:-2e+20
T:19.0, L:4, R:81, C:87, Head:-2e+20
T:20.0, L:4, R:81, C:87, Head:-2e+20

Process finished with exit code 0
langevin-usgs commented 8 years ago

Did you check the develop branch? I don't think we have released this as master yet. I think we will be doing so this week.

rabbl commented 8 years ago

Thanks... ;)