saradata.py 10.5 KB
Newer Older
1
2
3
4


import sqlite3
import csv
5
from shapely import wkt
6
7
import geopandas
from pygeotile.tile import Tile
8
9
import numpy

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

COQUIMBO_LONGITUDE = [-71.5, -71.0]
COQUIMBO_LATITUDE = [-30.25, -29.25]


def start_sara_database(database_filepath):

    # Create SQLite database for the exposure data
    print("Database filepath: " + database_filepath)
    conn = sqlite3.connect(database_filepath)
    print("Database created")
    conn.enable_load_extension(True)
    sql_statement = "SELECT load_extension('/usr/lib64/mod_spatialite.so.7');"
    conn.execute(sql_statement)

    # For debug
    rs = conn.execute('SELECT sqlite_version(), spatialite_version()')
    for row in rs:
        msg = "> SQLite v%s Spatialite v%s" % (row[0], row[1])
        print(msg)

    sql_statement = "SELECT InitSpatialMetaData(1);"
    conn.execute(sql_statement)
    print("Spatialite loaded")

    return conn


def create_sara_database(database_filepath):

    # Create SQLite database for the exposure data
    print("Database filepath: " + database_filepath)
    conn = sqlite3.connect(database_filepath)
    print("Database created")
    conn.enable_load_extension(True)
    sql_statement = "SELECT load_extension('/usr/lib64/mod_spatialite.so.7');"
    conn.execute(sql_statement)

    # For debug
    rs = conn.execute('SELECT sqlite_version(), spatialite_version()')
    for row in rs:
        msg = "> SQLite v%s Spatialite v%s" % (row[0], row[1])
        print(msg)

    sql_statement = "SELECT InitSpatialMetaData(1);"
    conn.execute(sql_statement)
    print("Spatialite loaded")

    # Create LOCATIONS table
    conn.execute('''CREATE TABLE LOCATIONS
                   (IDX   INTEGER PRIMARY KEY AUTOINCREMENT,
                    NAME  TEXT);''')
    sql_statement = "SELECT AddGeometryColumn('LOCATIONS', 'geometry', 4326, 'POINT', 'XY');"
    conn.execute(sql_statement)
64
65
    # sql_statement = "SELECT AddGeometryColumn('LOCATIONS', 'voronoi', 4326, 'POLYGON', 'XY');"
    # conn.execute(sql_statement)
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
    print("Table LOCATIONS created")


    # sql_statement = "CREATE INDEX coordinate_index ON LOCATION;"
    # conn.execute(sql_statement)
    # print("Index coordinate_index created")
    conn.execute('''CREATE TABLE ASSETS
                   (LOCATION    INTEGER,
                    TAXONOMY    TEXT,
                    NUMBER      REAL,
                    STRUCTURAL  REAL,
                    NIGHT       REAL);''')
    print("Table ASSETS created")
    sql_statement = "CREATE INDEX location_index ON ASSETS(LOCATION);"
    conn.execute(sql_statement)
    print("Index location_index created")

    # # Create Voronoi table
    # conn.execute("""CREATE TABLE Voronoi (id INTEGER NOT NULL PRIMARY KEY)""")
    # conn.execute("""SELECT AddGeometryColumn('Voronoi', 'geometry', 4326, 'POLYGON', 'XY')""")
    # conn.commit()
    # print("Table Voronoi created")

    return conn


def read_sara_exposure_into_database(conn, sara_exposure_filepath):

    with open(sara_exposure_filepath) as csv_file:
        csv_reader = csv.reader(csv_file, delimiter=',')
        next(csv_reader)  # Skip first line
        last_longitude = 1000
        last_latitude = 1000
        location_count = 0
