-
-
Notifications
You must be signed in to change notification settings - Fork 6
/
build_cryptids_database.py
135 lines (119 loc) · 4.68 KB
/
build_cryptids_database.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
try:
import pysqlite3 as sqlite3
except ImportError:
import sqlite3
from functools import partial, wraps
import json
import os
import pyproj
from shapely.geometry.multipolygon import MultiPolygon
from shapely.geometry import shape
from shapely.ops import transform
BASE_DIR = os.path.dirname(__file__)
try_these = (
"mod_spatialite",
"/usr/local/lib/mod_spatialite.dylib",
"/usr/lib/x86_64-linux-gnu/mod_spatialite.so",
"/usr/lib/x86_64-linux-gnu/libspatialite.so.5",
"/usr/lib/x86_64-linux-gnu/libspatialite.so.7",
)
def insert_record(conn, id, geojson):
poly = shape(geojson['geometry'])
# Buffer by 5 miles - if you are within 5 miles of a lake you should
# still hear about the lake monster!
buffered_poly = buffer_polygon_in_meters(poly, 5 * 1600)
if poly.geom_type == "Polygon":
poly = MultiPolygon([poly])
if buffered_poly.geom_type == "Polygon":
buffered_poly = MultiPolygon([buffered_poly])
conn.execute("""
INSERT INTO cryptids (
id, name, wikipedia_url, additional_url, description, copyright, first_sighted, last_sighted, geom, range)
VALUES(?, ?, ?, ?, ?, ?, ?, ?, GeomFromText(?, 4326), GeomFromText(?, 4326))
""", (
id,
geojson["properties"]["name"],
geojson["properties"].get("wikipedia_url"),
geojson["properties"].get("additional_url"),
geojson["properties"]["description"],
geojson["properties"].get("copyright"),
geojson["properties"].get("first_sighted") or "",
geojson["properties"].get("last_sighted") or "",
buffered_poly.wkt,
poly.wkt,
))
def build_database(dbpath="cryptids.db"):
conn = sqlite3.connect(dbpath)
# Enable Spatialite
conn.enable_load_extension(True)
spatialite_found = False
for extension in try_these:
try:
conn.load_extension(extension)
conn.execute('select InitSpatialMetadata(1)')
spatialite_found = True
except:
continue
assert spatialite_found, 'Did not find spatialite - looked in {}'.format(
repr(try_these)
)
# Create the database table
conn.execute('''
create table cryptids (
id integer primary key,
name text,
wikipedia_url text,
additional_url text,
description text,
copyright text,
first_sighted integer,
last_sighted integer
)
''')
conn.execute("SELECT AddGeometryColumn('cryptids', 'geom', 4326, 'MULTIPOLYGON', 2);")
conn.execute("SELECT CreateSpatialIndex('cryptids', 'geom');")
conn.execute("SELECT AddGeometryColumn('cryptids', 'range', 4326, 'MULTIPOLYGON', 2);")
conn.execute("SELECT CreateSpatialIndex('cryptids', 'range');")
# Loop through and insert the records
for filename in os.listdir(os.path.join(BASE_DIR, 'cryptids')):
filepath = os.path.join(BASE_DIR, 'cryptids', filename)
if filepath.endswith(".geojson"):
id = filename.split(".")[0]
insert_record(conn, id, json.load(open(filepath)))
conn.commit()
conn.close()
wgs84 = pyproj.Proj("+init=EPSG:4326")
def process_wgs84_in_meters(fn):
# Given a function which takes a polygon as the
# first argument, this decorator ensures that
# the polygon will first be transformed to a custom
# aeqd projection centered around the center of that
# polygon, which operates in meters. This newly
# transformed polygon will be passed to the real
# function, which can then operate on it using meters
# as the units. It can return a polygon-in-meters
# which will then be transformed back to wgs84
@wraps(fn)
def wrapper(polygon, *args, **kwargs):
# Custom aeqd projection centered on the polygon centroid
centroid = polygon.centroid
custom_projection = pyproj.Proj(
"+proj=aeqd +lat_0=%s +lon_0=%s +units=m" % (centroid.y, centroid.x)
)
project = partial(pyproj.transform, wgs84, custom_projection)
polygon_in_m = transform(project, polygon)
# Actually run the decorated function
resulting_polygon_in_m = fn(polygon_in_m, *args, **kwargs)
# Convert back to wgs84 again
reverse_project = partial(pyproj.transform, custom_projection, wgs84)
return transform(reverse_project, resulting_polygon_in_m)
return wrapper
@process_wgs84_in_meters
def buffer_polygon_in_meters(polygon, meters):
# Given a shapely polygon with its points
# defined using lat/lon, return a polygon
# representing the input expanded in all
# directions by X meters.
return polygon.buffer(meters)
if __name__ == "__main__":
build_database()