1 | #!/bin/env python
|
---|
2 |
|
---|
3 | # Extract Bhutan L3s (Chiwogs) from OSM & write as Sahana import file
|
---|
4 | # Designed to be run in Web2Py context with the BT_L1.csv & BT_L2.csv already imported
|
---|
5 |
|
---|
6 | # L1 (Dzonkhags) are admin_level = '4' in OSM
|
---|
7 | # L2 (Gewogs) are admin_level = '6' in OSM
|
---|
8 | # L3 (Chiwogs) are admin_level = '8' in OSM
|
---|
9 |
|
---|
10 | from json import loads as json_loads
|
---|
11 |
|
---|
12 | from pydriosm.reader import GeofabrikReader
|
---|
13 |
|
---|
14 | from shapely.geometry import point
|
---|
15 | from shapely.geometry import shape as geojson_loads
|
---|
16 | from shapely.wkt import dumps as wkt_dumps
|
---|
17 | from shapely.wkt import loads as wkt_loads
|
---|
18 |
|
---|
19 | geofabrik_reader = GeofabrikReader(data_dir = None)
|
---|
20 | # NB This takes a while!
|
---|
21 | pbf_raw = geofabrik_reader.read_osm_pbf("Bhutan")
|
---|
22 | multipolygons = pbf_raw["multipolygons"]["multipolygons"]
|
---|
23 |
|
---|
24 | db = current.db
|
---|
25 | table = current.s3db.gis_location
|
---|
26 | base_query = (table.level.belongs(("L1", "L2"))) & \
|
---|
27 | (table.deleted == False)
|
---|
28 |
|
---|
29 | output = "BT_L3.csv"
|
---|
30 | outputFile = open(output, "w", encoding="utf-8")
|
---|
31 | write = outputFile.write
|
---|
32 | write("Country,L1,L2,L3,WKT\n")
|
---|
33 | for feature in multipolygons:
|
---|
34 | feature = json_loads(feature)
|
---|
35 | p = feature["properties"]
|
---|
36 | if p["admin_level"] != "8":
|
---|
37 | continue
|
---|
38 | #name = '"%s"' % p["name"].split(" Chiwogs")[0].split(" Chiwog")[0]
|
---|
39 | name = p["name"].split(" Chiwogs")[0].split(" Chiwog")[0]
|
---|
40 | geom = geojson_loads(feature["geometry"])
|
---|
41 | centroid = geom.centroid
|
---|
42 | lat = centroid.y
|
---|
43 | lon = centroid.x
|
---|
44 | test = point.Point(lon, lat)
|
---|
45 | query = base_query & (table.lat_min < lat) & \
|
---|
46 | (table.lat_max > lat) & \
|
---|
47 | (table.lon_min < lon) & \
|
---|
48 | (table.lon_max > lon)
|
---|
49 | rows = db(query).select(table.name,
|
---|
50 | table.level,
|
---|
51 | table.wkt,
|
---|
52 | )
|
---|
53 | L1 = L2 = None
|
---|
54 | for row in rows:
|
---|
55 | shape = wkt_loads(row.wkt)
|
---|
56 | ok = test.intersects(shape)
|
---|
57 | if ok:
|
---|
58 | if row.level == "L1":
|
---|
59 | L1 = row.name
|
---|
60 | if L2:
|
---|
61 | break
|
---|
62 | elif row.level == "L2":
|
---|
63 | L2 = row.name
|
---|
64 | if L1:
|
---|
65 | break
|
---|
66 | WKT = '"%s"' % wkt_dumps(geom)
|
---|
67 | write("BT,%s,%s,%s,%s\n" % (L1, L2, name, WKT))
|
---|
68 |
|
---|
69 | outputFile.close()
|
---|