how to plot sympy Lines Circles etc.

[Erhy]

Hallo!
Habe Darstellende Geometrie Beispiel mit sympy nachcodiert
und würde die Elemente gerne in einem Plot betrachten.

Wie geht das am komfortabelsten

Danke für Euren Tip
Erhy

[Erhy]

Dass es für solch einfache Probleme keine Standardlösungen gibt, wundert mich.

Die sympy Plot Funktion wird in Python 3 nicht unterstützt. (Plot mit Großem 'P')

Werde versuchen MatPlotLib für 2D zu verwenden,
wo es mit

import matplotlib.lines as lines
from matplotlib import artist

äquivalente Ansätze gibt.

Wenn es mir gelingt, werde ich die Routinen hier oder auch in GitHub bekannt machen.

Erhy

[__blackjack__]

Das ist wohl weniger eine Python 3 Frage sondern eher das Pyglet als Plotting-Backend nicht mehr unterstützt wird (auch wenn da anscheinend noch Code vorhanden ist den keiner mehr pflegt und der gemütlich vor sich hin rottet). Und für das Matplotlib-Backend hat anscheinend noch niemand Code geschrieben der die geometrischen Objekte unterstützt. Wenn das so ein einfaches Problem ist, kannst Du ja eine Lösung dafür schreiben und einen „pull request“ beim `sympy`-Projekt machen. Die freuen sich sicher.

[Erhy]

mein Beispiel, um geometrische 2D Elemente, erstellt mit sympy
mit MatPlotLib zu visualisieren.

Code: Alles auswählen

# PlotSympyGeometryWithMatPlotLibExample.py  by Erhy
import matplotlib.lines as mlines
import matplotlib.pyplot as plt
from matplotlib import patches
import sympy.geometry as geom
from sympy import Rational, pi
# import numpy as np #is numpy used here ?

plotArtists = [] #MatPlotLib artists to plot

plot_xlim = (0, .300)
plot_ylim = (0, .255)

plot_xlimR = ( Rational(str(plot_xlim[0])), Rational(str(plot_xlim[1])) )
plot_ylimR = ( Rational(str(plot_ylim[0])), Rational(str(plot_ylim[1])) )

# plot area for conversion from sympy.geometry elements to MatPlotLib artists
rect0 = geom.Polygon((plot_xlimR[0], plot_ylimR[0]),
                     (plot_xlimR[1], plot_ylimR[0]),
                     (plot_xlimR[1], plot_ylimR[1]),
                     (plot_xlimR[0], plot_ylimR[1]))

# # # conversion from sympy.geometry elements to Matplotlib artists # # #
def addGeomForPlot(sympyZs):

    rect0j = rect0 # plot area

    def twoIntersectionWithBorder( pois ) :
        p1 = [pois[0].args[0], pois[0].args[1]]
        p2 = [pois[1].args[0], pois[1].args[1]]
        resul = mlines.Line2D([float(p1[0]), float(p2[0])],
                               [float(p1[1]), float(p2[1])])
        return resul

    for sympyZ in sympyZs:
        # if type(sympyZ) is geom.line.Line2D :
        typstr = str(type(sympyZ))
        if typstr == "<class 'sympy.geometry.line.Segment2D'>":
            pois = geom.intersection(rect0j, sympyZ)
            if len(pois) == 2:  # two points
                result = twoIntersectionWithBorder( pois )
                result.set_linewidth( 2 )
                plotArtists.append(result)
            elif len(pois) == 1:  # one point
                # find the other point of the segment, which is inside the borders
                if rect0j.encloses_point(sympyZ.args[0]) :
                    poi2 = sympyZ.args[0]
                else:
                    poi2 = sympyZ.args[1]

                p1 = [pois[0].args[0], pois[0].args[1]]
                p2 = [poi2.args[0], poi2.args[1]]
                result = mlines.Line2D([float(p1[0]), float(p2[0])],
                                      [float(p1[1]), float(p2[1])],
                                      linewidth=2)
                plotArtists.append(result)
            else:
                if rect0j.encloses_point(sympyZ.args[0]) and rect0j.encloses_point(sympyZ.args[1]):
                    # inside the borders
                    result = mlines.Line2D([float(sympyZ.args[0].args[0]), float(sympyZ.args[1].args[0])],
                                           [float(sympyZ.args[0].args[1]), float(sympyZ.args[1].args[1])],
                                           linewidth=2)
                    plotArtists.append(result)

        elif typstr == "<class 'sympy.geometry.line.Line2D'>":
            rect0j = rect0
            # pois = rect0j.intersection(sympyZ)
            pois = geom.intersection(rect0j, sympyZ)
            if len(pois) == 2:  # two points
                result = twoIntersectionWithBorder( pois )
                result.set_linewidth( 1 )
                plotArtists.append(result)

        elif typstr == "<class 'sympy.geometry.line.Ray2D'>":
            rect0j = rect0
            pois = geom.intersection(rect0j, sympyZ)
            if len(pois) == 2:  # two points
                result = twoIntersectionWithBorder( pois )
                result.set_linewidth( 1 )
                plotArtists.append(result)
            elif len(pois) == 1:
                p1 = [float(pois[0].args[0]), float(pois[0].args[1])]
                pSource = sympyZ.source
                p2 = [float(pSource.args[0]), float(pSource.args[1])]
                result = mlines.Line2D([p1[0], p2[0]],
                                       [p1[1], p2[1]], linewidth=1)
                plotArtists.append(result)

        elif typstr == "<class 'sympy.geometry.ellipse.Circle'>":
            result = patches.Circle((float(sympyZ.args[0].args[0]), float(sympyZ.args[0].args[1]),
                                     ), float(sympyZ.args[1]), linewidth=1, fill=False)
            plotArtists.append(result)