100
        asset_count = 0
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
        for row in csv_reader:
            # Create location in database
            if (float(row[0]) < COQUIMBO_LONGITUDE[0]) or (float(row[0]) > COQUIMBO_LONGITUDE[1]) or \
               (float(row[1]) < COQUIMBO_LATITUDE[0]) or (float(row[1]) > COQUIMBO_LATITUDE[1]):
                continue
            if not ((last_longitude == row[0]) and (last_latitude == row[1])):
                # Create new location entry
                location_count += 1
                location = "GeomFromText('POINT(%f %f)', 4326)" % (float(row[0]), float(row[1]))
                sql_statement = "INSERT INTO LOCATIONS (geometry, NAME) VALUES (%s, '%s')" % (location, row[3])
                print("SQL statement: " + sql_statement)
                conn.execute(sql_statement)
                conn.commit()
                print("Location created")
                last_longitude = row[0]
                last_latitude = row[1]
            else:
118
119
120
121
122
123
                sql_statement = 'INSERT INTO ASSETS (LOCATION, TAXONOMY, NUMBER, STRUCTURAL, NIGHT) VALUES (' \
                                + str(location_count) + ', """' + row[5] + '""", ' + row[6] + ', ' + row[7] + ', ' + row[8] + ')'
                print("SQL statement: " + sql_statement)
                conn.execute(sql_statement)
                conn.commit()
                print("Asset created")
124
125
                asset_count += 1
    print("Number of assets: " + str(asset_count))
126

Danijel Schorlemmer's avatar
Danijel Schorlemmer committed
127

128
def add_voronoi_cells(conn):
129
130
131

    cur = conn.cursor()

132
    # Create the Voronoi table
133
134
    sql = 'CREATE TABLE Voronoi ('
    sql += 'loc_id INTEGER)'
135
136
    cur.execute(sql)
    # creating a POLYGON Geometry column
137
    sql = "SELECT AddGeometryColumn('Voronoi', "
138
139
    sql += "'geom', 4326, 'POLYGON', 'XY')"
    cur.execute(sql)
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
    print("Table Voronoi created")

    cur.execute("""SELECT ST_AsText(ST_VoronojDiagram(ST_Collect(geometry))) FROM LOCATIONS""")
    conn.commit()
    diag = cur.fetchall()
    multipolygon = diag[0][0]
    print(multipolygon)
    p1 = wkt.loads(multipolygon)
    for polygon in p1:  # same for multipolygon.geoms
        print(polygon)
        location = "GeomFromText('%s', 4326)" % polygon
        sql_statement = "INSERT INTO Voronoi (geom) VALUES (%s)" % location
        print("SQL statement: " + sql_statement)
        conn.execute(sql_statement)
        conn.commit()

    # Loop over points in Locations
    sql_statement = "SELECT idx, ST_AsText(geometry) FROM LOCATIONS;"
    print(sql_statement)
    cur.execute(sql_statement)
    rows = cur.fetchall()
    for row in rows:
        print("  ")
        print("Point")
        print(row)
        location_index = row[0]
        print("Location index: " + str(location_index))
        print("Point coordinates: " + str(row[1]))
        sql_statement = "SELECT ST_AsText(geom) as voro FROM voronoi where within(GeomFromText('%s', 4326), voronoi.geom);" % row[1]
        print(sql_statement)
        cur.execute(sql_statement)
        polys = cur.fetchall()
        print(polys[0][0])
        # Add Voronoi polygon to LOCATIONS table
        sql_statement = "UPDATE Voronoi SET loc_id = %d WHERE geom = GeomFromText('%s', 4326);" % (location_index, polys[0][0])
        print(sql_statement)
        cur.execute(sql_statement)
177

Danijel Schorlemmer's avatar
Danijel Schorlemmer committed
178
    conn.commit()
179

180

181
def add_buildings(conn, building_filepath):
182

183
184
185
    cur = conn.cursor()
    # Create the Buildings table
    sql = 'CREATE TABLE Buildings ('
186
    sql += 'IDX   INTEGER PRIMARY KEY AUTOINCREMENT,'
187
188
    sql += 'loc_id INTEGER,'
    sql += 'quadkey TEXT);'
189
190
    cur.execute(sql)
    # creating a POLYGON Geometry column
191
192
    sql = "SELECT AddGeometryColumn('Buildings', "
    sql += "'geom', 4326, 'MULTIPOLYGON', 'XY')"
193
    cur.execute(sql)
194
    print("Table Buildings created")
195

196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
    buildings = geopandas.read_file(building_filepath)
    for index, row in buildings.iterrows():  # Looping over all points
        sql_statement = "SELECT loc_id, geom FROM voronoi where within(GeomFromText('%s', 4326), voronoi.geom);" % row.geometry.centroid
        # print(sql_statement)
        cur.execute(sql_statement)
        polys = cur.fetchall()
        try:
            location_index = polys[0][0]
        except:  # If centroid lies on the voronoi boundary
            centroid = "POINT(%f %f)" % (row.geometry.centroid.x + 0.000001, row.geometry.centroid.y + 0.000001)
            print(centroid)
            sql_statement = "SELECT loc_id, geom FROM voronoi where within(GeomFromText('%s', 4326), voronoi.geom);" % centroid
            cur.execute(sql_statement)
            polys = cur.fetchall()
            location_index = polys[0][0]

        print(location_index)

        quadkey = Tile.for_latitude_longitude(longitude=float(row.geometry.centroid.x), latitude=float(row.geometry.centroid.y), zoom=18).quad_tree
        building_geometry = "GeomFromText('%s', 4326)" % (row.geometry)
        sql_statement = "INSERT INTO Buildings (loc_id, quadkey, geom) VALUES (%d, %s, %s)" % (location_index, quadkey, building_geometry)
        print("SQL statement: " + sql_statement)
        conn.execute(sql_statement)
        conn.commit()
220
221


222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
def classify_buildings(conn):

    # Create building classification table
    cur = conn.cursor()
    # Create the Buildings table
    # sql = 'CREATE TABLE Exposure ('
    # sql += 'building_id INTEGER, '
    # sql += 'taxonomy TEXT, '
    # sql += 'NUMBER   REAL, '
    # sql += 'STRUCTURAL  REAL, '
    # sql += 'NIGHT       REAL);'
    # cur.execute(sql)
    # print("Table Exposure created")

    # Simple version
    # Loop over all exposure-data locations (Voronoi cells)
    sql_statement = "SELECT idx FROM Locations;"
    print(sql_statement)
    cur.execute(sql_statement)
    rows = cur.fetchall()
    for row in rows:
        print("Location: " + str(row[0]))
        sql_statement = "SELECT Taxonomy, Number, Structural, Night FROM Assets WHERE Location = %d" % row[0]
        cur.execute(sql_statement)
        assets = cur.fetchall()
        for asset in assets:
            print("Asset: %s, %f, %f, %f" % (asset[0], asset[1], asset[2], asset[3]))
        taxonomies = numpy.array(assets)[:, 0]
        print(taxonomies)
        a = numpy.asarray(assets)[:, 1:].astype(numpy.float)
        print(a)
        sum = a.sum(axis=0)
        print(sum)
        proportions = a/sum[None, :]
        print(proportions)
        proportions[numpy.isnan(proportions)] = 0
        print(proportions)
        # Select all buildings in the Voronoi cell
        sql_statement = "SELECT idx FROM Buildings WHERE loc_id = %d" % row[0]
        print(sql_statement)
        cur.execute(sql_statement)
        buildings = cur.fetchall()
        for building_index in buildings:
            # Write exposure data for the building
            print(building_index[0])
            for counter in range(taxonomies.shape[0]):
                print(taxonomies[counter])
                print(proportions[counter,:])

                sql_statement = "INSERT INTO Exposure (building_id, taxonomy, number, structural, night) VALUES (%d, %s, %f, %f, %f)" % (building_index[0], taxonomies[counter], proportions[counter,0], proportions[counter,0], proportions[counter,0])
                cur.execute(sql_statement)
        conn.commit()




278
def main():
279
280
281
282
283
284
285
    # conn = create_sara_database('../skr/testdb.sqlite')
    conn = start_sara_database('../skr/testdb.sqlite')
    # read_sara_exposure_into_database(conn, '../data/Exposure_Ind_Chile.csv')
    # add_voronoi_cells(conn)
    # add_buildings(conn, '../data/obm.building.small.gpkg')
    classify_buildings(conn)

286

287
288

main()