| 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()
|
|---|
