davegimo

Untitled

Feb 21st, 2021
40
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import magnet as rip
  2. import psycopg2
  3. import googlemaps
  4. from geomet import wkt
  5. import json
  6.  
  7. gmaps = googlemaps.Client(key='AIzaSyAthjChcZHYkj8y7j5_drctNkvwBktzfj4')
  8.  
  9.  
  10. n_clusters_civici = 4
  11. n_clusters = 50 # 20
  12. iterazioni = 300
  13. init = 10
  14.  
  15.  
  16.  
  17. con = psycopg2.connect(database="carto", user="pcu_master", password="pcu_master", host="rdsppaldc.ci6trjgque6a.eu-central-1.rds.amazonaws.com", port="5432")
  18. print("Database opened successfully")
  19.  
  20. cur = con.cursor()
  21. cur.execute(''' select streetnum_id, prefix, street_name, address_num, comune, geometry, dx, dy
  22. from enet.t_civici_dave
  23. where cft = 'D61D159'
  24. ''')
  25.  
  26. rows = cur.fetchall()
  27.  
  28. dict = {}
  29. dict["type"] = "FeatureCollection"
  30. dict["features"] = []
  31. lista = []
  32.  
  33. for row in rows:
  34. feature = {}
  35.  
  36. x = wkt.loads(row[5])
  37. s = "(" + str(x["coordinates"][0]) + ", " + str(x["coordinates"][1]) + ")"
  38.  
  39. feature["civico_id"] = int(row[0])
  40. feature["pref_via"] = row[1]
  41. feature["nome_via"] = row[2]
  42. feature["num_civico"] = row[3]
  43. feature["comune"] = row[4]
  44. feature["dx"] = row[6]
  45. feature["dy"] = row[7]
  46. lista.append([row[6], row[7]])
  47. feature["coord_wgs84"] = s
  48.  
  49.  
  50. dict["features"].append(feature)
  51.  
  52.  
  53. print("Operation done successfully")
  54. con.close()
  55.  
  56.  
  57. city = row[4]
  58. province = "RC"
  59. aggiornati = 0 #se sono molti gli aggiornati, allora mando in postgres
  60.  
  61. for p in dict["features"]:
  62. print(p)
  63. if p["dx"] != 0:
  64. print("gia fatto")
  65. continue
  66. s = rip.pref_via(p['pref_via'])
  67. if s == "":
  68. continue
  69. else:
  70. aggiornati += 1
  71. via = s + p['nome_via'] + " " + str(p['num_civico']) + ", " + city + " " + province
  72. wgs84 = rip.parse_tuple(p['coord_wgs84'])
  73. lat = wgs84[1]
  74. lng = wgs84[0]
  75. geocode_result = gmaps.geocode(via)
  76. if geocode_result == []:
  77. continue
  78. delta_x = geocode_result[0]["geometry"]["location"]['lng'] - lng
  79. delta_y = geocode_result[0]["geometry"]["location"]['lat'] - lat
  80.  
  81. print(delta_x, delta_y)
  82.  
  83. p["dx"] = delta_x
  84. p["dy"] = delta_y
  85.  
  86.  
  87.  
  88. with open('civici.json', 'w') as fp:
  89. json.dump(dict, fp, indent=4)
  90. print("Civici salvati")
  91.  
  92.  
  93.  
  94.  
  95. if aggiornati > 300:
  96. lista = []
  97. print("Aggiornamento . . .", aggiornati)
  98. con = psycopg2.connect(database="carto", user="pcu_master", password="pcu_master",
  99. host="rdsppaldc.ci6trjgque6a.eu-central-1.rds.amazonaws.com", port="5432")
  100. cur = con.cursor()
  101. sql_update_query = """ Update enet.t_civici_dave set dx = %s, dy = %s where streetnum_id = %s """
  102.  
  103. for p in dict["features"]:
  104. cur.execute(sql_update_query, (p["dx"], p["dy"], p["civico_id"]))
  105. con.commit()
  106. count = cur.rowcount
  107. print(count, "Record Updated successfully")
  108. lista.append([p["dx"], p["dy"]])
  109. con.close()
  110.  
  111. else:
  112. print("inutile")
  113.  
  114.  
  115.  
  116. cluster = rip.find_delta_k_means(lista, n_clusters, iterazioni, init)
  117.  
  118.  
  119.  
  120. print("Aggiornamento rete. . .")
  121. con = psycopg2.connect(database="carto", user="pcu_master", password="pcu_master",
  122. host="rdsppaldc.ci6trjgque6a.eu-central-1.rds.amazonaws.com", port="5432")
  123. cur = con.cursor()
  124. sql_update_query = """ Update enet.t_ramo_bt_dave set newgeom = ST_Translate(geom, %s, %s) where cft = 'D61D159' """
  125. cur.execute(sql_update_query, (cluster[0], cluster[1]))
  126. con.commit()
  127. count = cur.rowcount
  128. print(count, "Record Updated successfully")
  129. con.close()
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136.  
  137.  
  138.  
  139.  
  140.  
  141.  
  142. '''
  143.  
  144.  
  145. city = "Reggio Calabria"
  146. province = "RC"
  147.  
  148. civici = "coords/civici/" + city + ".json"
  149. file_in = "coords/AGUI/" + city + ".json"
  150. file_out = "coords/NEW/" + city + ".json"
  151. file_cluster = "coords/cluster/" + city + ".json"
  152.  
  153. html_agui = "build/" + city + "/old.html"
  154. html_new = "build/" + city + "/new.html"
  155. html_civici = "build/" + city + "/civici.html"
  156. html_cluster = "build/" + city + "/civici_cluster.html"
  157. html_linee = "build/" + city + "/linee_cluster.html"
  158. html_new_cluster = "build/" + city + "/new_cluster.html"
  159.  
  160. delta = "delta/" + city
  161.  
  162. n_clusters_civici = 4
  163. n_clusters = 50 # 20
  164. iterazioni = 300
  165. init = 10
  166.  
  167. colors = []
  168. for i in range(n_clusters):
  169. colors.append(rip.random_color())
  170.  
  171. # clustering dei civici e della linea
  172. print("Clustering civici")
  173. lista_civici = rip.plot_civici(civici, html_civici, city, province)
  174. clustered_list, centroids, cluster_assignments = rip.civici_k_means(lista_civici, n_clusters_civici, iterazioni, init)
  175. rip.plot_cluster_civici(clustered_list, html_cluster, city, province, n_clusters_civici, colors)
  176. clustered_lines, line_assignment = rip.cluster_assignment(file_in, centroids, n_clusters_civici)
  177. rip.plot_cluster_linee(clustered_lines, html_linee, city, province, n_clusters_civici, colors)
  178.  
  179. rip.plot_on_map(file_in, html_agui, city, province)
  180.  
  181. if rip.exists(delta):
  182. print("Lista dei delta già presente")
  183. lista = rip.load_list(delta)
  184.  
  185.  
  186. else:
  187. print("Lista dei delta da calcolare")
  188. # rip.plot_on_map(file_in, html_agui, city, province)
  189. lista = rip.create_delta_list(civici, city, province)
  190. print("Lista creata")
  191. rip.save_list(lista, delta)
  192.  
  193. print("elementi da clusterizzare: " + str(len(lista)))
  194. cluster = rip.find_delta_k_means(lista, n_clusters, iterazioni, init)
  195. rip.apply_delta(file_in, file_out, cluster[0], cluster[1])
  196. rip.plot_on_map(file_out, html_new, city, province)
  197.  
  198. print("CLUSTER = " + str(cluster))
  199.  
  200. ################# NUOVA PARTE
  201. cluster_delta = []
  202. for i in range(n_clusters_civici):
  203. cluster_delta.append([])
  204.  
  205. size_for_range = len(cluster_assignments)
  206. if size_for_range > len(lista):
  207. size_for_range = len(lista)
  208.  
  209. if len(lista) != len(cluster_assignments):
  210. print("DIVERSI!")
  211. else:
  212. print("UGUALI!")
  213.  
  214. # smistaggio da rivedere
  215. for i in range(size_for_range):
  216. cluster_delta[cluster_assignments[i]].append(lista[i])
  217.  
  218. deltas = []
  219. for i in range(len(cluster_delta)):
  220. c = rip.best_number(len(cluster_delta[i]))
  221. cluster = rip.find_delta_k_means(cluster_delta[i], c, iterazioni, init)
  222. deltas.append(cluster)
  223.  
  224. rip.apply_deltas(file_in, file_cluster, deltas, line_assignment)
  225. rip.plot_on_map(file_cluster, html_new_cluster, city, province, "black")
  226. # rip.plot_cluster_linee(clustered_lines,html_new_cluster,city, province, n_clusters_civici, colors, flag=False)
  227.  
  228.  
  229.  
  230. line assignment mi serve per la lista dei delta da suddividere dai delta calcolati
  231.  
  232. creo un elenco di n_cluster delta, per ognuno di essi calcolo il delta ottimo con kmeans
  233. lo applico su un singolo cluster
  234. nuova funzione apply delta per ritornare il json
  235.  
  236.  
  237.  
  238. '''
  239.  
  240.  
  241. ################
  242.  
  243.  
  244.  
  245. import googlemaps
  246. import json
  247. import ast
  248. import numpy as np
  249. from matplotlib import pyplot as plt
  250. from sklearn.cluster import KMeans
  251. from collections import Counter
  252. import folium
  253. import pickle
  254.  
  255.  
  256. gmaps = googlemaps.Client(key='AIzaSyAthjChcZHYkj8y7j5_drctNkvwBktzfj4')
  257.  
  258.  
  259. def getCoords(city, province):
  260. geocode_result = gmaps.geocode(city + ", " + province + " Italia")
  261. return (geocode_result[0]["geometry"]["location"]['lat'], geocode_result[0]["geometry"]["location"]['lng'])
  262.  
  263. def inserisciMarker(mappa, lat, lng, popup, icon, color):
  264. folium.Marker(location=[lat, lng], tooltip=popup, icon=folium.Icon(icon=icon, color=color)).add_to(mappa)
  265.  
  266.  
  267. def inserisciLinea(mappa, points, color, w, opacity, tooltip):
  268. folium.PolyLine(points, tooltip=tooltip, color=color, weight=w, opacity=opacity).add_to(mappa)
  269.  
  270.  
  271. def inserisciCerchio(mappa, lat, lng, radius):
  272. folium.CircleMarker([lat, lng], radius=radius * 10 / 2, fill=True, opacity=0.1).add_to(mappa)
  273.  
  274.  
  275. def inserisciPunto(mappa, lat, lng, txt):
  276. folium.CircleMarker([lat, lng], radius=5, fill=True, color="black",opacity=0.2, popup=txt).add_to(mappa)
  277.  
  278. def inserisciPuntoCustom(mappa, lat, lng, txt, color):
  279. folium.CircleMarker([lat, lng], radius=5, fill=True, color=color,opacity=0.8, popup=txt).add_to(mappa)
  280.  
  281. def inserisciPuntoR(mappa, lat, lng, txt):
  282. folium.CircleMarker([lat, lng], radius=5, fill=True, color="red", popup=txt).add_to(mappa)
  283.  
  284.  
  285.  
  286. def random_color():
  287. r = np.random.randint(255)
  288. g = np.random.randint(255)
  289. b = np.random.randint(0,100)
  290. rgb = (r,g,b)
  291. return '#%02x%02x%02x' % rgb
  292.  
  293.  
  294. def disegnaLinee(lista,m,color):
  295. for i in range(len(lista) - 1):
  296. l = []
  297. l.append((lista[i][1], lista[i][0]))
  298. l.append((lista[i + 1][1], lista[i + 1][0]))
  299. inserisciLinea(m, l, color, 2, 0.8, "")
  300.  
  301. def parse_tuple(string):
  302. try:
  303. s = ast.literal_eval(str(string))
  304. if type(s) == tuple:
  305. return s
  306. return
  307. except:
  308. return
  309.  
  310. def pref_via(s):
  311. if "STR" in s:
  312. return ""
  313. if "PZA" in s:
  314. return "Piazza"
  315. if "CON" in s:
  316. return "Contrada"
  317. return "Via"
  318.  
  319.  
  320.  
  321. def create_delta_list(file, city, province):
  322.  
  323. with open(file) as json_file:
  324. json_data = json.load(json_file)
  325.  
  326. lista = []
  327.  
  328. for p in json_data['features']:
  329. s = pref_via(p['pref_via'])
  330. if s == "":
  331. continue
  332. else:
  333. via = s + p['nome_via'] + " " + str(p['num_civico']) + ", " + city + " " + province
  334. wgs84 = parse_tuple(p['coord_wgs84'])
  335. lat = wgs84[1]
  336. lng = wgs84[0]
  337. geocode_result = gmaps.geocode(via)
  338. if geocode_result == []:
  339. continue
  340. delta_x = geocode_result[0]["geometry"]["location"]['lng'] - lng
  341. delta_y = geocode_result[0]["geometry"]["location"]['lat'] - lat
  342. lista.append([delta_x, delta_y])
  343. #print(delta_x, delta_y)
  344.  
  345. return lista
  346.  
  347. def save_list(lista, outfile):
  348. with open(outfile, 'wb') as fp:
  349. pickle.dump(lista, fp)
  350.  
  351. def load_list(file):
  352. with open(file, 'rb') as fp:
  353. itemlist = pickle.load(fp)
  354. return itemlist
  355.  
  356. def exists(file):
  357. try:
  358. f = open(file)
  359. flag = True
  360. except IOError:
  361. flag = False
  362. finally:
  363. return flag
  364.  
  365. def find_delta_k_means(lista, n_clusters, iterazioni, init):
  366. X = np.array(lista)
  367. #print(X)
  368.  
  369. kmeans = KMeans(n_clusters=n_clusters, init='k-means++', max_iter=iterazioni, n_init=init, random_state=0)
  370. pred_y = kmeans.fit_predict(X)
  371.  
  372. # plt.scatter(X[:, 0], X[:, 1])
  373. # plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], s=300, c='red')
  374. # plt.show()
  375.  
  376. #print(pred_y)
  377. l = list(pred_y)
  378.  
  379. data = Counter(l)
  380. #print(data.most_common(1))
  381. most_common = data.most_common(1)[0]
  382. cluster = most_common[0]
  383. lunghezza = most_common[1]
  384. #print(cluster)
  385. print(lunghezza)
  386.  
  387. print(kmeans.cluster_centers_[cluster])
  388.  
  389. return kmeans.cluster_centers_[cluster]
  390.  
  391. def civici_k_means(lista, n_clusters, iterazioni, init):
  392. X = np.array(lista)
  393. kmeans = KMeans(n_clusters=n_clusters, init='k-means++', max_iter=iterazioni, n_init=init, random_state=0)
  394. pred_y = kmeans.fit_predict(X)
  395.  
  396. # plt.scatter(X[:, 0], X[:, 1])
  397. # plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], s=300, c='red')
  398. # plt.show()
  399.  
  400. centroids = kmeans.cluster_centers_
  401. print(centroids)
  402.  
  403. size = len(lista)
  404. clustered_list = []
  405.  
  406. for e in range(n_clusters):
  407. clustered_list.append([])
  408.  
  409. cluster_assignment = []
  410. for e in range(size):
  411. clustered_list[pred_y[e]].append(lista[e])
  412. cluster_assignment.append(pred_y[e])
  413.  
  414. return clustered_list, centroids, cluster_assignment
  415.  
  416. def best_number(n):
  417. print("Cluster = " + str(n))
  418. if n < 20:
  419. return 6
  420. if n < 50:
  421. return 8
  422. if n < 180:
  423. return 10
  424. if n < 300:
  425. return 15
  426. if n < 500:
  427. return 25
  428. return 50
  429.  
  430. def kmeans_assignment(centroids, points):
  431. points = np.array(points)
  432.  
  433. num_centroids, dim = centroids.shape
  434. num_points, _ = points.shape
  435.  
  436. # Tile and reshape both arrays into `[num_points, num_centroids, dim]`.
  437. centroids = np.tile(centroids, [num_points, 1]).reshape([num_points, num_centroids, dim])
  438. points = np.tile(points, [1, num_centroids]).reshape([num_points, num_centroids, dim])
  439.  
  440. # Compute all distances (for all points and all centroids) at once and
  441. # select the min centroid for each point.
  442. distances = np.sum(np.square(centroids - points), axis=2)
  443. return np.argmin(distances, axis=1)
  444.  
  445. def cluster_assignment(file, centroids, n_clusters):
  446.  
  447. with open(file) as json_file:
  448. json_data = json.load(json_file)
  449. #print(json_data)
  450.  
  451. clustered_list = []
  452. line_assignment = []
  453. for e in range(n_clusters):
  454. clustered_list.append([])
  455.  
  456. for p in json_data['features']:
  457. if p["geometry"]["type"] == "LineString":
  458. x = kmeans_assignment(centroids,p["geometry"]["coordinates"])
  459. clustered_list[np.bincount(x).argmax()].append(p["geometry"]["coordinates"])
  460. line_assignment.append(np.bincount(x).argmax())
  461.  
  462. return clustered_list,line_assignment
  463.  
  464.  
  465.  
  466. def apply_delta(file_in, file_out, delta_x, delta_y):
  467. with open(file_in) as json_file:
  468. json_data = json.load(json_file)
  469. #print(json_data)
  470. for p in json_data['features']:
  471. if p["geometry"]["type"] == "Point":
  472. p["geometry"]["coordinates"][0] += delta_x
  473. p["geometry"]["coordinates"][1] += delta_y
  474. elif p["geometry"]["type"] == "LineString":
  475. for c in p["geometry"]["coordinates"]:
  476. c[0] += delta_x
  477. c[1] += delta_y
  478.  
  479. with open(file_out, 'w') as outfile:
  480. json.dump(json_data, outfile, sort_keys=True, indent=4,
  481. ensure_ascii=False)
  482.  
  483. def apply_deltas(file_in, file_out, deltas, line_assignment):
  484. with open(file_in) as json_file:
  485. json_data = json.load(json_file)
  486. i = 0
  487. for p in json_data['features']:
  488.  
  489. delta_x = deltas[line_assignment[i]][0]
  490. delta_y = deltas[line_assignment[i]][1]
  491.  
  492. if p["geometry"]["type"] == "Point":
  493. p["geometry"]["coordinates"][0] += delta_x
  494. p["geometry"]["coordinates"][1] += delta_y
  495. elif p["geometry"]["type"] == "LineString":
  496. for c in p["geometry"]["coordinates"]:
  497. c[0] += delta_x
  498. c[1] += delta_y
  499. i += 1
  500.  
  501. with open(file_out, 'w') as outfile:
  502. json.dump(json_data, outfile, sort_keys=True, indent=4,
  503. ensure_ascii=False)
  504.  
  505. def plot_on_map(file,file_out,city, province,color="blue"):
  506.  
  507. coords = getCoords(city, province)
  508. m = folium.Map(location=[coords[0], coords[1]], tiles='OpenStreetMap', zoom_start=16)
  509.  
  510. with open(file) as json_file:
  511. json_data = json.load(json_file)
  512. #print(json_data)
  513.  
  514. for p in json_data['features']:
  515. if p["geometry"]["type"] == "Point":
  516. inserisciPunto(m, p["geometry"]["coordinates"][1], p["geometry"]["coordinates"][0], "test")
  517. elif p["geometry"]["type"] == "LineString":
  518. disegnaLinee(p["geometry"]["coordinates"], m,color)
  519.  
  520. m.save(file_out)
  521.  
  522.  
  523. def plot_civici(file,file_out,city, province):
  524.  
  525. coords = getCoords(city, province)
  526. m = folium.Map(location=[coords[0], coords[1]], tiles='OpenStreetMap', zoom_start=16)
  527. l = []
  528. with open(file) as json_file:
  529. json_data = json.load(json_file)
  530.  
  531. for p in json_data['features']:
  532. punto = parse_tuple(p["coord_wgs84"])
  533. inserisciPunto(m, punto[1], punto[0], "test")
  534. l.append(punto)
  535.  
  536. m.save(file_out)
  537.  
  538.  
  539. return l
  540.  
  541.  
  542. def plot_cluster_civici(lista,file_out,city, province, n, colors):
  543.  
  544. coords = getCoords(city, province)
  545. m = folium.Map(location=[coords[0], coords[1]], tiles='OpenStreetMap', zoom_start=16)
  546.  
  547.  
  548. for i in range(n):
  549. for e in lista[i]:
  550. inserisciPuntoCustom(m,e[1],e[0], "cluster " + str(i), colors[i])
  551.  
  552. m.save(file_out)
  553.  
  554. def swap(l):
  555. for e in l:
  556. temp = e[0]
  557. e[0] = e[1]
  558. e[1] = temp
  559. return l
  560.  
  561. def plot_cluster_linee(lista,file_out,city, province, n, colors, flag=True):
  562.  
  563. coords = getCoords(city, province)
  564. m = folium.Map(location=[coords[0], coords[1]], tiles='OpenStreetMap', zoom_start=16)
  565.  
  566. for i in range(n):
  567. for e in lista[i]:
  568. if flag:
  569. e = swap(e)
  570. inserisciLinea(m, e, colors[i], 2, 0.7, "cluster " + str(i))
  571.  
  572. m.save(file_out)
RAW Paste Data