Defining and plotting lake grid when doing the preprocessing of bathymetric data in the tutorial

[Elias Pentti]

I just started learning how to use PITLAKQ and am also a newcomer as a Python user. I installed miniconda3 this week and use PITLAKQ with Python version 3.9.16. I'm going through the Tutorial on the PITLAKQ webpage, and at the "Preprocess the bathymetry data" stage, started to get unexpected results. I created a simple lake contour and in reservoir.txt, defined lake grid columns as

Columns RW
1 200
2 250
3 325
4 400
5 500
6 600
7 700
8 800
9 860
10 940

Then I ran the preprocessing and show_bath programs and got a plot like this

which is not as I think it should be. There are gaps between the boxes and the grid starts from -50 m in the horizontal dimension instead of 0. From bath.nc I managed to read that the segment_lenght list is

"100, 50, 75, 75, 100, 100, 100, 100, 60, 80, 100" and after making show_bath.py print the value of "pos", got an output

[  0. 100. 150. 225. 300. 400. 500. 600. 700. 760. 840.]

So, it seems to me that first, the width of the first column in the lake grid is 100 m, not 200 m as I had defined, and second, show_bath had drawn the rectangles so that their middle points are in the horizontal position where actually the left edge should be, and third, the final column is too narrow to reach the edge of my contour grid, which is 1000 m wide just like the one in the tutorial.

Could someone explain why the program behaves like this? Is it buggy or have I done something wrong?

[Elias Pentti]

Now I found the discussion in July 2021 about the default setting issue concerning pylab.bar and added "align='edge'" in show_bath.py, which solved the plotting error. The other issue still remains, the first and last columns in the lake grid become 100 m wide regardless of what I define in reservoir.txt.

[Dr. Mike Mueller]

Hi Elias,

Thanks for finding this. Looks like matplotlib changed a default value.
It is a display problem, the bathymetry data are fine.

I fixed it and created the enclosed picture with this command:

show_bath.main("preprocessing/output/bath.nc", bottom=150,
properties={'xoffset': 100}
)

The gray cells are inactive. The first and last segments always have a width
of 100 m. The values for "Columns RW" are x coordinates of the left border of
a segment.

I need a few days to make a bug-fix release. For now, you can fix it yourself.

In your Notebook type:

show_bath

this will show the full file path:

<module 'pitlakq.postprocessing.show_bath' from
'/path/to/pitlakq/postprocessing/show_bath.py'>

Open this file and search for this code (should be close to line 84):

pylab.bar(pos, height, self.dxs, bottom=bottom, color=color,
linewidth=self.linewidth, edgecolor=self.edgecolor)

and add `align='edge'`, i.e. it should look like this:

pylab.bar(pos, height, self.dxs, bottom=bottom, color=color,
linewidth=self.linewidth, edgecolor=self.edgecolor,
align='edge')

Let me know if it works.

Best,
Mike

Am 25.05.23 um 12:08 schrieb Elias Pentti:

> Now I found the discussion in July 2021 about the default setting issue
> concerning pylab.bar and added "align='edge'" in show_bath.py, which solved the
> plotting error. The other issue still remains, the first and last columns in
> the lake grid become 100 m wide regardless of what I define in reservoir.txt.
>
> torstai 25. toukokuuta 2023 klo 12.54.55 UTC+3 Elias Pentti kirjoitti:
>
> I just started learning how to use PITLAKQ and am also a newcomer as a
> Python user. I installed miniconda3 this week and use PITLAKQ with Python
> version 3.9.16. I'm going through the Tutorial on the PITLAKQ webpage, and
> at the "Preprocess the bathymetry data" stage, started to get unexpected
> results. I created a simple lake contour and in reservoir.txt, defined lake
> grid columns as
> Columns RW
> 1200
> 2250
> 3325
> 4400
> 5500
> 6600
> 7700
> 8800
> 9860
> 10940
> Then I ran the preprocessing and show_bath programs and got a plot like this

> lakegrid.png
> which is not as I think it should be. There are gaps between the boxes and
> the grid starts from -50 m in the horizontal dimension instead of 0. From

> bath.nc <http://bath.nc> I managed to read that the segment_lenght list is

> "100, 50, 75, 75, 100, 100, 100, 100, 60, 80, 100" and after making
> show_bath.py print the value of "pos", got an output
> [  0. 100. 150. 225. 300. 400. 500. 600. 700. 760. 840.]
> So, it seems to me that first, the width of the first column in the lake
> grid is 100 m, not 200 m as I had defined, and second, show_bath had drawn
> the rectangles so that their middle points are in the horizontal position
> where actually the left edge should be, and third, the final column is too
> narrow to reach the edge of my contour grid, which is 1000 m wide just like
> the one in the tutorial.
> Could someone explain why the program behaves like this? Is it buggy or
> have I done something wrong?
>

> --
> You received this message because you are subscribed to the Google Groups
> "pitlakq-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email
> to pitlakq-user...@googlegroups.com
> <mailto:pitlakq-user...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/pitlakq-users/b4606b52-2021-4d5d-8567-ec9377c7dfaen%40googlegroups.com <https://groups.google.com/d/msgid/pitlakq-users/b4606b52-2021-4d5d-8567-ec9377c7dfaen%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Elias Pentti]

