22
votes
Techonnology and methods used:
R – tidyverse, lubridate, snakecase, sp, raster, spData, sf, leaflet, mapview, ggplot2, shiny, maps, devtools, geojsonio, rgdal, leaflet.esri, leaflet.extras
Python – s3fs, pandas, numpy, matplotlib, plotly
1. Business Understanding
Air pollution beyond the norms is a common problem in many locations. Examining the causes behind and being able to predict it would help control and reduce pollution, facilitate solving environmental issues and providing a better life to citizens.
Air pollution seems to be recurring problem in our city – Sofia. This matter can be affected by climate, topography and anthropogenic factors as well. Our aim with this project is to explain how different human activities affect the pollution and our health.
2. Data Understanding
The provided datasets consist of static data and dynamic data for 20 day in November 2016. Data is constructed as following:
- Weather Stability data – Data on the state of the upper atmosphere from radiosonde probes of Sofia
- Topological data – containg GPS coordinates for different points in Sofia.
- Official meteorological data for 20 days for Sofia city
- Industrial pollutants – 71 industrial emission devices; their location, height and annual debit of PM10.
- Construction sites – 389 construction sights; location and types of construction
- Household pollutants – 39944 households and questionnaire data determining what source of heating the households use.
- Data from official measurement stations – data from 4 stations located in Nadezhda, Hipodruma, Druzhba, Pavlovo for 20 days in November 2016.
- Lots of Miscellaneous tables (from the Methodology of the case) which we use to calculate the pollution based on the Dispersion model.
3. Data Preparation
Firstly, we decided to determine how the Dispersion model should apply to separate pollutants. In order to do so, we should use the Weather stability and wind speed from the meteorological data. Those factors helped us to determine how the pollution from separate emission devices will spread in the atmosphere of our city. Based on Weather stability we calculated the vertical temperature gradient in the upper atmosphere above Sofia for each of the 20 days. This helped us to determine that the Pasquill stability class of the atmosphere is type E (Slightly Stable) for the whole period.
Next step was to calculate sigmas – sigma y and sigma z, which represent the vertical and horizontal spread of the plum rise for each pollutant.
In order to visualize how different emission devices spread the pollution we choose a radial approach. We plot separate circles with different concentration of particles around the industrial pollutants. Those are plotted on a interactive map (.html) which you can find in our GitHub repo.
A snapshot of initial pollution with all parameters:
A snapshot of initial pollution with pollution spread in a radius of 1000 meters:
A snapshot of initial pollution with pollution spread in a radius of 10000 meters:
A snapshot of predicted pollution for 1 day ahead:
4. Modeling
To start modelling we measured the distances between industrial pollutants and official measuring stations. The data was used to define how concentration of different particles from industrial stations affects the PM10 concentrations measured in Sofia.
We succeeded to define what is the correlation between Industrial emission stations and to project what the pollution will be in the succeeding five days.
Our estimations for the industrial pollution:
Official measurement station pollution:
Due to some technical issues and inability to post our code and viz in the article you can find our work in the repo below:
Here you can find a repo with our code on R and Python.
The work continues!!! Check our progress in the repo.
Afterwards we proceed with plotting the industrial pollution so that we account for the concentration and the spread of the pollution. We load the libraries...¶
In [2]:
library(tidyverse)
library(lubridate)
library(sp)
# install.packages("sf")
# library(sf)
library(raster)
library(leaflet)
library(leaflet.esri)
library(leaflet.extras)
library(snakecase)
library(sf)
library(spData)
library(leaflet) # for interactive maps
library(mapview) # for interactive maps
library(ggplot2) # tidyverse vis package
library(shiny) # for web applications
library(maps)
library(devtools)
library(geojsonio)
library(rgdal)
and the data¶
In [9]:
# Load data
# 1) Industrial
InitIndustrial <- read.csv("industrial_pollution_v2.csv") %>%
rename(Latitude = 'X.', Longitude = 'Y.', DeviceHeight = m, AnnualPM = 't.y')
# 2) Atmosphere profile
InitAtmoProfile <- read.csv("atmosphere_profile_train.csv") %>%
mutate_at(vars(Date), funs(as_date(.)))
# 3) Construction Sites
InitConstrSites <- read.csv("construction_sites.csv")
# 4) Sofia Topography
InitSofiaTopo <- read.csv("sofia_topo.csv")
# 5) Stations
InitStations <- read.csv("stations_data_train.csv")
# 6) Weather
InitWeather <- read.csv("weather_train.csv")
# 7) HH data
#InitHH <- read.csv("household_heating.csv")
# 8) Stations
StationLoc <- read.csv("Stations.csv")
Some transformations¶
In [10]:
# Transform geolocation format in data on Industrial pollutants
Industrial <- InitIndustrial %>%
mutate_at(vars(Longitude, Latitude), funs(as.character)) %>%
mutate_at(vars(Longitude, Latitude), funs(sub('"', '', .))) %>%
mutate_at(vars(Longitude, Latitude), funs(sub('""', '\" ', .))) %>%
mutate_at(vars(Longitude, Latitude), funs(sub("''", "'", .))) %>%
mutate_at(vars(Longitude, Latitude), funs(substring(., 2))) %>%
mutate_at(vars(Longitude, Latitude), funs(str_sub(., 1, str_length(.)-2))) %>%
mutate_at(vars(Longitude, Latitude), funs(as.numeric(char2dms(., chd = "°", chm = "'", chs='"'))))
In [11]:
# Load Concentrations for day 1
# Load concentrations
newNames <- paste0(rep("Spread_", 10),
c(1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000))
ConcInd <- read.csv("C1.csv") %>%
dplyr::select(-X) %>%
rename_all(~newNames) %>%
mutate(Latitude = Industrial$Latitude,
Longitude = Industrial$Longitude)
A simple raster of topographic data¶
In [62]:
rast <- rasterFromXYZ(InitSofiaTopo[, c(2,1,3)])
plot(rast)
Mapping¶
Create icons to embed in the map
In [12]:
industrialIcon <- makeIcon(
iconUrl = "https://cdn3.iconfinder.com/data/icons/buildings-and-real-estates/30/factory-pollution-512.png",
iconWidth = 55, iconHeight = 70,
iconAnchorX = 22, iconAnchorY = 94
)
stationIcon <- makeIcon(
iconUrl = "https://cdn3.iconfinder.com/data/icons/universal-web-5/83/201-512.png",
iconWidth = 55, iconHeight = 70,
iconAnchorX = 22, iconAnchorY = 94
)
In [13]:
head(Industrial)
In [65]:
# PollutionMap <- leaflet() %>%
#
# # Base groups
# addProviderTiles(providers$Esri.WorldGrayCanvas, group = "Basic Map") %>%
# addProviderTiles(providers$Esri.WorldTopoMap, group = "Topo Map") %>%
#
# # Overlay groups
#
# # ------- industrial polluters
# addMarkers(data = Industrial, ~Longitude, ~Latitude,
# clusterOptions = markerClusterOptions(),
# icon = industrialIcon, group = "Industrial Pollution") %>%
#
# # ------ stations
# addMarkers(data = StationLoc, ~Longitude, ~Latitude,
# # clusterOptions = markerClusterOptions(),
# icon = stationIcon, group = "Official Stations") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_1000, gradient = "Greys",
# blur = 30, radius = 27,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 1000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_2000, gradient = "Greys",
# blur = 30, radius = 29,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 2000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_3000, gradient = "Greys",
# blur = 30, radius = 31,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 3000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_4000, gradient = "Greys",
# blur = 30, radius = 33,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 4000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_5000, gradient = "Greys",
# blur = 30, radius = 35,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 5000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_6000, gradient = "Greys",
# blur = 30, radius = 37,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 6000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_7000, gradient = "Greys",
# blur = 30, radius = 39,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 7000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_8000, gradient = "Greys",
# blur = 30, radius = 41,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 8000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_9000, gradient = "Greys",
# blur = 30, radius = 43,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 9000 m") %>%
#
# addHeatmap(data = ConcInd, lng = ~Longitude, lat = ~Latitude,
# intensity = ~Spread_10000, gradient = "Greys",
# blur = 30, radius = 45,
# minOpacity = 0.40,
# group = "Industrial Pollution: Spread 10000 m") %>%
#
# addHeatmap(data = StationLoc, lng = ~Longitude, lat = ~Latitude,
# intensity = ~2016-11-01, gradient = "Reds",
# blur = 30, radius = 25,
# minOpacity = 0.05,
# group = "Pollution: Official Measurements") %>%
#
# addLayersControl(
# baseGroups = c("Topo Map", "Basic Map"),
# overlayGroups = c("Official Stations", "Industrial Pollution",
# "Pollution: Official Measurements",
# "Industrial Pollution: Spread 1000 m", "Industrial Pollution: Spread 2000 m",
# "Industrial Pollution: Spread 3000 m", "Industrial Pollution: Spread 4000 m",
# "Industrial Pollution: Spread 5000 m", "Industrial Pollution: Spread 6000 m",
# "Industrial Pollution: Spread 7000 m", "Industrial Pollution: Spread 8000 m",
# "Industrial Pollution: Spread 9000 m", "Industrial Pollution: Spread 10000 m"),
# options = layersControlOptions(collapsed = FALSE)
# )
Our articles on particles would be nothing without a proper map of the data!¶
Our approach was to create interactive maps of air pollution in Sofia and for that we used R and the libraries leaflet and leaflet.extra.
Our solution includes a map showcasing a few key aspects of air pollution in Sofia:
- Official measurement stations in Sofia
- Industrial pollutants
- Pollution concentration as measured by official stations in red
- Pollution concentration for industrial units in black
The intensity of the color reflects the degree of pollution:
- The bigger the circle, the farther the pollution spreads
- The darker the colour, the higher the concentration of particles
The map is interactive and allows you, our curious readers, to explore the pollution in Sofia in several ways!
- You may choose to look at a topographic map or a basic one
- You may view all elements or only some of them
- The official stations button allows you to locate the measurement stations on the map
- Industrial pollutants may be identified as clusters of heavy polluters (depicted with counts) or as their specific locations once you zoom in
- The degree of pollution may be examined by the measure of spread - ranging from 1000 to 10,000 meters in radius from the pollutant
This map represents the first day of our historical data.
In [6]:
from IPython.display import IFrame
IFrame('https://www.datasciencesociety.net/wp-content/uploads/2019/04/Initial_Pollution.html','100%','360px')
Out[6]:
In [5]:
from IPython.display import IFrame
IFrame('https://www.datasciencesociety.net/wp-content/uploads/2019/04/Pollution_2016.11.21.html','100%','360px')
Out[5]:
https://www.datasciencesociety.net/wp-content/uploads/2019/04/Initial_Pollution.html https://www.datasciencesociety.net/wp-content/uploads/2019/04/Pollution_2016.11.21.html
In [5]:
import pandas
import plotly
import numpy as np
import matplotlib.pyplot as plt
In [6]:
atm = pandas.read_csv('https://s3.us-east-2.amazonaws.com/datacases/gd-2019/sofia_air_2/atmosphere_profile_train.csv')
construction = pandas.read_csv('https://s3.us-east-2.amazonaws.com/datacases/gd-2019/sofia_air_2/construction_sites.csv')
household = pandas.read_csv('https://s3.us-east-2.amazonaws.com/datacases/gd-2019/sofia_air_2/household_heating.csv')
industrial = pandas.read_csv('https://s3.us-east-2.amazonaws.com/datacases/gd-2019/sofia_air_2/industrial_pollution.csv', encoding='cp1252')
sof_topo = pandas.read_csv('https://s3.us-east-2.amazonaws.com/datacases/gd-2019/sofia_air_2/sofia_topo.csv')
stations = pandas.read_csv('https://s3.us-east-2.amazonaws.com/datacases/gd-2019/sofia_air_2/stations_data_train.csv')
weather = pandas.read_csv('https://s3.us-east-2.amazonaws.com/datacases/gd-2019/sofia_air_2/weather_lbsf_20161101-20161130_IP_train.csv')
In [7]:
atm.head()
Out[7]:
In [8]:
plt.figure(0)
plt.plot(atm['HGHT(m)'])
plt.legend('height')
plt.show()
plt.figure(1)
plt.plot(atm['TEMP(C)'])
plt.legend('temp')
plt.show()
In [9]:
Ht=atm['HGHT(m)']
t=atm['TEMP(C)']
atm['diffs'] = atm.groupby(['Date'])['HGHT(m)'].transform(lambda x: x.diff())
atm['diffs1'] = atm.groupby(['Date'])['TEMP(C)'].transform(lambda x: x.diff())
atm['delta'] = atm['diffs1']/atm['diffs']
In [10]:
atm1 = atm[(atm['HGHT(m)'] <= 1000)]
In [11]:
atm1
Out[11]:
In [12]:
atm2 = atm1[['Date', 'delta']]
In [13]:
atm3 = atm2.dropna().groupby('Date').mean()
In [14]:
atm3.describe()
# results show that we have only stability class E
Out[14]:
In [15]:
atm3
Out[15]:
In [16]:
# generate x and y for the calculation of C
# x=y
import numpy as np
xs = [i for i in np.arange(1000,11000,1000)]
In [17]:
# calculate sigmas only for stability class E
sig_y, sig_z = [],[]
for x in xs:
if x < 1000:
s_y = 50.5*(x**0.894)
s_z = 22.8*(x**0.678)-1.3
sig_y.append(s_y)
sig_z.append(s_z)
else:
s_y = 50.5*(x**0.894)
s_z = 55.4*(x**0.305)-34.0
sig_y.append(s_y)
sig_z.append(s_z)
In [18]:
# convert list to dataframe
sigma = pandas.DataFrame({'x': xs})
sigma['sigma_y'] = sig_y
sigma['sigma_z'] = sig_z
In [19]:
sigma
Out[19]:
In [21]:
# import transformed industrial data
industrial2 = pandas.read_csv('Industrial_Transformed.csv')
In [22]:
industrial2['DailyPM(kg)'] = industrial2['AnnualPM'].apply(lambda x: (x/365)*1000)
In [23]:
industrial2['DailyPM(g/d)'] = industrial2['DailyPM(kg)'].apply(lambda x: x*1000)
In [ ]:
#industrial2 = industrial2.drop('DailyPM(g/s)', axis=1)
In [25]:
industrial2.head()
Out[25]:
In [26]:
C_list = []
for d in range(len(atm3)):
C = np.zeros((len(industrial2),len(xs)))
C = pandas.DataFrame(C)
C.columns = xs # rename columns
u = weather.iloc[d,13]*1000/360 # average wind speed per day in meters
for index, row in industrial2.iterrows():
Q = industrial2.iloc[index, 5]
z = industrial2.iloc[index, 2]
for i,x in enumerate(xs):
#print (index, row[3], x)
#print ((Q/(2.*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-(x**2.)/(2.*sigma.iloc[i,1]**2.))*(np.exp(-(z**2.)/(2.*sigma.iloc[i,2]**2.)))**2.)
#print ((Q/(2*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-x**2/2*sigma.iloc[i,1]**2)*(np.exp(-z**2/2*sigma.iloc[i,2]**2)**2))
C.iloc[index, i] = (Q/(2.*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-(x**2.)/(2.*sigma.iloc[i,1]**2.))*(np.exp(-(z**2.)/(2.*sigma.iloc[i,2]**2.)))**2.
C = C*1000000
C_list.append(C)
In [27]:
# calculate distance to stations
Geo = pd.DataFrame({'ind':industrial2.index,'Lat':industrial2['Latitude'],'Long':industrial2['Longitude']})
Geo.head()
#official measurment
Lat1 = (42.666508,42.680558,42.669797,42.732292)
Long1=(23.400164,23.296786,23.268403,23.310972)
OM = ('STA-BG0052A','STA-BG0050A','STA-BG0073A','STA-BG0040A')
geoloc = pd.DataFrame(
{'station':OM,'latitude':Lat1,
'longitude':Long1})
geoloc.head()
#br = ({'STA-BG0052A',23.400164,42.666508},{'STA-BG0050A',23.296786,42.680558},{'STA-BG0073A',23.268403,42.669797},{'STA-BG0040A',23.310972,42.732292})
def getDist(lat1, long1, lat2, long2):
earth_radius = 6371
dLat = np.deg2rad(lat2 - lat1)
dLon = np.deg2rad(long2 - long1)
a = np.sin(dLat/2) * np.sin(dLat/2) + np.cos(np.deg2rad(lat1)) * np.cos(np.deg2rad(lat2)) * np.sin(dLon/2) * np.sin(dLon/2)
c = 2 * np.arcsin(np.sqrt(a))
d = earth_radius * c
return(d)
D = np.zeros((len(industrial2),4))
D = pd.DataFrame(D)
D.columns = OM # rename columns
# calculate distances
for i in range(4):
for j in range(len(industrial2)):
D.iloc[j,i] = getDist(geoloc.iloc[i,1], geoloc.iloc[i,2], Geo.iloc[j,1], Geo.iloc[j,2])
D.describe()
Out[27]:
In [28]:
D_list = []
D1 = round(D)
D1 = D1.astype(int)*1000
for d in range(len(atm3)): # 20 days
D_c = np.zeros((len(industrial2),4))
D_c = pandas.DataFrame(D_c)
D_c.columns = OM # rename columns
u = weather.iloc[d,13]*1000/360 # average wind speed per day in meters
for index, row in industrial2.iterrows():
Q = industrial2.iloc[index, 5]
z = industrial2.iloc[index, 2]
for i in range(4):
s_y = 50.5*(D1.iloc[index,i]**0.894)
s_z = 55.4*(D1.iloc[index,i]**0.305)-34.0
#print (index, row[3], x)
#print ((Q/(2.*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-(x**2.)/(2.*sigma.iloc[i,1]**2.))*(np.exp(-(z**2.)/(2.*sigma.iloc[i,2]**2.)))**2.)
#print ((Q/(2*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-x**2/2*sigma.iloc[i,1]**2)*(np.exp(-z**2/2*sigma.iloc[i,2]**2)**2))
D_c.iloc[index, i] = (Q/(2.*np.pi*u*s_y*s_z))*np.exp(-(D1.iloc[index,i]**2.)/(2.*s_y**2.))*(np.exp(-(z**2.)/(2.*s_z**2.)))**2.
D_c = D_c*1000000
D_list.append(D_c)
In [29]:
F = np.zeros((len(atm3),4))
F = pandas.DataFrame(F)
F.columns = OM # rename columns
for i in range(len(atm3)): # 20 days
F.iloc[i,0] = D_list[i]['STA-BG0052A'].sum()
F.iloc[i,1] = D_list[i]['STA-BG0050A'].sum()
F.iloc[i,2] = D_list[i]['STA-BG0073A'].sum()
F.iloc[i,3] = D_list[i]['STA-BG0040A'].sum()
In [30]:
# our estimations
plt.figure(0)
plt.plot(F)
plt.legend(F.columns)
plt.show()
# official measurements
plt.figure(1)
stations.plot()
plt.show()
In [31]:
from sklearn import linear_model
# define predictors
expl = pd.DataFrame(F, columns=['STA-BG0052A'])
# define target
target = pd.DataFrame(stations, columns=['STA-BG0052A'])
In [32]:
# fit model
lm = linear_model.LinearRegression()
model = lm.fit(expl,target)
In [33]:
# see R2
print (lm.score(expl,target))
In [34]:
# see intercept and coefficient
print (lm.intercept_)
print(lm.coef_)
In [35]:
predictions = lm.predict(expl)
predictions
Out[35]:
In [37]:
plt.scatter(target,predictions)
plt.show()
In [38]:
# save to csv
C_list[0].to_csv('C1.csv')
In [39]:
# import test sets
weather_t = pandas.read_csv('weather_lbsf_20161101-20161130_IP_test.csv')
stations_t = pandas.read_csv('stations_data_test.csv')
In [40]:
# calculate next 5 points
D_test = []
#D1 = round(D)
#D1 = D1.astype(int)*1000
for d in range(len(weather_t)): # 5 days
D_t = np.zeros((len(industrial2),4))
D_t = pandas.DataFrame(D_t)
D_t.columns = OM # rename columns
u = weather_t.iloc[d,13]*1000/360 # average wind speed per day in meters
for index, row in industrial2.iterrows():
Q = industrial2.iloc[index, 5]
z = industrial2.iloc[index, 2]
for i in range(4):
s_y = 50.5*(D1.iloc[index,i]**0.894)
s_z = 55.4*(D1.iloc[index,i]**0.305)-34.0
#print (index, row[3], x)
#print ((Q/(2.*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-(x**2.)/(2.*sigma.iloc[i,1]**2.))*(np.exp(-(z**2.)/(2.*sigma.iloc[i,2]**2.)))**2.)
#print ((Q/(2*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-x**2/2*sigma.iloc[i,1]**2)*(np.exp(-z**2/2*sigma.iloc[i,2]**2)**2))
D_t.iloc[index, i] = (Q/(2.*np.pi*u*s_y*s_z))*np.exp(-(D1.iloc[index,i]**2.)/(2.*s_y**2.))*(np.exp(-(z**2.)/(2.*s_z**2.)))**2.
D_t = D_t*1000000 # this result is not devided by 1225g/m3 (air density)
D_test.append(D_t)
In [41]:
# final results from test set
F_t = np.zeros((len(weather_t),4))
F_t = pandas.DataFrame(F_t)
F_t.columns = OM # rename columns
for i in range(len(weather_t)): # 5 days
F_t.iloc[i,0] = D_test[i]['STA-BG0052A'].sum()
F_t.iloc[i,1] = D_test[i]['STA-BG0050A'].sum()
F_t.iloc[i,2] = D_test[i]['STA-BG0073A'].sum()
F_t.iloc[i,3] = D_test[i]['STA-BG0040A'].sum()
In [42]:
F_t
Out[42]:
In [43]:
# export csv
D_test[0].to_csv('industrial_polluters_test.csv') # predictions for each polluter on day 21
F_t.to_csv('final_t.csv') # predictions for official stations
In [44]:
# our estimations for the test set
plt.figure(0)
plt.plot(F_t)
plt.legend(F_t.columns)
plt.show()
# official measurements from the test set
plt.figure(1)
stations_t.plot()
plt.show()
In [45]:
print(lm.predict(pd.DataFrame(F_t, columns=['STA-BG0052A'])))
In [46]:
print(stations_t['STA-BG0052A'])
In [47]:
# get results for different distances for the test set
C1_list = []
for d in range(len(weather_t)):
C = np.zeros((len(industrial2),len(xs)))
C = pandas.DataFrame(C)
C.columns = xs # rename columns
u = weather_t.iloc[d,13]*1000/360 # average wind speed per day in meters
for index, row in industrial2.iterrows():
Q = industrial2.iloc[index, 5]
z = industrial2.iloc[index, 2]
for i,x in enumerate(xs):
#print (index, row[3], x)
#print ((Q/(2.*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-(x**2.)/(2.*sigma.iloc[i,1]**2.))*(np.exp(-(z**2.)/(2.*sigma.iloc[i,2]**2.)))**2.)
#print ((Q/(2*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-x**2/2*sigma.iloc[i,1]**2)*(np.exp(-z**2/2*sigma.iloc[i,2]**2)**2))
C.iloc[index, i] = (Q/(2.*np.pi*u*sigma.iloc[i, 1]*sigma.iloc[i,2]))*np.exp(-(x**2.)/(2.*sigma.iloc[i,1]**2.))*(np.exp(-(z**2.)/(2.*sigma.iloc[i,2]**2.)))**2.
C = C*1000000
C1_list.append(C)
In [48]:
# export to csv
C1_list[0].to_csv('C1_t.csv')
4 thoughts on “Articles on Particles: Air Sofia Case by Wonder Gang”
Great job! Really good map (with the blur visualization it rocks). Your paper shows a great combo of R and Python in one task. Keep working and calculate the other factors.
Thank you for the encouraging words. We will continue working on the case.
I do like of approach of “blurring” over map for pollution spreading. It would be even better (I know you didn’t had time) to take into account surounding buildings and how they impact pollution spread to get more custom “shape” of spread instead of combination of circles but this visually tells amazing story. Keep good work of “story” presentation of other factors as well.
Thank you for your comment. Actually, we were thinkinkg to use topogrphic data to improve the shape of the spread visualization. Looking at the surrounding objects is definetly a good idea!