How do I select a set of uniformly random points from some finite portion of the upper half plane model $$\mathbb{H}= \{ (x ,y) \mid y > 0; x, y \in \mathbb{R} \}$$ with the usual hyperbolic metric? Usually, in flat space, I can select $x,y$ uniformly at random, but now I need a non-uniform density.
From what I can see (in the comments), the upper half plane has a measure of $$\frac{\mathrm{d}x \mathrm{d}y}{y^2}$$ so I can find the area of some region $R$ (such as a discrete cell of the upper-half plane) via $$ A(R) = \iint_{R}\frac{\mathrm{d}x \mathrm{d}y}{y^2} $$
From what I can see here, the heptagons are clearly of different areas as they move away from the lower boundary, so if I want a Poisson point process in the hyperbolic plane, each area of equal Lebesgue measure (e.g. two equal area heptagons) should have the same expected number of points in, so points should be much more dense near the boundary $y=0$.
You did not specify details about the shape of your finite area. For now I'll assume that there exists an axis aligned bounding box in the upper half plane model containing you shape. In that case you can sample uniformly from that rectangle, then repeat sampling if your sample is outside your actual shape. The tighter the bounding box the lower the expected number of samples needed till you get one inside your shape.
So if that axis-aligned rectangle is from $x_1$ to $x_2$ in one direction and $y_1$ to $y_2$ in the other, you want to sample uniformly, with density proportional to the area element
$$\frac{\mathrm dx\,\mathrm dy}{y^2}$$
Since the formula doesn't depend on $x$ you can choose the $x$ coordinate using simple uniform sampling. We need to focus on the $y$ coordinate as the more complicated aspect. Let's just integrate the area element in that single dimension.
$$\int_{y_1}^{y_2}\frac{\mathrm dy}{y^2}=\frac1{y_1}-\frac1{y_2}$$
The cumulative probability function has to be proportional to that. To be more precise, the probability that your random $Y$ coordinate is below a specific threshold $y$ is
$$\operatorname{Pr}(Y\le y)=\frac{\frac1{y_1}-\frac1{y}}{\frac1{y_1}-\frac1{y_2}}$$
Now you can reverse the game. Pick a random value $u$ uniformly from the interval $[0,1]$. Use that as the value of the cumulative distribution function and solve for $y$.
\begin{align*} u&=\frac{\frac1{y_1}-\frac1{y}}{\frac1{y_1}-\frac1{y_2}}\\ u\left(\frac1{y_1}-\frac1{y_2}\right)&=\frac1{y_1}-\frac1{y}\\ \frac1{y_1}-u\left(\frac1{y_1}-\frac1{y_2}\right)&=\frac1{y}\\ \frac1{\frac1{y_1}-u\left(\frac1{y_1}-\frac1{y_2}\right)}&=y \end{align*}
In terms of hyperbolic geometry, the axis-aligned rectangle is actually a pair of concentric horocycles (the horizontal lines) and a pair of limit parallel geodesics orthogonal to these (the vertical lines). Not every finite portion of the hyperbolic plane will allow for such a bounding box. In some cases you can deal with a single ideal point by using $\lim_{y_2\to\infty}$ which makes all those $\frac1{y_2}$ simply disappear. This drops the inner horocycle and works as long as the bounded $x$ range still applies. Other shapes might require more elaborate approaches to avoid dealing with infinite areas.