Hello Mike!

The correction to show_bath worked fine, thank you!

Thanks also for the information that the first and last segments of the lake grid must always be 100 m wide. Maybe you could tell that in the tutorial also.

Best regards,

   Elias

[Dr. Mike Mueller]

Hi Elias,

Am 26.05.23 um 10:02 schrieb Elias Pentti:

> Hello Mike!
>
> The correction to show_bath worked fine, thank you!

Nice.

> Thanks also for the
> information that the first and last segments of the lake grid must always be
> 100 m wide. Maybe you could tell that in the tutorial also.

That is not totally correct. The first and last segments are always inactive.
So the width does not matter at all. If you multiply a number by zero the
result is always zero no matter what number it is. ;) The width of 100 m is set
automatically. The user has no way to influence it. It is just for display.
This works for typical segments with a width of 50 to 1000 m. It will look
strange for if the segments are very small like 1 or 2 m. Maybe I can improve
the algorithm and set the inactive segments to the same width as the smallest
active cell. There is always room for improvement. ;)

Best,
Mike

>
> Best regards, Elias
>
> torstai 25. toukokuuta 2023 klo 18.52.48 UTC+3 mmueller kirjoitti:
>
> Hi Elias,
>
> Thanks for finding this. Looks like matplotlib changed a default value. It
> is a display problem, the bathymetry data are fine.
>
> I fixed it and created the enclosed picture with this command:
>

> show_bath.main("preprocessing/output/bath.nc <http://bath.nc>", bottom=150,

>> instead of 0. From bath.nc <http://bath.nc> <http://bath.nc

>> <http://bath.nc>> I managed to
> read that the segment_lenght list is
>> "100, 50, 75, 75, 100, 100, 100, 100, 60, 80, 100" and after making
>> show_bath.py print the value of "pos", got an output [ 0. 100. 150. 225.
>> 300. 400. 500. 600. 700. 760. 840.] So, it seems to me that first, the
>> width of the first column in the lake grid is 100 m, not 200 m as I had
>> defined, and second, show_bath had drawn the rectangles so that their
>> middle points are in the horizontal position where actually the left edge
>> should be, and third, the final column is too narrow to reach the edge of
>> my contour grid, which is 1000 m wide just like the one in the tutorial.
>> Could someone explain why the program behaves like this? Is it buggy or
>> have I done something wrong?
>>
>> -- You received this message because you are subscribed to the Google
>> Groups "pitlakq-users" group. To unsubscribe from this group and stop
>> receiving emails from it, send
> an email
>> to pitlakq-user...@googlegroups.com
>> <mailto:pitlakq-user...@googlegroups.com>. To view this discussion on the
>> web visit
>>
> https://groups.google.com/d/msgid/pitlakq-users/b4606b52-2021-4d5d-8567-ec9377c7dfaen%40googlegroups.com
> <https://groups.google.com/d/msgid/pitlakq-users/b4606b52-2021-4d5d-8567-ec9377c7dfaen%40googlegroups.com>

> <https://groups.google.com/d/msgid/pitlakq-users/b4606b52-2021-4d5d-8567-ec9377c7dfaen%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pitlakq-users/b4606b52-2021-4d5d-8567-ec9377c7dfaen%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
> -- You received this message because you are subscribed to the Google
> Groups "pitlakq-users" group. To unsubscribe from this group and stop
> receiving emails from it, send an email to
> pitlakq-user...@googlegroups.com
> <mailto:pitlakq-user...@googlegroups.com>. To view this discussion
> on the web visit

> https://groups.google.com/d/msgid/pitlakq-users/09af4bcf-3002-453c-8f23-413135d4e584n%40googlegroups.com
> <https://groups.google.com/d/msgid/pitlakq-users/09af4bcf-3002-453c-8f23-413135d4e584n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Elias Pentti]

Hello!

Now I got to the first simple model runs. First I had to do the fix to numpy.zeros in the file ncoutput.py that you gave last month to another user, because I encountered the exact same error. Now, the model runs, but when I follow the tutorial to visualize the resulting development of temperature, I get unexpected results again. Below is the top part of the plot I get. The upper edge of the white gap in the middle corresponds to the calculated lake level, which gradually rises from 195 m to over 196 m because I'm only considering precipitation. In reservoir.txt, my layer boundaries are at 140, 145, ... , 195, 200 m with even 5 m intervals, but when I read the "zu" variable from out.nc (supposed to be the bottom elevations of the layers), it's {201, 196, 191, 186, 181, 176, 171, 166, 161, 156, 151, 146, 141, 140}.

After this curious result, I experimented a little and made the layer structure uneven, setting the layers at 140, 147, 153, 158, 166, 178, 189, 193, 197.5 and 200 m. This time "zu" in out.nc is {201, 194, 188, 183, 175, 163, 152, 148, 143.5, 141, 140}. If you look at the intervals between these numbers, they are 7, 6, 5, 8, 12, 11, 4, 4.5, 2.5 and 1, the same sequence as in the levels I had set in reservoir.txt, but in the opposite order, as if the program had somehow flipped the layer heights from up to down at some point. Is this another consequence of a version incompatibility or something else?

Elias

[Dr. Mike Mueller]

Hi Elias,

Looks like you changed the number of layers.
Did adjust the values in input/w2/w2.yaml accordingly?
These are the ones from the tutorial model:

general:
titel:
value: A simple PITLAKQ test model
bounds:
number_of_segments:
value: 10
number_of_layers:
value: 52
number_of_constituents_segments:
value: 10
number_of_constituents_layers:
value: 52
times:
start:
value: 01.01.1998
end:
value: 31.01.2009
---
!File
name: bath.nc
format: netcdf
group_name: bathymetry
---
bathymetry:
starting_water_level:
value: 195.0
---
branch_geometry:
branch_name:
value:
- test lake # no w2 variable
branch_upstream_segments:
value:
- 2
branch_downstream_segments:
value:
- 8

Many values relate to the ones described in detail in the CE-QUL-W2 manual:
https://gitlab.com/hydrocomputing/pitlakq/-/blob/master/doc/w2v2-manual.pdf

Best,
Mike

Am 30.05.23 um 15:13 schrieb Elias Pentti:

> Hello!
>
> Now I got to the first simple model runs. First I had to do the fix to
> numpy.zeros in the file ncoutput.py that you gave last month to another
> user, because I encountered the exact same error. Now, the model runs, but
> when I follow the tutorial to visualize the resulting development of
> temperature, I get unexpected results again. Below is the top part of the
> plot I get. The upper edge of the white gap in the middle corresponds to the
> calculated lake level, which gradually rises from 195 m to over 196 m
> because I'm only considering precipitation. In reservoir.txt, my layer
> boundaries are at 140, 145, ... , 195, 200 m with even 5 m intervals, but
> when I read the "zu" variable from out.nc (supposed to be the bottom
> elevations of the layers), it's {201, 196, 191, 186, 181, 176, 171, 166,

> 161, 156, 151, 146, 141, 140}. T_result_top.png After this curious result, I