# # example for visual check of angle computation # #

poiA = [Rational('0.04'), Rational('0.015')]
poiB = [Rational('0.135'), Rational('0.05')]
poiC = [Rational('0.210'), Rational('0.140')]
p1, p2, p3 = geom.Point(poiA), geom.Point(poiB), geom.Point(poiC)

r0 = Rational('0.0') #zero
l0 = geom.Line(geom.Point([r0, r0]), geom.Point([r0 + 1, r0])) # horizontal line

segm1, segm2 = geom.Segment(p1, p2), geom.Segment(p2, p3)
segmMinLen = min([segm1.length, segm2.length])
# circle
ci = geom.Circle(p2, segmMinLen)

# line for visual test:
lMittel = geom.Line(
    segm1.intersection(ci)[0], segm2.intersection(ci)[0]). \
    parallel_line(p2)

# angle computation
angSegm1 = l0.smallest_angle_between(segm1)
angSegm2 = l0.smallest_angle_between(segm2)
angStretch = angSegm1 + ((angSegm2 - angSegm1) / Rational('2.0'))
# degStretch = np.rad2deg(float(angStretch))
forGradAng = pi + angSegm1 + ((angStretch - angSegm1) / Rational('2.0'))
# degGradAng = np.rad2deg(float(forGradAng))
rayForGradient = geom.Ray(p2, angle=forGradAng)

addGeomForPlot([segm1, segm2, ci, lMittel, rayForGradient])

fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(1, 1, 1, aspect=1)

ax.set_xlim(plot_xlim[0], plot_xlim[1])
ax.set_ylim(plot_ylim[0], plot_ylim[1])

for aj in plotArtists: # plot all artists
    ax.add_artist(aj)

plt.show()

[__blackjack__]

@Erhy: Namen werden in Python klein_mit_unterstrichen geschrieben. Ausnahmen sind Konstanten (KOMPLETT_GROSS) und Klassen (MixedCase).

Auf Modulebene gehört nur Code der Konstanten, Funktionen, und Klassen definiert. Das Hauptprogramm steht üblicherweise in einer Funktion die `main()` heisst.

Namen sollten nicht abgekürzt oder nummeriert werden. Wenn man `rect0` gleich `plot_area` nennen würde, bräuchte man keinen Kommentar schreiben was `rect0` eigentlich bedeuten sollte, sondern hätte überall wo der Name verwendet wird die Bedeutung deutlich lesbar im Namen stehen. Ähnliches für `r0` vs. `zero`. Wobei man da auch ein 0 hätte nehmen können. Oder `ci` vs. `circle`. `resul` statt `result`? Tatsächlich um *ein* Zeichen zu ”sparen”?

Funktionen und Methoden sollten alles was sie ausser Konstanten benötigen als Argument(e) übergeben bekommen. Sonst sind sie nicht unabhängig von ihrer Umgebung test- und wiederverwendbar. Auch das verschieben in ein anderes Modul wird so zu einer unnötig aufwändigen Aufgabe.

