Incorrect outputs of `xs_utils.CrossSection.plot_xs_field()` method

[ebus-brgm]

Hello Pierre, I've been testing this new feature recently and it is awesome. Thanks much for this contribution!
However i am encoutering a strange issue using the .plot_xs_field() function
I am using a single-layer, regular grid model. I am plotting the cross section by defining coordinates, e.g.:
xs = CrossSection(mm, cross_section_line=[(292, 340), (292,400) ])
I am plotting the chasim.out ("%SATURAT" parameter), using both plot_xs_field() and add_surface()
i am also plotting the water table ("CHARGE") using add_surface()

When I am plotting a cross section, the plot_xs_field() does not provide the right value per cell. As illustrated below:
I plotted the saturation field (0-100%) and the saturation surface (thus displaying between 0 and 100 on the vertical axis) so that the two can be compared (and should be identical).
we should see some red in the cells (0% sat) where the simulated water table is below or close to the bottom; but it is not the case in the cross-section. However the same field displayed using the .add_surface() function is correct (red line)

I should take the time to look into the code...but if it rings any bells... Many thanks!!

Originally posted by @ebus-brgm in #19 (reply in thread)

[pmatran]

Hi @ebus-brgm 👋

I'm glad that you like this feature and find it useful 🙏. First of all, I'd like to thank you very much for your feedback and for the quality of your request. We understand the problem straight away 'cause you've explained the context and the different functions you're targeting.

So, what is happening here exactly 😦?

The CrossSection instance store the input cross section line as shapely LineString object which will be used to intersect grid cells (referenced by node ids). We already faced an strange issue about that: if the cross section line correspond excatly to the cell edges, the intersection will either fail or give duplicated intersected nodes. We fixed that issue by applying an infinitesimal translation (managed by .offset argument) to the original cross section line. Since the .offset translation is calculated by default as $10$% of the tiniest cell resolution, this will absolutly not affect the intersection node ids whatsoever!

Then, your example is the perfect example of distingish 3 types of intersected nodes in the CrossSection class:

• xsnodes : all cell ids that are intersected by the cross section line in model extension
• valid_xsnodes: all active cell ids that are intersected by cross section line which have a top/bottom value
• plottable_xsnodes: all valid nodes for plotting purpose

As a result, we can affirm that:

➡️ len(xsnodes) ≥ len(valid_xsnodes) ≥ len(plottable_xsnodes)

Unfortunately, the translation made on the original cross section line was only taken into account when constructing the .xs_df object. So, the nodes used to create your graph with .plot_xs_field() are correct but most of them are duplicated!

We have recently proposed a solution to this problem by calculating the translation directly from the cross section line definition and storing it in a new attribute .xsline_translated which must be preferred to .xsline during intersection processes.

We also add the add_translated boolean argument in Crosssection.plot() method to actually plot the translated line on xy-plot if user needs to see it!

These changes have been compiled and proposed in #25. Would you mind to test it on your example and merge the pull request into dev if everything is ok ? 😊

Thanks for all,

@pmatran

[ebus-brgm]

Hello @pmatran, I tested the plot_xs_field (modified from #24), which works flawlessly! thanks much!

I also tested the xs.plot() with the add_translated option but this one I could not make it work. I am running into an error I can't really decode:

I can still merge the pull request as it fixes the .plot_xs_field()?

[pmatran]

Hi @ebus-brgm ,

Thank's for your reply! Next time, remember not to merge a bug tagged branch until all errors/issues have been fixed 😜.

Apparently the problem doesn't seem to come from the newest add_translated option, but from the unary_union of active cell geometries. I'm pretty sure the problem is based on shapely library wich is progressivly passing from 1.x to 2.x ! I already update some issues recently #22 since geometries are not considered iterables anymore and become immutable in newest versions. Since shapely is an optional librairy, we haven't fixed the used version for pymarthe yet and, enven if some changes in the code will require a few changes, shapely is only used in optional and non-priority features so I really think we can afford to move to the most recent stable version 2.0.1 and fix that version in requirements.txt. What do you think @apryet 💡?

@ebus-brgm can you confirm which shapely version you're using please 🙏?

Thank you 👋,

@pmatran

[ebus-brgm]

hi @pmatran, hehe sorry got too excited here
I am using shapely version '1.8.5.post1'

[apryet]

Hello @pmatran and @ebus-brgm,
Thank you for this tricky but constructive discussion on that amazing feature ! I haven't tested it yet, will be back asap.

Regarding the opportunity to include shapely 2.x in the requirements.... So far I remember the installation of this library is difficult on Windows ? If that precludes the installation of PyMarthe for users not necessarily in need for PyMarthe features derived from Shapefly, we could rather make it an optional dependency ?
Cheers,
Alex

[pmatran]

Hi @apryet @ebus-brgm ,

As far as I know, the shapely package installation does not cause any special issues on Windows 👍. The only tricky part is on the shapely.speedups module.

It is a feature designed to improve the performance of geometric operations using specific optimizations according to the platform on which shapely is run. That's in fact a C-extension module. However, there is also true that on Windows the speedups module is automatically available and causes some troubles (especialy with fiona io package).

📝 To sum up, we can include it as an optional dependency as @apryet suggested and force the speedups module to be disabled:

import shapely
shapely.speedups.disable()

As a result, the behavior will be the same on all platforms! 🤩

Let’s keep in touch,

@pmatran

[smagne-brgm]

Hi everyone,

To keep the issue tracker up to date, this issue will be closed soon
unless someone has additional comments or information to add.

Please note that the current development branch of pymarthe now
explicitly requires shapely >= 2.1.2.

Best regards,