%matplotlib inline
import requests
import pandas as pd
import numpy as np
import json
import time
r = requests.get('http://api.gios.gov.pl/pjp-api/rest/station/findAll')
allStations = pd.io.json.json_normalize(r.json())
pollutants = ["PM2.5", "PM10", "C6H6", "NO2", "SO2", "O3", "CO"]
def get_sensors(Id):
stationId = Id
r = requests.get('http://api.gios.gov.pl/pjp-api/rest/station/sensors/' + str(stationId))
sensors = pd.io.json.json_normalize(r.json())
return sensors[["param.paramCode","id"]].set_index("param.paramCode").to_dict(orient = "index")
allStations["gegrLat"] = allStations["gegrLat"].astype("float")
allStations["gegrLon"] = allStations["gegrLon"].astype("float")
%%time
allStations["sensors"] = allStations["id"].apply(get_sensors)
CPU times: user 3.13 s, sys: 91 ms, total: 3.22 s Wall time: 6min 43s
allStations2 = pd.concat([allStations, allStations['sensors'].str.join('|').str.get_dummies()], axis=1, sort=False)
allStations2[pollutants].sum()
PM2.5 64 PM10 136 C6H6 54 NO2 144 SO2 126 O3 101 CO 77 dtype: int64
#points = allStations2[allStations2["PM2.5"] == 1][["gegrLat","gegrLon"]].values
#points = allStations2[allStations2["PM10"] == 1][["gegrLat","gegrLon"]].values
#points = allStations2[allStations2["C6H6"] == 1][["gegrLat","gegrLon"]].values
#points = allStations2[allStations2["NO2"] == 1][["gegrLat","gegrLon"]].values
#points = allStations2[allStations2["SO2"] == 1][["gegrLat","gegrLon"]].values
#points = allStations2[allStations2["O3"] == 1][["gegrLat","gegrLon"]].values
points = allStations2[allStations2["CO"] == 1][["gegrLat","gegrLon"]].values
points
array([[51.129378, 17.02925 ], [51.086225, 17.012689], [51.204503, 16.180513], [50.768729, 16.269677], [51.150391, 15.008175], [50.913433, 15.765608], [51.119011, 15.275539], [53.121764, 17.987906], [53.134083, 17.995708], [51.75805 , 19.529786], [51.775411, 19.4509 ], [51.856692, 19.421231], [51.754613, 19.434925], [51.404406, 19.696956], [51.067439, 19.448694], [52.080625, 21.111186], [50.349608, 18.236575], [50.024242, 22.010575], [50.040675, 22.004656], [50.529892, 22.112467], [54.305908, 22.307681], [53.694628, 19.968892], [52.398175, 16.959519], [52.420319, 16.877289], [51.749053, 18.048389], [52.225633, 18.269036], [53.154408, 16.759572], [52.449331, 16.999683], [50.329111, 19.231222], [50.246795, 19.019469], [50.3165 , 18.772375], [50.111181, 18.516139], [50.029416, 18.689527], [49.802075, 19.04861 ], [50.817676, 19.117426], [50.836389, 19.130111], [49.738136, 18.639069], [50.007629, 18.455548], [50.878998, 20.633692], [50.429014, 21.277367], [51.1211 , 20.880631], [53.789233, 20.486075], [54.167847, 19.410942], [50.057678, 19.926189], [50.069308, 20.053492], [50.018253, 20.992578], [53.126689, 23.155869], [53.859528, 23.00075 ], [54.353336, 18.635283], [54.400833, 18.657497], [54.560836, 18.493331], [54.328336, 18.557781], [54.431667, 18.579722], [54.380279, 18.620274], [54.463611, 17.046722], [54.120694, 17.975861], [54.031247, 19.032899], [54.546167, 17.746194], [54.104111, 18.182972], [50.159406, 19.477464], [49.293564, 19.960083], [52.219298, 21.004724], [52.556279, 19.687672], [52.550938, 19.709791], [51.399084, 21.147474], [51.83512 , 20.791556], [52.115725, 21.237297], [52.656866, 18.987368], [51.259431, 22.569133], [52.738214, 15.228667], [51.939783, 15.518861], [52.437722, 15.122444], [51.799722, 16.3175 ], [51.642656, 15.127808], [53.017628, 18.612808], [52.658467, 19.059314], [53.49355 , 18.762139]])
lat_max = points[:,0].max()
lat_min = points[:,0].min()
lon_max = points[:,1].max()
lon_min = points[:,1].min()
lat_max, lat_min, lon_max, lon_min
(54.560836, 49.293564, 23.155869, 15.008175)
from scipy.spatial import voronoi_plot_2d
from scipy.spatial import Voronoi
vor = Voronoi(points)
vertices = pd.DataFrame(vor.vertices, columns = ["Lat", "Lon"])
import matplotlib.pyplot as plt
voronoi_plot_2d(vor)
plt.show()
vertices[(vertices["Lat"] < lat_max) & (vertices["Lat"] > lat_min)]
Lat | Lon | |
---|---|---|
2 | 49.897723 | 18.583786 |
3 | 49.854969 | 18.830590 |
4 | 49.308101 | 16.986539 |
5 | 51.003759 | 17.803593 |
6 | 54.479894 | 18.121208 |
7 | 54.375205 | 18.275268 |
8 | 53.296104 | 22.147586 |
9 | 54.075733 | 15.687073 |
10 | 52.667340 | 16.069961 |
11 | 52.708752 | 18.240208 |
12 | 52.945606 | 18.279678 |
13 | 53.513054 | 17.383704 |
14 | 50.008111 | 18.573507 |
15 | 50.193205 | 18.660792 |
16 | 50.159064 | 18.315523 |
17 | 50.384333 | 18.507644 |
18 | 51.782235 | 17.330251 |
19 | 49.849356 | 17.101874 |
20 | 50.846444 | 17.701989 |
21 | 54.413209 | 18.412363 |
24 | 52.597859 | 16.056563 |
25 | 52.582287 | 16.043141 |
26 | 51.935533 | 16.790836 |
27 | 51.761505 | 16.953055 |
30 | 52.429693 | 15.623297 |
31 | 53.628493 | 18.040671 |
32 | 54.107948 | 17.443307 |
33 | 54.466625 | 18.107778 |
34 | 53.721794 | 17.300789 |
35 | 53.615311 | 17.385385 |
... | ... | ... |
108 | 50.137116 | 18.855150 |
109 | 49.633657 | 19.552235 |
110 | 50.036401 | 19.216659 |
111 | 50.047572 | 19.218790 |
112 | 52.164610 | 19.034313 |
113 | 53.257792 | 19.421552 |
114 | 53.134216 | 19.074293 |
115 | 53.059749 | 19.014394 |
117 | 49.850550 | 21.502581 |
118 | 49.970607 | 21.499921 |
119 | 51.417238 | 20.348145 |
120 | 51.421510 | 20.322266 |
121 | 50.737898 | 21.014177 |
122 | 50.447929 | 20.811464 |
123 | 51.401440 | 18.963387 |
124 | 51.408946 | 19.337980 |
125 | 50.648311 | 19.509588 |
126 | 50.618141 | 19.612016 |
127 | 50.899488 | 20.029468 |
128 | 52.185398 | 19.609821 |
129 | 52.294766 | 19.322655 |
130 | 53.123465 | 19.836333 |
131 | 53.213456 | 19.472058 |
132 | 50.719729 | 20.000882 |
133 | 50.334253 | 20.538827 |
134 | 49.622034 | 20.500106 |
135 | 50.624694 | 19.818841 |
136 | 50.678968 | 19.933613 |
137 | 49.660411 | 19.600233 |
138 | 49.679250 | 20.024944 |
130 rows × 2 columns
proper_vertices = pd.DataFrame(vertices[(vertices["Lon"] < lon_max) &
(vertices["Lon"] > lon_min) &
(vertices["Lat"] < lat_max) &
(vertices["Lat"] > lat_min)], copy = True)
import geopy.distance
for index, row in proper_vertices.iterrows():
list_of_dists = []
for point in points:
#list_of_dists.append(np.sqrt(np.square(row["Lat"] - point[0]) + np.square(row["Lon"] - point[0])))
coords_1 = (row["Lat"], row["Lon"])
coords_2 = (point[0], point[1])
#print(coords_1,coords_2)
list_of_dists.append(geopy.distance.vincenty(coords_1, coords_2).km)
#print(list_of_dists)
proper_vertices.at[index, 'dist'] = min(list_of_dists)
proper_vertices.sort_values("dist").tail(n=1)
Lat | Lon | dist | |
---|---|---|---|
4 | 49.308101 | 16.986539 | 128.84905 |
import folium
from folium.plugins import MarkerCluster
boulder_coords = [lat_min + (lat_max - lat_min)/2, lon_min + (lon_max - lon_min)/2]
my_map = folium.Map(location = boulder_coords)
my_map.fit_bounds(bounds=[(lat_min, lon_min), (lat_max, lon_max)])
for point in points:
folium.Marker(point).add_to(my_map)
from geopy.geocoders import Nominatim
geolocator = Nominatim(user_agent="jakbadacdane.pl")
for index, row in proper_vertices.sort_values("dist", ascending=False).iterrows():
location = geolocator.reverse("{}, {}".format(row["Lat"], row["Lon"]), timeout = 30)
#print(location.raw)
#print(row["Lat"], row["Lon"])
if location.raw == {'error': 'Unable to geocode'}:
print("Out of space")
continue
try:
is_poland = location.raw["address"]["country"] == "RP"
except KeyError:
print(location)
break
if is_poland:
print(location)
folium.Circle(radius=row["dist"]*1000, location=[row["Lat"], row["Lon"]], fill = True).add_to(my_map)
break
else:
print(location.raw["address"]["country"])
Česko Świelubie, gmina Dygowo, powiat kołobrzeski, zachodniopomorskie, 78-113, RP
#for index, row in proper_vertices.sort_values("dist").tail(n=1).iterrows():
# #print(row["dist"]*1000)
# folium.Circle(radius=row["dist"]*1000, location=[row["Lat"], row["Lon"]], fill = True).add_to(my_map)
my_map