I've seen many example of loops, but I can't get that one to work.

I work for a regional transit agency and each year, we conduct several surveys. I geocode everything to the house level for internal use, but I've been asked to anonymize our results so I could be sent to our partners while preserving the anonymity of the persons who answered.

We have an address point table containing every street addresses on our territory. For each of them, we have the number of housing units.

To do it, I use PostGIS. For each geocoded address, I select the 50 nearest points, and then select base on the weight of the number of units. Basically, it's weighted random n nearest neighbor selection.

If I do my script for a single point, it's perfect. Maybe it's not the most elegant solution, but it works. I've made several tests and the weighted random select from the selected 50 nearest points is good.

The only thing is that I can't do this using a function (didn't say it can't be done, just that I can't ;) )

So, I created a python script, that loops using a serial value field and it's working. The problem is that I'm the only one that uses python so I wanted to do it entirely within PostGIS.

Here's my script:

conn = psycopg2.connect("dbname=anon user=postgres password=######")

#-- Cursor
cur = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)

#-- Add necessary fields
cur.execute("""alter table od2008 add column val serial;
               alter table od2008 add column aqidadr varchar(20);
               alter table od2008 add column xanon numeric;
               alter table od2008 add column yanon numeric;
               alter table od2008 add column distance numeric;
               alter table od2008 add column sm2 numeric;""")

#-- Number of rows
cur.execute("SELECT max(val) as maxval FROM od2008")
rec = cur.fetchone()
nb = rec['maxval']

val = 1
while val <= nb:

    #-- Select row "i"
    cur.execute("SELECT * FROM od2008 WHERE val = %i" %(val))
    rec = cur.fetchone()

    #-- (1) Select the 50 nearest points
    #-- (2) Calculating total weight
    #-- (3) Weighted random selection

    #-- Table od2008 --> to anonymize
    #- Table aq_adr --> Address points to select from

       -- (1)
       CREATE TABLE near50 AS
        SELECT aq_adr.*, ST_Distance(od2008.geom, aq_adr.geom) AS distance
        FROM aq_adr, od2008 
            WHERE ST_DWithin(od2008.geom, aq_adr.geom, 1000) AND
            od2008.val = %i
            AND od2008.sm = aq_adr.sm ORDER BY distance LIMIT 50;

       -- (2)
   CREATE TABLE weight AS 
            SELECT sum(nbunit) AS weight
        FROM near50;

       UPDATE od2008 
            SET aqidadr = (SELECT idadr
                FROM near50, weight
                ORDER BY Random()*(1/(nbunit*100/weight.weight)) LIMIT 1)
                WHERE od2008.val = %i;
   """%(val, val))
    val += 1

    #-- Delete temp tables
    cur.execute("DROP TABLE near50, weight;")


Any tips on how I could loop through each geometry and create my 50 nearest points selection?


  • 1
    I hate to be vague (lack of time), but maybe you could do a bigger join and get all the matches to your points. Then use "window functions" to grab the top sorted address for each point. I would also research "common table expressions" and "cursors". Hope it helps, sorry for handwaving. – forkandwait Aug 12 '13 at 21:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.