1

i created a Random Point Generator in python using the algorithm described here by @whuber. The algorithm is implemented using Shapely, PostGIS and psycopg2. To boost performance pythons multiprocessing and threading libraries are used. The source code can be found in my GitHub Repository. The code generates a total of 57.990.970 Points within 6388 Polygons in ~3517 sec on my Xeon-1231v3 with 16GB Ram.

Procedure SimpleRandomSample(P:Polygon, N:Integer) {
    U = Sorted list of N independent uniform values between 0 and 1
    Return SRS(P, BoundingBox(P), U)
}

Procedure SRS(P:Polygon, B:Rectangle, U:List) {
    N = Length(U)
    If (N == 0) {Return empty list}
    aP = Area(P)
    If (aP <= 0) {
        Error("Cannot sample degenerate polygons.")
        Return empty list
    }
    t = 2
    If (aP*t < Area(B)) {
        # Cut P into pieces
        If (Width(B) > Height(B)) {
            B1 = Left(B); B2 = Right(B)
        } Else {
            B1 = Bottom(B); B2 = Top(B)
        }
        P1 = Clip(P, B1); P2 = Clip(P, B2)
        K = Search(U, Area(P1) / aP)
        V = Concatenate( SRS(P1, B1, U[1::K]), SRS(P2, B2, U[K+1::N]) )
    } Else {
        # Sample P
        V = empty list
        maxIter = t*N + 5*Sqrt(t*N)
        While(Length(V) < N and maxIter > 0) {
            Decrement maxIter
            Q = RandomPoint(B)
            If (Q In P) {Append Q to V}
        }
       If (Length(V) < N) {
            Error("Too many iterations.")
       }
    }
    Return V
}

But concerning the algorithm I have an open question. It mentions a U = Sorted list of N independent uniform values between 0 and 1 which i could build using numpy.uniform(0, 1, numberOfValues). But then the algorithm will generate weird results, like one part of the polygon will not have any points in it. I work around this issue by weighting the points to be created by the recursive call using k = int(round(n * (p1.area / polygon.area))).

So my question is why this List of uniform values between 0 and 1 should be used in the first place.

Improve this question
4
  • Bill could answer more completely, but I'd think the discrete set of intervals (randomly chosen) would force the distribution to be more complete in less time. Can you quantify these "weird results" you'd expect? Commented Nov 18, 2014 at 21:25
  • Well the weird result appeared once t was set so something smaller than 2. It started splitting the Polygon and basicly one part of it had no points in it. The problem was that K becomes a value between 0 and 1, but after splitting U like U[1::K] the max value in U might become smaller than K which forces U[K+1::N] to have a length of 0, thus no points are generated in this recursive call. Commented Nov 18, 2014 at 21:32
  • K is an integer between 0 and N. Commented Nov 18, 2014 at 22:09
  • Ehm stupid me, I wanted to say the area of the sub-polygon p1 divided by the whole area is something between 0-1. So Search does look in U where the index of the element is which is smaller than the result. But again if U was split before, lets say at 0.5 and the result is 0.7, K will be the last index of U, thus in p2 no points will be created. It only would make sense to search U if you recreate it every time SRS is run. But then also you would not need to pass it to the function. Just the number of points to be generated would be sufficient. Commented Nov 18, 2014 at 22:24

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.