Shortlisting last-mile depot sites in Sydney

cat > /workspace/scratch/run_shortlist.py <<'EOF'
import numpy as np, pyarrow as pa, zn
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import dijkstra
con=zn.lake_con()
PP=zn.read('parcel_profile')
# Candidate pool: user criteria before catchment threshold
where="""zone_code IN ('E3','E4') AND area_m2 BETWEEN 3000 AND 10000 AND lot <> ''
AND minutes_truck_to_syd_airport <= 30 AND minutes_truck_to_motorway < 5
AND residential_within_200m IS FALSE AND flood_planning IS NOT TRUE
AND clm_record_site_in_cell = 0"""
cand=con.execute(f"SELECT parcel_id, lot, plan_label, suburb, lga_name, postcode, zone_code, zone_name, area_m2, lat, lng, h3_r9, minutes_truck_to_syd_airport, minutes_truck_to_motorway, minutes_truck_to_kfr, flood_planning, freight_minutes_from_neighbour_cell FROM {PP} WHERE {where} ORDER BY parcel_id").fetchall()
cols=['parcel_id','lot','plan_label','suburb','lga_name','postcode','zone_code','zone_name','area_m2','lat','lng','h3_r9','minutes_truck_to_syd_airport','minutes_truck_to_motorway','minutes_truck_to_kfr','flood_planning','freight_minutes_from_neighbour_cell']
print('candidate_pool',len(cand))
# ordered graph nodes, car profile for population catchment
n=con.execute(f"SELECT node_id, lat, lng, h3_r9 FROM {zn.read('drive_time',file='graph_nodes')} ORDER BY node_id").fetch_arrow_table()
e=con.execute(f"SELECT u,v,minutes_car AS w FROM {zn.read('drive_time',file='graph_edges')}").fetch_arrow_table()
assert int(n['node_id'][-1].as_py())==len(n)-1
u=e['u'].to_numpy(zero_copy_only=False); v=e['v'].to_numpy(zero_copy_only=False); w=e['w'].to_numpy(zero_copy_only=False)
graph=coo_matrix((np.r_[w,w],(np.r_[u,v],np.r_[v,u])),shape=(len(n),len(n))).tocsr()
lat=n['lat'].to_numpy(zero_copy_only=False); lng=n['lng'].to_numpy(zero_copy_only=False); cells=np.array(n['h3_r9'].to_pylist())
starts=[]; snaps=[]
for r in cand:
    la=float(r[9]); lo=float(r[10]); c=np.cos(np.radians(la)); i=int(np.argmin((lat-la)**2+((lng-lo)*c)**2)); starts.append(i); snaps.append(111320*float(np.hypot(lat[i]-la,(lng[i]-lo)*c)))
rows=[]
for j in range(0,len(starts),40):
    dist=dijkstra(graph,directed=True,indices=starts[j:j+40],limit=30)
    for r,di,snap in zip(cand[j:j+40],dist,snaps[j:j+40]):
        reach=np.unique(cells[np.isfinite(di)])
        rows.append((r[0],len(reach),round(snap),reach))
# register unique entity cell list and metadata
con.register('catchment',pa.table({'parcel_id':pa.array([pid for pid,_,_,cs in rows for _ in cs]),'h3_r9':pa.array([c for _,_,_,cs in rows for c in cs])}))
con.register('cand_meta',pa.table({
 'parcel_id':pa.array([r[0] for r in cand]),
 'catchment_cells':pa.array([x[1] for x in rows],type=pa.int64()),
 'snap_m':pa.array([x[2] for x in rows],type=pa.int64())}))
POP=zn.read('population',file='population_h3r9')
# full scored pool and thresholded top20, then join profile, roster and sales
# sales parser follows legal lot/plan normalization, latest date then highest price same date
S=zn.read('property_sales')
R=zn.read('hazard_flood_extents',file='flood_extents_roster')
sql=f"""
WITH scored AS (
 SELECT m.parcel_id,m.catchment_cells,m.snap_m,round(sum(p.persons)) AS catchment_persons
 FROM cand_meta m LEFT JOIN catchment c USING(parcel_id) LEFT JOIN {POP} p USING(h3_r9)
 GROUP BY 1,2,3
), ranked AS (
 SELECT *, row_number() OVER (ORDER BY catchment_persons DESC, parcel_id) AS rank
 FROM scored WHERE catchment_persons >= 2000000
), tokens AS (
 SELECT sale_id,sale_date,price,unnest(str_split(regexp_replace(upper(trim(dp_lot)), '\\s*,\\s*', ',', 'g'),' ')) AS token
 FROM {S} WHERE dp_lot IS NOT NULL AND trim(dp_lot)<>'' AND sale_date IS NOT NULL
), split AS (SELECT sale_id,sale_date,price,str_split(token,'/') AS piece FROM tokens WHERE token<>''), shaped AS (
 SELECT sale_id,sale_date,price,CASE WHEN len(piece) IN (2,3) THEN piece[1] END AS lots,CASE WHEN len(piece)=3 THEN piece[2] ELSE '' END AS section,CASE WHEN len(piece) IN (2,3) THEN piece[-1] END AS plan FROM split), expanded AS (
 SELECT sale_id,sale_date,price,unnest(str_split(lots,',')) AS lot,section,plan FROM shaped WHERE lots IS NOT NULL
), matched AS (
 SELECT DISTINCT sale_id,sale_date,price,lot||'/'||section||'/'||CASE WHEN regexp_matches(plan,'^[0-9]+$') THEN 'DP'||plan ELSE plan END AS parcel_id
 FROM expanded WHERE regexp_matches(lot,'^[0-9]+[A-Z]?$|^[A-Z]$') AND regexp_matches(section,'^[0-9A-Z]*$') AND regexp_matches(plan,'^(DP|SP)?[0-9]+$')
), agg AS (
 SELECT parcel_id,max(sale_date) AS last_sale_date,arg_max(price,(sale_date,price)) AS last_sale_price,count(*) AS sales_n FROM matched GROUP BY 1
), last_sale AS (
 SELECT a.parcel_id,a.last_sale_date,a.last_sale_price,a.sales_n,count(*) AS same_day_sales FROM agg a JOIN matched m ON m.parcel_id=a.parcel_id AND m.sale_date=a.last_sale_date GROUP BY 1,2,3,4
), out AS (
 SELECT r.rank,r.parcel_id,p.lot,p.plan_label,p.suburb,p.lga_name,p.postcode,p.zone_code,p.zone_name,round(p.area_m2) AS area_m2,p.lat,p.lng,
 round(p.minutes_truck_to_syd_airport,1) AS truck_min_sydney_airport,round(p.minutes_truck_to_motorway,1) AS truck_min_motorway_ramp,round(p.minutes_truck_to_kfr,1) AS truck_min_key_freight_route,
 r.catchment_persons,r.catchment_cells,r.snap_m, p.flood_planning,
 CASE WHEN p.flood_planning IS FALSE THEN 'assessed: outside' ELSE coalesce('unassessed: EPI layer does not assess this LGA ('||fr.coverage_status||')','unassessed: EPI layer does not assess this LGA') END AS flood_status,
 ls.last_sale_date,ls.last_sale_price,ls.sales_n,ls.same_day_sales,p.freight_minutes_from_neighbour_cell
 FROM ranked r JOIN {PP} p USING(parcel_id) LEFT JOIN {R} fr ON fr.lga_name=p.lga_name LEFT JOIN last_sale ls USING(parcel_id)
 WHERE r.rank<=20 ORDER BY r.rank
)
SELECT * FROM out
"""
zn.write_out(con,'depot_shortlist_top20',sql)
# summaries
print(con.execute(f"SELECT count(*) AS pool, count(*) FILTER (catchment_persons>=2000000) AS clears FROM (SELECT m.parcel_id,sum(p.persons) AS catchment_persons FROM cand_meta m LEFT JOIN catchment c USING(parcel_id) LEFT JOIN {POP} p USING(h3_r9) GROUP BY 1)").fetchall())
print(con.execute("SELECT flood_status,count(*) FROM read_csv_auto('/workspace/out/depot_shortlist_top20.csv') GROUP BY 1").fetchall())
EOF
python3 /workspace/scratch/run_shortlist.py