In [1]:
%matplotlib inline
import requests
import pandas as pd
import numpy as np
import json
In [2]:
import time
In [3]:
r = requests.get('http://api.gios.gov.pl/pjp-api/rest/station/findAll')
allStations = pd.io.json.json_normalize(r.json())
In [4]:
pollutants = ["PM2.5", "PM10", "C6H6", "NO2", "SO2", "O3", "CO"]
In [5]:
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")
In [6]:
allStations["gegrLat"] = allStations["gegrLat"].astype("float")
allStations["gegrLon"] = allStations["gegrLon"].astype("float")
In [7]:
%%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
In [8]:
allStations2 = pd.concat([allStations, allStations['sensors'].str.join('|').str.get_dummies()], axis=1, sort=False)
allStations2[pollutants].sum()
Out[8]:
PM2.5     64
PM10     136
C6H6      54
NO2      144
SO2      126
O3       101
CO        77
dtype: int64
In [9]:
#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
In [10]:
points
Out[10]:
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]])
In [11]:
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
Out[11]:
(54.560836, 49.293564, 23.155869, 15.008175)
In [12]:
from scipy.spatial import voronoi_plot_2d
from scipy.spatial import Voronoi

vor = Voronoi(points)
In [13]:
vertices = pd.DataFrame(vor.vertices, columns = ["Lat", "Lon"])
In [14]:
import matplotlib.pyplot as plt
voronoi_plot_2d(vor)
plt.show()
In [15]:
vertices[(vertices["Lat"] < lat_max) & (vertices["Lat"] > lat_min)]
Out[15]:
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

In [16]:
proper_vertices = pd.DataFrame(vertices[(vertices["Lon"] < lon_max) & 
         (vertices["Lon"] > lon_min) & 
         (vertices["Lat"] < lat_max) & 
         (vertices["Lat"] > lat_min)], copy = True)
In [17]:
import geopy.distance
In [18]:
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)
In [19]:
proper_vertices.sort_values("dist").tail(n=1)
Out[19]:
Lat Lon dist
4 49.308101 16.986539 128.84905
In [20]:
import folium
from folium.plugins import MarkerCluster
In [21]:
boulder_coords = [lat_min + (lat_max - lat_min)/2, lon_min + (lon_max - lon_min)/2]
In [22]:
my_map = folium.Map(location = boulder_coords)
my_map.fit_bounds(bounds=[(lat_min, lon_min), (lat_max, lon_max)])
In [23]:
for point in points:
    folium.Marker(point).add_to(my_map)
In [24]:
from geopy.geocoders import Nominatim
geolocator = Nominatim(user_agent="jakbadacdane.pl")
In [25]:
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
In [26]:
#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)
In [27]:
my_map
Out[27]:
In [ ]: