Commit d6865cdb authored by Eliot Blennerhassett's avatar Eliot Blennerhassett

Initial code

parent 951055b9
*.osm
*.osc
*.gpkg
*.sqlite
This is to support import of LINZ address database into OSM
Inputs
nz-street-address.gpkg obtained from LINZ
Run xapi-osm-nz-addr.sh to obtain nz_address-osm_raw.sqlite
Run spatialite osmlinzaddr.sqlite < init_osmlinzaddr_spatialite.sql
to initialize working database
Run spatialite osmlinzaddr.sqlite < import-xref.sql
to import the LINZ and OSM datasets, and create a cross reference
of duplicate addresses.
Run osmlinzaddr.py to generate OSM changesets
(currently hardcoded to subdirectory bysuburb)
-- Examples
-- spatialite osmlinzaddr.sqlite < example-queries.sql
-- Get linz addresses with ids that aren't in osm_linz
SELECT
linz_addr.linz_id,
linz_addr.full_address_number as housenumber,
linz_addr.full_road_name as street,
X(linz_addr.geom) as x,
Y(linz_addr.geom) as y,
linz_addr.suburb_locality as s,
linz_addr.town_city as c
FROM linz_addr
LEFT OUTER JOIN osm_linz
on (linz_addr.linz_id = osm_linz.linz_id)
WHERE
osm_linz.linz_id is NULL
and linz_addr.addr_type != "Water"
ORDER BY linz_addr.town_city, linz_addr.suburb_locality
LIMIT 25;
-- exact match on number, street, and close proximity
-- 39047/46858
create table temp.exact as
SELECT
o.osm_id as osm_id,
l.linz_id as linz_id,
Distance(o.geom, l.geom) as D
FROM osm_addr as o
JOIN linz_addr as l on
o.housenumber = l.full_address_number
and ((o.street = l.full_road_name))
and D < 0.001
-- remaining OSM entries with streets
-- 6539/46858
create table temp.inexact as
select type, osm_addr.osm_id, housenumber, street, geom
from osm_addr
left outer join exact
on osm_addr.osm_id = exact.osm_id
where exact.osm_id is null and street is not null
order by street, housenumber
-- and without streets
-- 1272/46858
create table temp.streetnull as
select type, osm_addr.osm_id, housenumber, street, geom
from osm_addr
left outer join exact
on osm_addr.osm_id = exact.osm_id
where exact.osm_id is null and street is null
order by street, housenumber
-- .open "osmlinzaddr.sqlite"
-- as main
.trace ON
ATTACH DATABASE "nz-street-address.gpkg" as "linz";
--SELECT AutoGPKGStart();
-- gpkg doesn't contain this one?!
-- INSERT into linz.spatial_ref_sys SELECT * FROM ol."spatial_ref_sys" WHERE srid = 4326
-- Convert LINZ geometry from NZ2000 to WGS84 SRID 4326
INSERT INTO linz_addr
SELECT
address_id as linz_id,
change_id,
address_type as addr_type,
suburb_locality,
town_city,
full_address_number,
full_road_name,
MakePoint(gd2000_xcoord, gd2000_ycoord, 4326) as geom
FROM linz.nz_street_address;
DETACH linz;
ATTACH DATABASE "nz_address-osm_raw.sqlite" AS "osm";
-- Gathering the OSM addresses as points
INSERT INTO osm_addr
SELECT
'way' as type,
w.way_id as osm_id,
t1.v as housenumber,
t2.v as street,
ST_Centroid(ST_Collect(n.Geometry)) as geom
FROM osm.osm_ways as w
INNER JOIN osm_way_tags as t1 on (w.way_id = t1.way_id and t1.k = "addr:housenumber")
LEFT JOIN osm.osm_way_tags as t2 on (w.way_id = t2.way_id and t2.k = "addr:street")
INNER JOIN osm.osm_way_refs as r on (r.way_id = w.way_id)
INNER JOIN osm.osm_nodes as n on (r.node_id = n.node_id)
GROUP BY w.way_id
;
INSERT INTO osm_addr
SELECT
'node' as type,
n.node_id as osm_id,
t1.v as housenumber,
t2.v as street,
Geometry as geom
FROM osm.osm_nodes as n
INNER JOIN osm.osm_node_tags as t1 on (n.node_id = t1.node_id and t1.k = "addr:housenumber")
LEFT JOIN osm.osm_node_tags as t2 on (n.node_id = t2.node_id and t2.k = "addr:street")
;
DETACH osm;
-- ---------------------------------------------------------------------
-- matching linz vs OSM nodes on housenumber and empirically determined distance
-- (Will take quite a long time!)
--
INSERT INTO osm_linz
SELECT
o.osm_id as osm_id,
l.linz_id as linz_id,
Distance(o.geom, l.geom) as D
FROM osm_addr as o
JOIN linz_addr as l on
o.housenumber = l.full_address_number
-- and (o.street IS NULL or (o.street = l.full_road_name))
and D < 0.001
;
-- initialise osmlinzaddr.sqlite (as spatialite database)
-- Subset of LINZ address data,
CREATE TABLE linz_addr(
linz_id mediumint,
change_id mediumint,
addr_type text(20),
suburb_locality text(80),
town_city text(80),
full_address_number text(100),
full_road_name text(250)
);
SELECT AddGeometryColumn('linz_addr', "geom", 4326, "POINT", "XY", 1);
-- OSM address data
CREATE TABLE osm_addr(
type TEXT,
osm_id INT,
housenumber TEXT,
street TEXT
);
SELECT AddGeometryColumn('osm_addr', "geom", 4326, "POINT", "XY", 1);
-- Cross reference of items representing the same address
CREATE TABLE osm_linz(
'osm_id' int,
'linz_id' int,
'D' real -- distance between OSM and LINZ address (? degrees ?)
);
# Currently python2 because pyspatialite is 2 only
from __future__ import print_function, division, unicode_literals
from lxml import etree as ET
import itertools
import os
import pyspatialite.dbapi2 as db
def generate_osmchange(dbfile, outdir,
changeset = 1,
min_nodes = 1, max_nodes = 8000,
include_location = False):
try:
os.makedirs(outdir)
except OSError:
pass
def add_tag(node, k, v):
ET.SubElement(node, "tag", k=k, v=v)
con = db.connect(dbfile)
cur = con.cursor()
if False:
q = '''pragma table_info("osm_address")'''
rows = cur.execute(q)
print('osm_address: ',list(rows))
q = '''
select
linz_addr.linz_id,
linz_addr.full_address_number as housenumber,
linz_addr.full_road_name as street,
X(linz_addr.geom) as x,
Y(linz_addr.geom) as y,
linz_addr.suburb_locality as s,
linz_addr.town_city as c
from linz_addr
left outer join osm_linz
on (osm_id = linz_addr.linz_id)
where
osm_linz.linz_id is NULL
and linz_addr.addr_type != "Water"
order by linz_addr.town_city, linz_addr.suburb_locality
'''
rows = cur.execute(q)
finished = False
prev_suburb = None
num_nodes = 0
block = 0
osmchange = None
while True:
try:
r = rows.next()
suburb = r[5]
city = r[6]
num_nodes += 1
except StopIteration:
finished = True
if (not osmchange or finished
or ((suburb != prev_suburb) and (num_nodes >= min_nodes))
or (num_nodes >= max_nodes)):
if osmchange:
tree = ET.ElementTree(osmchange)
fn = outdir + "/linz_addr_%04d_%s.osc" % (block, fn)
print("Writing", num_nodes, "to", fn)
tree.write(fn, pretty_print=True)
if finished:
break
block += 1
if city is None or city == suburb:
fn = suburb
else:
fn = city + '-' + suburb
fn = fn.replace(' ', '_')
fn = fn.replace('/', '_')
num_nodes = 1
osmchange = ET.Element('osmChange',
version="0.6",
generator="linz_address.py",
attribution="https://wiki.openstreetmap.org/wiki/Contributors#LINZ",
source="https://data.linz.govt.nz/layer/3353-nz-street-address"
)
if suburb != prev_suburb:
print('city', city, 'suburb/hamlet', suburb)
prev_suburb = suburb
create = ET.SubElement(osmchange, "create")
dnode = {
'id': "-%d" % r[0],
'changeset': '-1',
'lat': "%8.8f" % r[4],
'lon': "%8.8f" % r[3],
'visible': "true",
'version': "1"
}
node = ET.SubElement(create, "node", dnode)
add_tag(node, "addr:housenumber", r[1])
add_tag(node, "addr:street", r[2])
if include_location:
if city is not None:
add_tag(node, "addr:city", city)
add_tag(node, "addr:suburb", suburb)
else:
add_tag(node, "addr:hamlet", suburb)
if finished:
break
if __name__ == '__main__':
generate_osmchange('osmlinzaddr.sqlite', 'bysuburb', min_nodes=4000)
#!/bin/sh
# Get all NZ objects that have an addr:housenumber tag
#wget -O nz_addr.osm "http://www.overpass-api.de/api/xapi_meta?*[addr:housenumber=*][bbox=157.5,-59.0,179.9,-25.5]"
#wget -O chathams_addr.osm "http://www.overpass-api.de/api/xapi_meta?*[addr:housenumber=*][bbox=-177,-44.5,-175.5,-43.5]"
spatialite_osm_raw -d nz_address-osm_raw.sqlite -o nz_addr.osm
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment