EDIT: I think I should emphasize that I have no graphics program for this and am not competent to make one. The diagrams below were done by hand, then scanned on my one-page home scanner as jpegs; those seem to work better on MSE than pdf's. My programs give a good idea how the diagram ought to look, also eliminate simple arithmetic errors; however, a user needs to read some rather cryptic output and then draw the diagram.
ORIGINAL: Not Sage, but I have written several programs either using or helping to draw the river for a Pell form. First, i put four related excerpts at http://zakuski.utsa.edu/~jagy/other.html with prefix indefinite_binary. Second, the book by Conway that introduced this diagram is available at http://www.maths.ed.ac.uk/~aar/papers/conwaysens.pdf and for sale as a real book.
Especially for Pell forms, i have come to prefer a hybrid diagram, one that emphasizes the automorphism group of the form $x^2 - n y^2.$ See recent answer at Proving a solution to a double recurrence is exhaustive and, in fact, many earlier answers.
I can tell you that actually drawing these things is what explains them...Conway deliberately leaves out the automorphisms, he wanted a brief presentation I guess, I really wanted to include that and show how the diagram displays the generator of that group. Also discussed in many number theory books, including my favorite, Buell.
You are welcome to email me, gmail is better (click on my profile and go to the AMS Combined Membership Listings link). I have many diagrams, programs in C++, what have you.
Here is the simpler of two diagrams I did for $x^2 - 8 y^2.$ All I mean by the automorphism group is the single formula
$$ (3x+8y)^2 - 8 (x+3y)^2 = x^2 - 8 y^2, $$ with the evident visual column vector $(3,1)^T$ giving a form value of $1$ and the column vector $(8,3)^T$ directly below it giving a form value of $-8,$ thus replicating the original form.
This is another pretty recent, the very similar $x^2 - 2 y^2,$ where I was emphasizing finding all solutions to $x^2 - 2 y^2 = 7,$ and how there is more than one "orbit" of the automorphism group involved, i.e. every other pair...
Well, why not. One should be aware that the Gauss-Lagrange method of cycles of "reduced" forms is part of the topograph, in fact one such cycle is the exact periodocity of Conway's river. Reduced forms, that is $a x^2 + b xy + c y^2$ with $ac < 0$ and $b > |a+c|,$ occur at what Weissman calls "riverbends," where the action switches sides of the river. Anyway, all the following information is automatically part of the diagram for $x^2 - 13 y^2.$ As a result, the diagram is quite large, it took me two pages. Generate solutions of Quadratic Diophantine Equation
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./Pell 13
0 form 1 6 -4 delta -1
1 form -4 2 3 delta 1
2 form 3 4 -3 delta -1
3 form -3 2 4 delta 1
4 form 4 6 -1 delta -6
5 form -1 6 4 delta 1
6 form 4 2 -3 delta -1
7 form -3 4 3 delta 1
8 form 3 2 -4 delta -1
9 form -4 6 1 delta 6
10 form 1 6 -4
disc 52
Automorph, written on right of Gram matrix:
109 720
180 1189
Pell automorph
649 2340
180 649
Pell unit
649^2 - 13 * 180^2 = 1
=========================================
Pell NEGATIVE
18^2 - 13 * 5^2 = -1
=========================================
4 PRIMITIVE
11^2 - 13 * 3^2 = 4
=========================================
-4 PRIMITIVE
3^2 - 13 * 1^2 = -4
=========================================
