GIS/OpenStreetMap: chiwogs.py

File chiwogs.py, 2.3 KB (added by Fran Boon, 2 years ago)

Script to export Bhutan L3s from OSM

Line 
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
10from json import loads as json_loads
11
12from pydriosm.reader import GeofabrikReader
13
14from shapely.geometry import point
15from shapely.geometry import shape as geojson_loads
16from shapely.wkt import dumps as wkt_dumps
17from shapely.wkt import loads as wkt_loads
18
19geofabrik_reader = GeofabrikReader(data_dir = None)
20# NB This takes a while!
21pbf_raw = geofabrik_reader.read_osm_pbf("Bhutan")
22multipolygons = pbf_raw["multipolygons"]["multipolygons"]
23
24db = current.db
25table = current.s3db.gis_location
26base_query = (table.level.belongs(("L1", "L2"))) & \
27 (table.deleted == False)
28
29output = "BT_L3.csv"
30outputFile = open(output, "w", encoding="utf-8")
31write = outputFile.write
32write("Country,L1,L2,L3,WKT\n")
33for 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
69outputFile.close()