Man sollte einer Funktion keine Liste übergeben die in der Funktion dann gefüllt wird. Falls man so etwas macht, sollte man sehr deutlich dokumentieren das die Liste in der Funktion verändert wird. Besser ist es die Funktion die Liste erstellen und als Rückgabewert liefern lassen. Oder eine Generatorfunktion zu schreiben, dann kann der Aufrufer entscheiden ob er tatsächlich eine Liste braucht, oder nur etwas zum iterieren.

Die Zuweisungen an `rect0j` machen keinen Sinn. Das führt ohne Not einfach nur noch einen weiteren Namen ein, der zudem noch nicht einmal zur Verständlichkeit beiträgt.

Funktionen sollte man nicht schachteln sofern man nicht tatsächlich ein Closure benötigt, oder die innere Funktion wirklich sehr kurz und simpel ist. Man kann die innere Funktion nicht separat testen und der Leser muss erst einmal schauen ob es sich um ein Closure handelt oder nicht.

`twoIntersectionWithBorder()` ist ein nicht wirklich offensichtlicher Name für eine Funktion die aus zwei Sympy-Punkten eine Matplotlib-Linie erstellt.

`p1` und `p2` lassen sich auch einfacher ermitteln, weil kein Grund besteht die Werte einzeln in neue Listen zu kopieren:

Code: Alles auswählen

    p1 = [points[0].args[0], points[0].args[1]]
    p2 = [points[1].args[0], points[1].args[1]]
    
    # ->
    
    p1 = points[0].args
    p2 = points[1].args

Die gesamte Funktion kann man ziemlich simpel auch so schreiben:

Code: Alles auswählen

def make_2d_line_from_points(points):
    return mlines.Line2D(
        [point.args[0] for point in points],
        [point.args[1] for point in points],
    )

Da direkt nach jedem Aufruf der Funktion die Linienstärke festgelegt wird, sollte man das vielleicht auch gleich mit *in* die Funktion nehmen.

Die Funktion wird auch nicht immer verwendet wenn man sie verwenden könnte, was zu sehr ähnlichem Code führt.

Datentypen mit `type()` und `str()` in Zeichenketten umzuwandeln um die dann zu vergleichen ist falsch. Es gibt keine wirkliche Garantie wie diese Zeichenketten konkret aussehen. Wenn man Typtests macht, dann mit `isinstance()` – was auch Vererbung berücksichtigt. Noch als Anmerkung dazu: Typtests sind hier okay, aber generell ein „code smell“ in objektorientierten Programmiersprachen.

Ein paar von den `args[]`-Zugriffen kann man verständlicher schreiben wenn man den dazugehörigen Attributnamen verwendet.

Fast alles an Werten in Listen umkopieren und in `float()` wandeln ist überflüssig.

An das Ende der ``if``/``elif``-Kaskade gehört ein ``else`` das auf unbekannte Typen reagiert.

Die Konvertierungsfunktion würde ich noch mit `functools.singledispatch` aufteilen und die `main()` ist für meinen Geschmack auch schon zu lang, aber das hier ist der Zwischenstand:

Code: Alles auswählen

#!/usr/bin/env python3
# PlotSympyGeometryWithMatPlotLibExample.py  by Erhy
from matplotlib import patches
from matplotlib import lines as mlines
from matplotlib import pyplot as plt
from sympy import Rational, pi as PI
from sympy import geometry as geom


def make_2d_line_from_points(points, linewidth):
    return mlines.Line2D(
        [point.x for point in points],
        [point.y for point in points],
        linewidth=linewidth,
    )


def convert_geometry_to_artists(plot_area, geometry_entities):
    #
    # TODO Use `functools.singledispatch` to break this up into smaller, self
    #   contained functions, so there's no chance of leaking values from one
    #   type/loop iteration to another.
    #
    for entity in geometry_entities:
        if isinstance(entity, geom.line.Segment2D):
            points = geom.intersection(plot_area, entity)
            if len(points) == 2:
                yield make_2d_line_from_points(points, 2)
            elif len(points) == 1:
                #
                # Find the other point of the segment, which is inside the
                # borders.
                #
                point_b = next(
                    point
                    for point in entity.args
                    if plot_area.encloses_point(point)
                )
                yield make_2d_line_from_points([points[0], point_b], 2)
            else:
                if all(map(plot_area.encloses_point, entity.args)):
                    # Completely inside the borders.
                    yield make_2d_line_from_points(entity.args, linewidth=2)

        elif isinstance(entity, geom.line.Line2D):
            points = geom.intersection(plot_area, entity)
            if len(points) == 2:
                yield make_2d_line_from_points(points, 1)

        elif isinstance(entity, geom.line.Ray2D):
            points = geom.intersection(plot_area, entity)
            if len(points) == 2:
                yield make_2d_line_from_points(points, 1)
            elif len(points) == 1:
                yield make_2d_line_from_points([points[0], entity.source], 1)

        elif isinstance(entity, geom.ellipse.Circle):
            yield patches.Circle(
                entity.center, entity.radius, linewidth=1, fill=False
            )

        else:
            raise ValueError(f"unknown geometry entity {entity!r}")