> <http://bath.nc <http://bath.nc>>", bottom=150,

> <http://bath.nc>> <http://bath.nc <http://bath.nc>

> <https://groups.google.com/d/msgid/pitlakq-users/09af4bcf-3002-453c-8f23-413135d4e584n%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pitlakq-users/09af4bcf-3002-453c-8f23-413135d4e584n%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
> -- You received this message because you are subscribed to the Google
> Groups "pitlakq-users" group. To unsubscribe from this group and stop
> receiving emails from it, send an email to
> pitlakq-user...@googlegroups.com
> <mailto:pitlakq-user...@googlegroups.com>. To view this discussion
> on the web visit

> https://groups.google.com/d/msgid/pitlakq-users/ffba33cf-cc75-4364-af44-8916f5baca65n%40googlegroups.com
> <https://groups.google.com/d/msgid/pitlakq-users/ffba33cf-cc75-4364-af44-8916f5baca65n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Elias Pentti]

Yes, I did update the number_of_layers and number_of_constituent_layers values in w2.yaml so that they correspond to my bathymetry grid, as well as bottom_elevation in waterbody_coordinates to equal the lowest layer boundary defined in reservoir.txt just like in the tutorial.

To be clear, my reservoir.txt reads

Layers ZU
1 140
2 147
3 153
4 158
5 166
6 178
7 189
8 193
9 197.5
10 200

so there are 11 layers, including the extra 1-metre layers below 140 m and above 200 m. show_bath plots the grid correctly, now that I've fixed the alignment:

To continue troubleshooting, I added lines "print(deltaz)" and "print(zu)" in the file ncoutput.py, and this is the output:

[ 1.   2.5  4.5  4.  11.  12.   8.   5.   6.   7.   1. ]
[201.0 194.0 188.0 183.0 175.0 163.0 152.0 148.0 143.5 141.0 140.0]

So, "deltaz" contains the layer heights from bath.nc, topmost layer first, but "zu" is a list of where the bottom elevations of the layers would be if the order of "deltaz" were from bottom upwards and if the bottom of the lowest layer were at 140 m.

In addition to the obviously incorrect order of layer heights, I wonder if I should have set the bottom_elevation value at 139 m instead on 140 m, not to the lowest layer boundary that I defined myself but to the actual bottom of the model grid, because that way the first and last two values in "zu" would at least coincide with the correct layer boundaries.

Regards,

  Elias

[Dr. Mike Mueller]

I tried to e-create your problem. With attached bathymetry I get the attached
temperature. Looks like expected. The thick layers create considerable
dispersion. In general I would recommend to keep the upper 10 m in 1-m layers.
Are you sure you copied the generated bath.nc to input/w2?


Am 31.05.23 um 11:03 schrieb Elias Pentti:

> Yes, I did update the number_of_layers and number_of_constituent_layers
> values in w2.yaml so that they correspond to my bathymetry grid, as well as
> bottom_elevation in waterbody_coordinates to equal the lowest layer
> boundary defined in reservoir.txt just like in the tutorial. To be clear, my
> reservoir.txt reads Layers ZU 1140 2147 3153 4158 5166 6178 7189 8193

> 9197.5 10200 so there are 11 layers, including the extra 1-metre layers

> below 140 m and above 200 m. show_bath plots the grid correctly, now that

