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 cur.execute(""" -- (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; --(3) 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;") conn.commit()
Any tips on how I could loop through each geometry and create my 50 nearest points selection?