def main():
    #
    # Plot area for conversion from sympy.geometry elements to MatPlotLib
    # artists.
    #
    plot_max_x = Rational("0.300")
    plot_max_y = Rational("0.255")
    plot_area = geom.Polygon(
        (0, 0), (plot_max_x, 0), (plot_max_x, plot_max_y), (0, plot_max_y)
    )
    #
    # Example for visual check of angle computation.
    #
    point_a = geom.Point([Rational("0.04"), Rational("0.015")])
    point_b = geom.Point([Rational("0.135"), Rational("0.05")])
    point_c = geom.Point([Rational("0.210"), Rational("0.140")])

    horizontal_line = geom.Line(geom.Point([0, 0]), geom.Point([1, 0]))

    segment_a = geom.Segment(point_a, point_b)
    segment_b = geom.Segment(point_b, point_c)
    min_segment_length = min(segment_a.length, segment_b.length)
    circle = geom.Circle(point_b, min_segment_length)

    # Line for visual test.
    middle_line = geom.Line(
        segment_a.intersection(circle)[0], segment_b.intersection(circle)[0]
    ).parallel_line(point_b)

    # Angle computation.
    segment_a_angle = horizontal_line.smallest_angle_between(segment_a)
    segment_b_angle = horizontal_line.smallest_angle_between(segment_b)
    angle_stretch = segment_a_angle + ((segment_b_angle - segment_a_angle) / 2)
    gradient_angle = (
        PI + segment_a_angle + ((angle_stretch - segment_a_angle) / 2)
    )
    gradient_ray = geom.Ray(point_b, angle=gradient_angle)

    plot_artists = convert_geometry_to_artists(
        plot_area, [segment_a, segment_b, circle, middle_line, gradient_ray]
    )

    figure = plt.figure(figsize=(7, 7))
    axis = figure.add_subplot(1, 1, 1, aspect=1)
    axis.set_xlim(0, float(plot_max_x))
    axis.set_ylim(0, float(plot_max_y))
    for artist in plot_artists:
        axis.add_artist(artist)
    plt.show()


if __name__ == "__main__":
    main()

[Erhy]

danke fürs lesen!
Aus den Ratschlägen kann ich viel lernen.
Werde beizeiten ein Update nachreichen.
Erhy

[__blackjack__]

So, die Konvertierungsfunktion mit Hilfe von `functools.singledispatch()` auf eine Funktion pro `GeometryEntity` aufgeteilt, die letzten `args`-Zugriffe durch spezifischere Attributnamen ersetzt, und einen Teil der `main()`-Funktion in eine eigene Funktion herausgezogen:

Code: Alles auswählen

#!/usr/bin/env python3
# PlotSympyGeometryWithMatPlotLibExample.py  by Erhy
from functools import singledispatch

from matplotlib import lines, patches, pyplot as plt
from sympy import geometry, pi as PI, Rational


def make_2d_line_from_points(points, linewidth):
    return lines.Line2D(
        [point.x for point in points],
        [point.y for point in points],
        linewidth=linewidth,
    )


@singledispatch
def convert_to_artist(geometry_entity, _plot_area):
    raise ValueError(f"unknown geometry entity {geometry_entity!r}")


@convert_to_artist.register(geometry.line.Segment2D)
def _convert_segment(segment, plot_area):
    points = plot_area.intersection(segment)
    if len(points) == 2:
        return make_2d_line_from_points(points, 2)
    elif len(points) == 1:
        #
        # Find the other point of the segment, which is inside the
        # borders.
        #
        point_b = next(filter(plot_area.encloses_point, segment.points))
        return make_2d_line_from_points([points[0], point_b], 2)
    elif all(map(plot_area.encloses_point, segment.points)):
        # Completely inside the borders.
        return make_2d_line_from_points(segment.points, linewidth=2)
    else:
        return None