> I've fixed the alignment: uneven_lakegrid.png To continue troubleshooting, I

> added lines "print(deltaz)" and "print(zu)" in the file ncoutput.py, and
> this is the output: [ 1. 2.5 4.5 4. 11. 12. 8. 5. 6. 7. 1.
> ] [201.0 194.0 188.0 183.0 175.0 163.0 152.0 148.0 143.5 141.0 140.0] So,
> "deltaz" contains the layer heights from bath.nc, topmost layer first, but
> "zu" is a list of where the bottom elevations of the layers would be if the
> order of "deltaz" were from bottom upwards and if the bottom of the lowest
> layer were at 140 m. In addition to the obviously incorrect order of layer
> heights, I wonder if I should have set the bottom_elevation value at 139 m
> instead on 140 m, not to the lowest layer boundary that I defined myself but
> to the actual bottom of the model grid, because that way the first and last
> two values in "zu" would at least coincide with the correct layer
> boundaries.
>
> Regards, Elias
>
> tiistai 30. toukokuuta 2023 klo 20.41.00 UTC+3 mmueller kirjoitti:
>
> Hi Elias,
>
> Looks like you changed the number of layers. Did adjust the values in
> input/w2/w2.yaml accordingly? These are the ones from the tutorial model:
>
> general: titel: value: A simple PITLAKQ test model bounds:
> number_of_segments: value: 10 number_of_layers: value: 52
> number_of_constituents_segments: value: 10 number_of_constituents_layers:
> value: 52 times: start: value: 01.01.1998 end: value: 31.01.2009 --- !File

> name: bath.nc <http://bath.nc> format: netcdf group_name: bathymetry ---

> bathymetry: starting_water_level: value: 195.0 --- branch_geometry:
> branch_name: value: - test lake # no w2 variable branch_upstream_segments:
> value: - 2 branch_downstream_segments: value: - 8
>
> Many values relate to the ones described in detail in the CE-QUL-W2 manual:
> https://gitlab.com/hydrocomputing/pitlakq/-/blob/master/doc/w2v2-manual.pdf
> <https://gitlab.com/hydrocomputing/pitlakq/-/blob/master/doc/w2v2-manual.pdf>
>
> Best, Mike
>
> Am 30.05.23 um 15:13 schrieb Elias Pentti:
>> Hello!
>>
>> Now I got to the first simple model runs. First I had to do the fix to
>> numpy.zeros in the file ncoutput.py that you gave last month to another
>> user, because I encountered the exact same error. Now, the model runs,
>> but when I follow the tutorial to visualize the resulting development of
>> temperature, I get unexpected results again. Below is the top part of the
>> plot I get. The upper edge of the white gap in the middle corresponds to
> the
>> calculated lake level, which gradually rises from 195 m to over 196 m
>> because I'm only considering precipitation. In reservoir.txt, my layer
>> boundaries are at 140, 145, ... , 195, 200 m with even 5 m intervals, but

>> when I read the "zu" variable from out.nc <http://out.nc> (supposed to

> be the bottom
>> elevations of the layers), it's {201, 196, 191, 186, 181, 176, 171, 166,
>> 161, 156, 151, 146, 141, 140}. T_result_top.png After this curious
> result, I
>> experimented a little and made the layer structure uneven, setting the
>> layers at 140, 147, 153, 158, 166, 178, 189, 193, 197.5 and 200 m. This
> time

>> "zu" in out.nc <http://out.nc> is {201, 194, 188, 183, 175, 163, 152,

>> <http://bath.nc <http://bath.nc> <http://bath.nc <http://bath.nc>>>",

>> <http://bath.nc <http://bath.nc>>> <http://bath.nc <http://bath.nc>

> <https://groups.google.com/d/msgid/pitlakq-users/ffba33cf-cc75-4364-af44-8916f5baca65n%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pitlakq-users/ffba33cf-cc75-4364-af44-8916f5baca65n%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
> -- You received this message because you are subscribed to the Google
> Groups "pitlakq-users" group. To unsubscribe from this group and stop
> receiving emails from it, send an email to
> pitlakq-user...@googlegroups.com
> <mailto:pitlakq-user...@googlegroups.com>. To view this discussion
> on the web visit

> https://groups.google.com/d/msgid/pitlakq-users/7d852ef4-8a07-41ca-af1b-43fc3c513238n%40googlegroups.com
> <https://groups.google.com/d/msgid/pitlakq-users/7d852ef4-8a07-41ca-af1b-43fc3c513238n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Elias Pentti]

Hello!

The plots you created do look fine, but did you check the "zu" list in out.nc? Because I'm sure that the program does not generate it correctly, the error becoming evident when layer heights are uneven. (Of course I'm going to use much thinner layers to improve accuracy once I'm ready to start actual modelling, but first I want to ensure that I understand how the program works and to be convinced that it does exactly what it should.) Despite my very basic knowledge of Python, I think I was able to fix the error with the following corrections to ncoutput.py:

1. Remove (or comment away) the line (#170 in my version of the file)

        bottom_index += 1

2. On lines #174 & #175, change

        zu[:bottom_index] = (numpy.add.accumulate(deltaz[:bottom_index])[::-1]
                             + bottom)
 into (differing parts in bold)

        zu[:bottom_index] = (numpy.add.accumulate(deltaz[bottom_index:0:-1])[::-1]
                             + bottom)

This way, the cumulative summation of layer heights in deltaz starts from the lake bottom upwards, not from the top of the model downwards as in the original code. Then, bottom_elevation value in waterbody_coordinates of w2.yaml must be set at the bottom elevation of the deepest active cell actually within the lake (not the bottom elevation of the lake grid if the lake does not reach it) in order to make the resulting "zu" equal to the user-defined layer bottom elevations.

The correction to the line between these two (#173) that you posted earlier in another thread, the correct line being

        zu = numpy.zeros_like(output.variables['zu'])

must also be done, at least in my environment.

With best regards,

    Elias

[Dr. Mike Mueller]

Hi Elias,

Thanks for being persistent. There were a few issues with non-equidistant layer
thicknesses. I published a new version 1.7.2. Please install this. It now
supports Python 3.9, 3.10, and 3.11.

I recommend to create a new environment, for example with Python 3.10:

conda create -n pitlakq_py310 python=3.10

(You can also use mamba instead of conda if you work with mamba.)

Activate:

conda activate pitlakq_py310

and install:

conda install pitlakq

You should update your .pitlakq accordingly for the templates and resources
paths.

Please re-do the preprocessing, creating the bath.nc. The recommendation for
the bottom level might be different. Please use the new bath.nc for the
calculation. Let me know if it works for or if you still have issues.

Best,
Mike

Am 05.06.23 um 13:15 schrieb Elias Pentti:

> Hello! The plots you created do look fine, but did you check the "zu" list
> in out.nc? Because I'm sure that the program does not generate it correctly,
> the error becoming evident when layer heights are uneven. (Of course I'm
> going to use much thinner layers to improve accuracy once I'm ready to start
> actual modelling, but first I want to ensure that I understand how the
> program works and to be convinced that it does exactly what it should.)
> Despite my very basic knowledge of Python, I think I was able to fix the
> error with the following corrections to ncoutput.py:
>
> 1. Remove (or comment away) the line (#170 in my version of the file)
> bottom_index += 1
>
> 2. On lines #174 & #175, change zu[:bottom_index] =

> (numpy.add.accumulate(deltaz[*:bottom_index*])[::-1] + bottom) into

> (differing parts in bold) zu[:bottom_index] =

> (numpy.add.accumulate(deltaz[*bottom_index:0:-1*])[::-1] + bottom) This way,

> the cumulative summation of layer heights in deltaz starts from the lake
> bottom upwards, not from the top of the model downwards as in the original
> code. Then, bottom_elevation value in waterbody_coordinates of w2.yaml must
> be set at the bottom elevation of the deepest active cell actually within
> the lake (not the bottom elevation of the lake grid if the lake does not
> reach it) in order to make the resulting "zu" equal to the user-defined
> layer bottom elevations. The correction to the line between these two (#173)
> that you posted earlier in another thread, the correct line being zu =
> numpy.zeros_like(output.variables['zu']) must also be done, at least in my
> environment.
>
> With best regards, Elias keskiviikko 31. toukokuuta 2023 klo 17.12.19 UTC+3
> mmueller kirjoitti:
>
> I tried to e-create your problem. With attached bathymetry I get the
> attached temperature. Looks like expected. The thick layers create
> considerable dispersion. In general I would recommend to keep the upper 10 m
> in 1-m layers. Are you sure you copied the generated bath.nc

> <http://bath.nc> to input/w2?

>
>
> Am 31.05.23 um 11:03 schrieb Elias Pentti:
>> Yes, I did update the number_of_layers and number_of_constituent_layers
>> values in w2.yaml so that they correspond to my bathymetry grid, as well
>> as bottom_elevation in waterbody_coordinates to equal the lowest layer
>> boundary defined in reservoir.txt just like in the tutorial. To be
> clear, my
>> reservoir.txt reads Layers ZU 1140 2147 3153 4158 5166 6178 7189 8193
>> 9197.5 10200 so there are 11 layers, including the extra 1-metre layers
>> below 140 m and above 200 m. show_bath plots the grid correctly, now that
>> I've fixed the alignment: uneven_lakegrid.png To continue
> troubleshooting, I
>> added lines "print(deltaz)" and "print(zu)" in the file ncoutput.py, and
>> this is the output: [ 1. 2.5 4.5 4. 11. 12. 8. 5. 6. 7. 1. ] [201.0 194.0
>> 188.0 183.0 175.0 163.0 152.0 148.0 143.5 141.0 140.0] So, "deltaz"

>> contains the layer heights from bath.nc <http://bath.nc>,

> topmost layer first, but
>> "zu" is a list of where the bottom elevations of the layers would be if
>> the order of "deltaz" were from bottom upwards and if the bottom of the
>> lowest layer were at 140 m. In addition to the obviously incorrect order
>> of layer heights, I wonder if I should have set the bottom_elevation value
>> at 139 m instead on 140 m, not to the lowest layer boundary that I defined
>> myself
> but
>> to the actual bottom of the model grid, because that way the first and
>> last two values in "zu" would at least coincide with the correct layer
>> boundaries.
>>
>> Regards, Elias
>>
>> tiistai 30. toukokuuta 2023 klo 20.41.00 UTC+3 mmueller kirjoitti:
>>
>> Hi Elias,
>>
>> Looks like you changed the number of layers. Did adjust the values in
>> input/w2/w2.yaml accordingly? These are the ones from the tutorial model:
>>
>> general: titel: value: A simple PITLAKQ test model bounds:
>> number_of_segments: value: 10 number_of_layers: value: 52
>> number_of_constituents_segments: value: 10 number_of_constituents_layers:
>> value: 52 times: start: value: 01.01.1998 end: value: 31.01.2009 ---

>> !File name: bath.nc <http://bath.nc> <http://bath.nc <http://bath.nc>>

> <http://out.nc <http://out.nc>> (supposed to

>> be the bottom
>>> elevations of the layers), it's {201, 196, 191, 186, 181, 176, 171,
>>> 166, 161, 156, 151, 146, 141, 140}. T_result_top.png After this curious
>> result, I
>>> experimented a little and made the layer structure uneven, setting the
>>> layers at 140, 147, 153, 158, 166, 178, 189, 193, 197.5 and 200 m. This
>> time

>>> "zu" in out.nc <http://out.nc> <http://out.nc <http://out.nc>> is {201,

> <https://groups.google.com/d/msgid/pitlakq-users/7d852ef4-8a07-41ca-af1b-43fc3c513238n%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pitlakq-users/7d852ef4-8a07-41ca-af1b-43fc3c513238n%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
> -- You received this message because you are subscribed to the Google
> Groups "pitlakq-users" group. To unsubscribe from this group and stop
> receiving emails from it, send an email to
> pitlakq-user...@googlegroups.com
> <mailto:pitlakq-user...@googlegroups.com>. To view this discussion
> on the web visit

> https://groups.google.com/d/msgid/pitlakq-users/1bcf6575-6907-4275-8028-e76669f457fbn%40googlegroups.com
> <https://groups.google.com/d/msgid/pitlakq-users/1bcf6575-6907-4275-8028-e76669f457fbn%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Elias Pentti]

Dear Mike,

The updated version works fine, thanks for your effort!

I moved on to the next challenge, adding branch outflow to the model. As there are definitions for the branch outflow (and inflow) files in the w2.yaml template, but no group for defining outlets at the downstream section, I added the following group
outlets:
    number_of_outlets:
        default:
            - 1
        w2_code_name: nout
        dimension:
            - number_of_branches
    outlet_layer_location:
        default:
            - 2
        w2_code_name: kout
        dimension:
            - number_of_outlet_structures

in the template, and the corresponding group

outlets:
    number_of_outlets:
        value:
            - 1
    outlet_layer_location:
        value:
            -  7

in the model's w2.yaml file. I also created an input file for branch outflow based on the template branch_outflow.txt and included it in the model file. When I run the model, it terminates at the following error message:

Error: W2 variable kout
Error: variable outlet_layer_location
[7]
[[0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]
 [0]]
has dimensions of (1,) but required are dimensions of (18, 1).

So, the program wants "outlet_layer_location" to be a table with 18 rows and one column (one row per model layer, apparently), although I had specified that its dimension should be the number of outlet structures (one in this case), as it should be according to the CE-QUAL-W2 documentation. If I change the input in the "outlets" group into a table with 18 elements as

    outlet_layer_location:
        value:
            - [ 7 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]
            - [ 0 ]

the model runs successfully and produces sensible results with the lake surface elevation changing in accordance with the tributary inflow and branch outflow I had defined. However, this kind of a workaround is not what one should be forced to use. Could you tell how to properly include branch outflow through an outlet or several outlets, in a PITLAKQ model?

I would also like to know if it's possible to use the selective withdrawal feature of CE-QUAL-W2. I have tried to define the necessary keywords and variables for that also, but at least the SWC parameter turned out to be problematic. It should be assigned values of ON or OFF for each branch of the model to switch selective withdrawal on or off in them. Or is it intended that surface outflow from a pit lake or other similar system through a spillway or discharge channel is modelled with the groundwater inflow functionality, defining a water level-dependent groundwater outflow to some cells close to ground level to imitate discharge on the surface? When modelling the formation of a pit lake, it would naturally be preferable to have a way to have an outflow that depends on the water level, not a predefined time series of branch outflow.

Elias

[Dr. Mike Mueller]

Dear Elias,

Am 09.06.23 um 14:59 schrieb Elias Pentti:

> Dear Mike,
>
> The updated version works fine, thanks for your effort!

Good to hear.

Regarding your use of branch outflow. It is not really implemented.
The files in the directory `templates` for branch in- and outflow are needed
by CE-QUAL-W2 but their values are never used. PITLAKQ just copies them
into internal work directory.

You can use a tributary with negative inflow instead. For example,
add these to your `w2.yaml`:

tributaries:
number_of_tributaries:
value: 2
tributary_names:
value:
- mytribname
- outflow
tributary_segment:
value:
- 5
- 7
tributary_inflow_placement:
value:
- 1 # density
- 0 # distributed
tributary_inflow_top_elevation:
value:
- 20.0
- 20.0
tributary_inflow_bottom_elevation:
value:
- 18.0
- 18.0


These are the options for `tributary_inflow_placement`:
distributed: 0
density: 1
specify: 2

Even though `distributed` and `density` don't use
`tributary_inflow_top_elevation` or `tributary_inflow_bottom_elevation`,
you still need to specify them to make dimensions correct.

The file `outflow_inflow.txt` has a negative inflow, i.e. an outflow:

date time tributary_inflow
01.01.1997 00:00 -0.0025
01.01.2011 00:00 0.0

The files `outflow_temperature.txt` and `outflow_concentration.txt` need to be
present, even though their values are not used.

Selective withdrawal is designed for reservoirs with a dam. It is not
implemented in PITLAKQ. SWC is always false. But you maybe able to get
what you want with a tributary with negative inflow and a
`tributary_inflow_placement` set to `specify` (2) with according values for
`tributary_inflow_top_elevation` and `tributary_inflow_bottom_elevation`.

Best,
Mike

>
> I moved on to the next challenge, adding*branch outflow* to the model. As

> possible to use the *selective withdrawal* feature of CE-QUAL-W2. I have

> tried to define the necessary keywords and variables for that also, but at
> least the SWC parameter turned out to be problematic. It should be assigned
> values of ON or OFF for each branch of the model to switch selective
> withdrawal on or off in them. Or is it intended that surface outflow from a
> pit lake or other similar system through a spillway or discharge channel is
> modelled with the groundwater inflow functionality, defining a water
> level-dependent groundwater outflow to some cells close to ground level to
> imitate discharge on the surface? When modelling the formation of a pit
> lake, it would naturally be preferable to have a way to have an outflow that
> depends on the water level, not a predefined time series of branch outflow.
>
> Elias
>

> -- You received this message because you are subscribed to the Google Groups
> "pitlakq-users" group. To unsubscribe from this group and stop receiving
> emails from it, send an email to pitlakq-user...@googlegroups.com
> <mailto:pitlakq-user...@googlegroups.com>. To view this
> discussion on the web visit

> https://groups.google.com/d/msgid/pitlakq-users/5ed0cac7-4909-4ce1-8bdf-657aa2a1c33en%40googlegroups.com
> <https://groups.google.com/d/msgid/pitlakq-users/5ed0cac7-4909-4ce1-8bdf-657aa2a1c33en%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Elias Pentti]

Dear Mike,

Thanks for your advice. I guess it's possible to use tributaries and negative inflow also, the layers just need to be dense enough because the "bottom_elevation" seems to be used only to choose which (active) layers the tributary removes water from, not so that the simulated lake level should actually be above the defined value.

Now I'm finally getting to the water quality calculations, but ran into problems again. I have two tributaries in my training model, one for inflow and one for outflow. I've defined inflow concentrations according to the tutorial, and also edited the files pitlakq.yaml and w2.yaml accordingly to enable phreeqc modelling.

The tutorial instructs to

"start our run with pitlakq -p pitlakq_tut_qual", but that command won't work (using the name of my own model of course). Pitlakq does not recognize the "-p" switch, and does nothing but prints an error message.

If I run the simulation with "pitlakq run modelname" as in the previous phases, it starts and I get messages referring to phreeqc. However, the resulting out.nc file does not contain any water quality results, in particular no pH so that tutorial's "do_show_ph" script cannot plot anything. If I look at the files in the "output\balance" folder, both chlorid_change.txt and sodium_change.txt suggest that the inflow tributary has not caused any change at all (all zeros in the corresponding column), although the tributary has positive concentrations for both species and the inflow has had the same effect on the lake level as before. The columns for the change related to the outflow tributary, however, show some negative figures. What am I doing wrong this time?

In addition, I would like to ask, just to be sure, which way does PITLAKQ handle the sign of longitude in the waterbody coordinates? I noticed that in CE-QUAL-W2, west longitude is positive and east is negative, probably because of the American origin of the program, but in your tutorial, you seem to be using the coordinates of a lake in western Australia if the positive longitude is interpreted to be east of Greenwich, not west.

Best regards,

Elias

[Elias Pentti]

I found out by myself why most water quality parameters were not treated and recorded. I hadn't updated activeconst.txt, but had used the version from the template folder. Now I have ones where they are supposed to be.

[Dr. Mike Mueller]

Dear Elias,

Am 16.06.23 um 09:09 schrieb Elias Pentti:
> Dear Mike,
>

> Thanks for your advice. I guess it's possible to use tributaries and
> negative inflow also, the layers just need to be dense enough because the
> "bottom_elevation" seems to be used only to choose which (active) layers the
> tributary removes water from, not so that the simulated lake level should
> actually be above the defined value.

I am not sure if I understand what you mean. If the lake water level is below
a layer that is in your specified range, no water will be taken from this
layer because there is none. So the simulated lake level is above the defined
value.

...

> The tutorial instructs to "start our run with pitlakq -p pitlakq_tut_qual",
> but that command won't work (using the name of my own model of course).
> Pitlakq does not recognize the "-p" switch, and does nothing but prints an
> error message.

Ups, that is a typo in the documentation. Correct:
pitlakq run pitlakq_tut_qual
Likely that it did not get updated in the past. Will fix.

> In addition, I would like to ask, just to be sure, which way does PITLAKQ
> handle the sign of longitude in the waterbody coordinates? I noticed that in
> CE-QUAL-W2, west longitude is positive and east is negative, probably
> because of the American origin of the program, but in your tutorial, you
> seem to be using the coordinates of a lake in western Australia if the
> positive longitude is interpreted to be east of Greenwich, not west.
>

Can you point me to source for the information that CE-QUAL-W2 treats the
west longitude as positive and east as negative.

Thanks.

Best,
Mike

> <https://groups.google.com/d/msgid/pitlakq-users/5ed0cac7-4909-4ce1-8bdf-657aa2a1c33en%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pitlakq-users/5ed0cac7-4909-4ce1-8bdf-657aa2a1c33en%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
> -- You received this message because you are subscribed to the Google
> Groups "pitlakq-users" group. To unsubscribe from this group and stop
> receiving emails from it, send an email to
> pitlakq-user...@googlegroups.com
> <mailto:pitlakq-user...@googlegroups.com>. To view this discussion
> on the web visit

> https://groups.google.com/d/msgid/pitlakq-users/a65db292-cd67-459b-81c1-d60ee64189ecn%40googlegroups.com
> <https://groups.google.com/d/msgid/pitlakq-users/a65db292-cd67-459b-81c1-d60ee64189ecn%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Elias Pentti]

User Manual for CE-QUAL-W2 version 2.0 from 1995 does not define the sign convention of longitude (or at least I can't find it there), but in the example Control Input File on page C185 of Appendix C (page 313 of the pdf file provided), the coordinates of De Gray Reservoir are LAT 34.2 and LONG 93.3, both positive. The reservoir is in the USA, so the positive longitude must be to the west.

On page 20 in Part 3 (Input and Output files) of version 4.5 User Manual it reads:

"The model was set up for LONG W and LAT N coordinates as being positive. Hence one would start at Greenwich as 0 and go 360 degrees toward the West or go negative toward the east"

[Elias Pentti]

Hello again!

I think I ran into some kind of Python version confllict after adding groundwater inflow in my lake model using the "gwh" option. Pitlakq terminates after printing the following error messages:

C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\pitlakq\submodels\gw\gwh.py:107: FutureWarning: In the future `np.bool` will be defined as the corresponding NumPy scalar.
  len(zone_distribution[0]))).astype(numpy.bool)

runtime: 0:00:02.690194
Traceback (most recent call last):
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\Scripts\pitlakq-script.py", line 9, in <module>
    sys.exit(main())
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\pitlakq\metamodel\running\run_pitlakq.py", line 517, in main
    args.func(args)
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\pitlakq\metamodel\running\run_pitlakq.py", line 450, in run_project
    run(args.project_name)
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\pitlakq\metamodel\running\run_pitlakq.py", line 531, in run
    model = Pitlakq(config,
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\pitlakq\metamodel\running\run_pitlakq.py", line 100, in __init__
    self.init_sub_models()
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\pitlakq\metamodel\running\run_pitlakq.py", line 305, in init_sub_models
    self.gwh.read_data()
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\pitlakq\submodels\gw\gwh.py", line 45, in read_data
    zone_mapping, self.array_mapping = self._make_zone_mapping(
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\pitlakq\submodels\gw\gwh.py", line 107, in _make_zone_mapping
    len(zone_distribution[0]))).astype(numpy.bool)
  File "C:\Users\hyh972\AppData\Local\miniconda3\envs\pitlakq_py310\lib\site-packages\numpy\__init__.py", line 305, in __getattr__
    raise AttributeError(__former_attrs__[attr])
AttributeError: module 'numpy' has no attribute 'bool'.
`np.bool` was a deprecated alias for the builtin `bool`. To avoid this error in existing code, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
    https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations. Did you mean: 'bool_'?

Best regards,

   Elias

[Dr. Mike Mueller]

Hi Elias,

Thanks for pointing this out. I fixed it. Please upgrade to the latest
pitlakq version:

conda upgrade pitlakq
or
mamba upgrade pitlakq

Best,
Mike

PS: The fix has been ready for several weeks, but I could not test it on a
Windows system to make sure it works for you. I am on MacOS and was
traveling without accesses to my Windows machine.

Am 12.07.23 um 14:41 schrieb Elias Pentti:

> <https://groups.google.com/d/msgid/pitlakq-users/a65db292-cd67-459b-81c1-d60ee64189ecn%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pitlakq-users/a65db292-cd67-459b-81c1-d60ee64189ecn%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
> -- You received this message because you are subscribed to the Google
> Groups "pitlakq-users" group. To unsubscribe from this group and stop
> receiving emails from it, send an email to
> pitlakq-user...@googlegroups.com
> <mailto:pitlakq-user...@googlegroups.com>. To view this discussion
> on the web visit

> https://groups.google.com/d/msgid/pitlakq-users/9a43d0e0-1b14-47c3-857f-9363943e3279n%40googlegroups.com
> <https://groups.google.com/d/msgid/pitlakq-users/9a43d0e0-1b14-47c3-857f-9363943e3279n%40googlegroups.com?utm_medium=email&utm_source=footer>.