@convert_to_artist.register(geometry.line.Line2D)
def _convert_line(line, plot_area):
    points = geometry.intersection(plot_area, line)
    return make_2d_line_from_points(points, 1) if len(points) == 2 else None


@convert_to_artist.register(geometry.line.Ray2D)
def _convert_ray(ray, plot_area):
    points = plot_area.intersection(ray)
    if len(points) == 2:
        return make_2d_line_from_points(points, 1)
    elif len(points) == 1:
        return make_2d_line_from_points([points[0], ray.source], 1)
    else:
        return None


@convert_to_artist.register(geometry.ellipse.Circle)
def _convert_circle(circle, _plot_area):
    return patches.Circle(
        circle.center, circle.radius, linewidth=1, fill=False
    )


def convert_geometries_to_artists(plot_area, geometry_entities):
    for entity in geometry_entities:
        artist = convert_to_artist(entity, plot_area)
        if artist is not None:
            yield artist


def create_geometry_entities():
    """Example for visual check of angle computation."""

    point_a = geometry.Point([Rational("0.04"), Rational("0.015")])
    point_b = geometry.Point([Rational("0.135"), Rational("0.05")])
    point_c = geometry.Point([Rational("0.210"), Rational("0.140")])

    segment_a = geometry.Segment(point_a, point_b)
    segment_b = geometry.Segment(point_b, point_c)
    min_segment_length = min(segment_a.length, segment_b.length)
    circle = geometry.Circle(point_b, min_segment_length)

    # Line for visual test.
    middle_line = geometry.Line(
        segment_a.intersection(circle)[0], segment_b.intersection(circle)[0]
    ).parallel_line(point_b)

    # Angle computation.
    horizontal_line = geometry.Line(
        geometry.Point([0, 0]), geometry.Point([1, 0])
    )
    segment_a_angle = horizontal_line.smallest_angle_between(segment_a)
    segment_b_angle = horizontal_line.smallest_angle_between(segment_b)
    angle_stretch = segment_a_angle + ((segment_b_angle - segment_a_angle) / 2)
    gradient_angle = (
        PI + segment_a_angle + ((angle_stretch - segment_a_angle) / 2)
    )
    gradient_ray = geometry.Ray(point_b, angle=gradient_angle)

    return [segment_a, segment_b, circle, middle_line, gradient_ray]


def main():
    #
    # Plot area for conversion from sympy.geometry elements to MatPlotLib
    # artists.
    #
    plot_max_x = Rational("0.300")
    plot_max_y = Rational("0.255")
    plot_area = geometry.Polygon(
        (0, 0), (plot_max_x, 0), (plot_max_x, plot_max_y), (0, plot_max_y)
    )

    plot_artists = convert_geometries_to_artists(
        plot_area, create_geometry_entities()
    )

    figure = plt.figure(figsize=(7, 7))
    axis = figure.add_subplot(1, 1, 1, aspect=1)
    axis.set_xlim(0, float(plot_max_x))
    axis.set_ylim(0, float(plot_max_y))
    for artist in plot_artists:
        axis.add_artist(artist)
    plt.show()


if __name__ == "__main__":
    main()

[Erhy]

Danke!
Sehr professionell, wie ich es mit einem halben Jahr Python Erfahrung nicht geschafft hätte.

Wie ich aus deiner Homepage http://blog.marc.rintsch.de/ ersehe, ähneln unsere EDV Lebensläufe.
Mein EDV Kontakt war 1976 mit Micro-Prozessoren.

Ich habe nichts dagegen, wenn du mit dem Code
unter deinen Namen einen „pull request“ beim `sympy`-Projekt machst.

Mir würden dabei sicherlich formale Fehler passieren.

Oder wir überlassen das einem jungen hoffnungsvollen Informatiker.

Erhy

[ggamauf]

...

[ThomasL]

@Erhy Ob ihm das so recht ist, von dir identitätsmäßig geoutet zu werden?

[snafu]

@Erhy Ob ihm das so recht ist, von dir identitätsmäßig geoutet zu werden?

Wie meinst du das? Der Blog in seinem Profil verlinkt. Sehe ich nicht so sehr als Outen, wenn man den auch hier im Thread nennt.

[Erhy]

der Administrator dieses Forums kann gerne meinen letzten Beitrag löschen oder die Webadresse ausblenden
Erhy

[__blackjack__]

Wie snafu schon schrieb ist das kein grosses Geheimnis. Allet jut.