{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Emission factor for each powerplant in Puerto Rico\n## This notebook will compute below two equations with linear regression model\n#### Simplified Emissions Factor = (NO2 emission from Powerplants) / (Quanity of electricity generated)\n#### NO2 emission from Powerplants = NO2 emission - NO2 emissions from rest(Nightlight + Precipitation + Road + Landuse)\n## In short, the purpose of this notebook is to disaggregate NO2 Emission caused from powerplant and caused from else than powerplant\n## At the end of notebook, *Each powerplant* will have Simplified Emission Factor for each month.\n\n\n#### Note that the output of this notebook are two (you can download at the end of the page):\nMunicipality.csv: Data aggregated by municipality level. Nothing scaled, all the values are average in municipality level\n\nresult.csv: Final output, all the list of plants and its Simplified Emission Factors.\n\n\n\n# 1. Data step\nImport library"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"### Library import\nlibrary(tidyverse) \nlibrary(rgdal)\nlibrary(ggplot2)\nlibrary(raster)\nlibrary(viridis)\nlibrary(stars)\nlibrary(data.table)\nlibrary(sf)\nlibrary(RColorBrewer)\nlibrary(raster)\nlibrary(scales)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Import NO2 emission data, 387 images, y-axis is flipped. "},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"#### Import NO2 emission data, 374 images, y-axis is flipped. Gather all in one variable. \nfiles = list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2\",pattern=\"*.tif$\", full.names=TRUE)\nrs <-stack(files,quick=TRUE) #bands=1\nrs1 <- sum(rs,na.rm = TRUE)\nrs1 <- flip(rs1, direction = \"y\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Import powerplants data to R"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"#### powerplants import\nplants=read.csv('../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv',stringsAsFactors = FALSE)\n\nplants$lon=0\nplants$lat=0\n\nfor (i in 1:dim(plants)[1])\n{\nplants$lon[i]=strsplit(gsub(\"\\\\[|\\\\]\", \"\", gsub(\"[\\\\(\\\\)]\", \"\", regmatches(plants[,'.geo'], gregexpr(\"\\\\[.*?\\\\]\", plants[,'.geo']))[[i]])),split=',')[[1]][1]\nplants$lat[i]=strsplit(gsub(\"\\\\[|\\\\]\", \"\", gsub(\"[\\\\(\\\\)]\", \"\", regmatches(plants[,'.geo'], gregexpr(\"\\\\[.*?\\\\]\", plants[,'.geo']))[[i]])),split=',')[[1]][2]\n    }\n\n\nplants_geo <- st_as_sf(plants, coords = c(\"lon\", \"lat\"),crs = 4326, agr = \"constant\")\n\n#ggplot() +geom_raster(aes(x = x, y = y, fill = layer,alpha=layer),data=rs2)+ scale_alpha(range =  c(0.2, 1.5), guide = \"none\") + geom_sf(data = world,alpha=0,col='black')+geom_sf(data=plants_geo,aes(size=capacity_mw))+ theme(legend.position = \"none\")\n\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Distribute all NO2 emission in municipality level "},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"##### Import municipality in Puerto Rico \n## Source: GADM https://gadm.org/maps/PRI.html\npri <- st_read('../input/puerto/gadm36_PRI_shp/gadm36_PRI_1.shp')\n\n### Aggregate total NO2 emission by municipality level\npri$NO2 <- extract(rs1, pri,fun=sum)\n\n### Aggregate NO2 emission per month by municipality level\n\n## Assign the data for each months\nfrequency <- table(substr(list.files(path = '../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2'),9,14))\nselect <- cumsum(frequency)\nbegin <- select-frequency+1\nend <- select\n\n## A loop to fill in NO2 emission per each month, each loop is each month\n\nfor (i in 1:12)\n{\n    files = list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2\",pattern=\"*.tif$\", full.names=TRUE)[begin[i]:end[i]]\n    rs <-stack(files,quick=TRUE) #bands=1\n    rs1 <- sum(rs,na.rm = TRUE)\n    rs1 <- flip(rs1, direction = \"y\")\n    pri[,paste0('Y',names(select)[i])] <- extract(rs1, pri,fun=mean)\n    pri[,paste0('SY',names(select)[i])] <- extract(rs1, pri,fun=sum)\n    #plot <- ggplot()+geom_sf(aes_string(fill=paste0('Y',names(select)[i])),data=pri)+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")+ggtitle(paste0('Y',names(select)[i]))\n    #print(plot)\n}\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Import nightlight data"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"######### import nightlight data from NOAA, it is one file per month. \n\nnight_files <- c(paste0('../input/nightlight-pri/',list.files(path = '../input/nightlight-pri')),paste0('../input/nightlight-pri-2019/',list.files(path = '../input/nightlight-pri-2019')))\n\n\nfor (i in 1:12)\n{untar(night_files[i])\n files = list.files(pattern='rade9h.tif$') \n rs1 <- raster(files)\n pri[,paste0('N',names(select)[i])] <- extract(rs1, pri,fun=mean)\n pri[,paste0('SN',names(select)[i])] <- extract(rs1, pri,fun=sum)\n #plot <- ggplot()+geom_sf(aes_string(fill=paste0('Y',names(select)[i])),data=pri)+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")+ggtitle(paste0('N',names(select)[i]))\n #print(plot) \n file.remove(list.files()[-1]) #Nightlight file is very big whereas kaggle's capacity limit is only 5gb, each loop keep remove all the files except note itself\n }\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Import wind and precipitation\n\n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"##gldas doesn't cover whole Puerto Rico \n#list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/gfs\",pattern=\"*.tif$\", full.names=TRUE)\n\n### Assign gfs data per month\nfrequency_noaa <- table(substr(list.files(path = '../input/ds4g-environmental-insights-explorer/eie_data/gfs'),5,10))\nselect_noaa <- cumsum(frequency_noaa)\nbegin_noaa <- select_noaa-frequency_noaa+1\nend_noaa <- select_noaa\n\n\n## loop for uwind\nfor (i in 1:12)\n{\n    files = list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/gfs\",pattern=\"*.tif$\", full.names=TRUE)[begin_noaa[i]:end_noaa[i]]\n    rs <-stack(files,bands=4)\n    rs1<-sum(rs,na.rm = TRUE)\n    pri[,paste0('uwind',names(select_noaa)[i])] <- extract(rs1, pri,fun=mean)\n    pri[,paste0('sum_uwind',names(select_noaa)[i])] <- extract(rs1, pri,fun=sum)\n}\n## loop for vwind\nfor (i in 1:12)\n{\n    files = list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/gfs\",pattern=\"*.tif$\", full.names=TRUE)[begin_noaa[i]:end_noaa[i]]\n    rs <-stack(files,bands=5) \n    rs1<-sum(rs,na.rm = TRUE)\n    pri[,paste0('vwind',names(select_noaa)[i])] <- extract(rs1, pri,fun=mean)\n    pri[,paste0('sum_vwind',names(select_noaa)[i])] <- extract(rs1, pri,fun=sum)\n}\n## loop for precipitation\nfor (i in 1:12)\n{\n    files = list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/gfs\",pattern=\"*.tif$\", full.names=TRUE)[begin_noaa[i]:end_noaa[i]]\n    rs <-stack(files,bands=6) \n    rs1<-sum(rs,na.rm = TRUE)\n    pri[,paste0('rain',names(select_noaa)[i])] <- extract(rs1, pri,fun=mean)\n    pri[,paste0('sum_rain',names(select_noaa)[i])] <- extract(rs1, pri,fun=sum)\n}\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"import OSM"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"pri_road <- st_read('../input/puerto-road/puerto_road/road_pri.shp')\n#ggplot()+geom_sf(aes(col='blue'),data=pri_road)+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")+ggtitle(paste0('Y',names(select)[i]))\n\nints <- st_intersection(pri_road, pri)\nroad <- tapply(st_length(ints), ints$NAME_1,sum)\npri$road[match(names(road),pri$NAME_1)] = road\n\n#ggplot()+geom_sf(aes(fill=road),data=pri)+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")+ggtitle(paste0('Y',names(select)[i]))\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Make a complete table"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"pri_industry <- st_read('../input/puerto-road/industrial_osm/industrial_osm.shp')\n#ggplot()+geom_sf(aes(col='red'),data=pri_industry)+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")+ggtitle(paste0('Y',names(select)[i]))\nints_ind <- st_intersection(pri_industry, pri)\nindustry <- tapply(st_area(ints_ind), ints_ind$NAME_1,sum)\npri$industry[match(names(industry),pri$NAME_1)] = industry\n\n#ggplot()+geom_sf(aes(fill=industry),data=pri)+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")+ggtitle(paste0('Y',names(select)[i]))\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Import Data from EIA.gov (From Vincent Larmet's notebook https://www.kaggle.com/vlarmet/an-r-notebook-for-no2-emission-factor#Exploratory-Data-Analysis)"},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"df18<-read.csv(\"../input/electricity-production-in-puerto-rico/EIA923_Schedules_puerto_rico_2018.csv\",stringsAsFactors=FALSE)\ndf17<-read.csv(\"../input/electricity-production-in-puerto-rico/EIA923_Schedules_puerto_rico_2017.csv\",stringsAsFactors=FALSE)\n\n#2017\nvars<-colnames(df17)[c(4,16,20:91)]\ndf17<-data.table::melt(df17[,vars],id.vars=c(\"Plant.Name\",\"AER.Fuel.Type.Code\"))\ndf17$value<-as.numeric(df17$value)\ndf17$variable<-as.character(df17$variable)\ndf17$month<-sapply(strsplit(df17$variable,\"\\\\.\"), function(x) x[2])\ndf17$variable<-sapply(strsplit(df17$variable,\"\\\\.\"), function(x) x[1])\nelm2<-setDT(df17)[variable==\"Netgen\",.(value=sum(value,na.rm=T)),.(month,AER.Fuel.Type.Code)]\nlabs<-unique(elm2$month)\nelm2$month<-(1:12)[match(elm2$month,labs)]\n\n                      \n#2018\nvars<-colnames(df18)[c(4,16,20:91)]\ndf18<-data.table::melt(df18[,vars],id.vars=c(\"Plant.Name\",\"AER.Fuel.Type.Code\"))\ndf18$value<-as.numeric(df18$value)\ndf18$variable<-as.character(df18$variable)\ndf18$month<-sapply(strsplit(df18$variable,\"\\\\.\"), function(x) x[2])\ndf18$variable<-sapply(strsplit(df18$variable,\"\\\\.\"), function(x) x[1])\nelm<-setDT(df18)[variable==\"Netgen\",.(value=sum(value,na.rm=T)),.(month,AER.Fuel.Type.Code)]\nlabs<-unique(elm$month)\nelm$month<-(1:12)[match(elm$month,labs)]\n                      \n#2017\n\nelm<-df17[variable==\"Netgen\" & AER.Fuel.Type.Code %in% c(\"COL\",\"RFO\",\"DFO\",\"NG\"),\n          .(value=sum(value,na.rm=T)),\n          .(AER.Fuel.Type.Code,Plant.Name)]\n\n\n#2018\nelm<-df18[variable==\"Netgen\" & AER.Fuel.Type.Code %in% c(\"COL\",\"RFO\",\"DFO\",\"NG\"),\n          .(value=sum(value,na.rm=T)),\n          .(AER.Fuel.Type.Code,Plant.Name)]\n\n                      eltot<-df18[variable==\"Netgen\",c(\"Plant.Name\",\"AER.Fuel.Type.Code\",\"month\",\"value\"),with=FALSE]\neltot<-eltot[,.(value=sum(value,na.rm=F)),.(Plant.Name,AER.Fuel.Type.Code,month)]\neltot2<-df17[variable==\"Netgen\",c(\"Plant.Name\",\"AER.Fuel.Type.Code\",\"month\",\"value\"),with=FALSE]\neltot2<-eltot2[,.(value=sum(value,na.rm=F)),.(Plant.Name,AER.Fuel.Type.Code,month)]\neltot<-merge(eltot,eltot2,by=c(\"Plant.Name\",\"AER.Fuel.Type.Code\",\"month\"),all.x=TRUE)\ncolnames(eltot)[4:5]<-c(\"y2018\",\"y2017\")\n\neltot$y2019<-eltot$y2018\neltot$y2019[eltot$month %in% c(\"January\",\"February\",\"March\",\"April\",\"May\",\"June\")]<-rowMeans(eltot[month %in% c(\"January\",\"February\",\"March\",\"April\",\"May\",\"June\"),c(\"y2018\",\"y2017\"),with=FALSE],na.rm=T)\neltot$y2019[is.na(eltot$y2019)]<-eltot$y2017[is.na(eltot$y2019)]\nhead(eltot)\nelm3<-eltot[,.(value=sum(y2019)),.(month,AER.Fuel.Type.Code)]\nlabs<-c(\"July\",\"August\",\"September\",\"October\",\"November\",\"December\",\"January\",\"February\",\"March\",\"April\",\"May\",\"June\")\nelm3$month<-(1:12)[match(elm3$month,labs)]\n\nto_merge <- data.frame(\n                      Plant.Name=c(\"Cambalache Plant\",\"HIMA San Pablo - Caguas\",\"Daguao\",\"Culebra\",\"AES ILUMINA\",\"AES Puerto Rico\",\"Jobos\",\"Janssen Ortho LLC\",\"Humacao Solar Project, LLC\",\"Yabucoa\",\"Oriana Energy Hybrid\",\"San Fermin Solar Farm\",\"Mayaguez Plant\",\"Punta Lima Wind Farm\",\"Costa Sur Plant\",\"EcoElectrica\",\"Aguirre Plant\",\"Central San Juan Plant\",\"Pattern Santa Isabel LLC\",\"Holsum de Puerto Rico, Inc.\",\"Palo Seco Plant\",\"Caonillas\",\"Vega Baja\",\"Vieques\"),\n                      Municipality=c(\"Arecibo\",\"Caguas\",\"Ceiba\",\"Culebra\",\"Guayama\",\"Guayama\",\"Guayama\",\"Gurabo\",\"Humacao\",\"Humacao\",\"Isabela\",\"Loiza\",\"Mayaguez\",\"Naguabo\",\"Penuelas\",\"Penuelas\",\"Salinas\",\"San Juan\",\"Santa Isabel\",\"Toa Baja\",\"Toa Baja\",\"Utuado\",\"Vega Baja\",\"Vieques\")\n                      )\n \n               \njoined <- merge(eltot,to_merge,by='Plant.Name')\ngeneration_2019 <- aggregate(y2019~Municipality+month ,data=joined,FUN=sum)\ngeneration_2018 <- aggregate(y2018~Municipality+month ,data=joined,FUN=sum)\njoined$Month <- match(joined$month,month.name)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 2. Seasonality "},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"names(pri)\nresult_table <- pri[,c(4,12:dim(pri)[2])]\nresult_table <- as.data.frame(result_table)\nresult_table$geometry <- NULL\n#names(result_env)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"#library(reshape2)\nmelted_table <- melt(result_table,id.vars='NAME_1')\nmelted_table <- c(melted_table,data.frame(do.call('rbind', strsplit(as.character(melted_table$variable),'20',fixed=TRUE))))\nmelted_table <- data.frame(melted_table)\nTotal_table <- melted_table[which(melted_table$X1%in%c('SY','SN','sum_uwind','sum_vwind','sum_rain')),]\nTotal_table$value <- abs(Total_table$value)\nTotal_result <- aggregate(data=Total_table,value~variable+X1+X2,sum)\nTotal_result$month <- substr(Total_result$X2,3,4)\nTotal_result$value <- round(Total_result$value,3)\n#dcast(Total_result\n#ggplot(data=nb) + geom_bar(aes(x=Month,y=Total_No2),stat='identity', fill=\"steelblue\")+ggtitle('Total No2 emission each month in Puerto Rico')+ theme_minimal()\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This section is to check the seasonality for the variable I am going to use for the analysis. \n\n#### 1. Total NO2 emission\n\nIt is the most important variable. From the below bar graph, you can see it is significantly high for March to June. "},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"Total_result$value <- factor(Total_result$value)\nggplot(data=subset(Total_result,X1=='SY')) + geom_bar(aes(x=month,y=value),stat='identity',fill=\"steelblue\")+ggtitle('Total NO2 emission each month in Puerto Rico')+theme_minimal()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 2. Total nightlight \n\nSince NO2 emission is mainly due to the burning fuel, Nightlight must be one of key variable to track down. \nNightlight also look very seasonal, March to June high, July to November low. "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ggplot(data=subset(Total_result,X1=='SN')) + geom_bar(aes(x=month,y=value),stat='identity', fill=\"steelblue\")+ggtitle('Total nightlight each month in Puerto Rico')+ theme_minimal()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 3. Total Precipitation \n\n\nSince rain affects a lot for floating NO2 emission and also makes nightlight blurry, it is important variable as well. \n\nYou can see overall in Puerto Rico, May to October there's more rain than the other months. It corresponds general knowledge on Puerto Rico that the rainy season is between May and October. \n\nIt will be discussed later, but depends on the region, rainy season is different. I will go more detail level of precipitation in later chapter."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ggplot(data=subset(Total_result,X1=='sum_rain')) + geom_bar(aes(x=month,y=value),stat='identity', fill=\"steelblue\")+ggtitle('Total precipitation in Puerto Rico')+ theme_minimal()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 4. Strength of wind\n\nFrom Global Forecast System (GFS) dataset, there are two bands for wind speed. \n\nthere are two indicators for wind, U component and V component. (m/s,wind 10m above ground)\n\nThe meteorological convention for winds is that U component is positive for a west to east flow (eastward wind) and the V component is positive for south to north flow (northward wind).\n\nYou can observe that overall in Puerto Rico, \n\n #### a. East & West wind is much stronger than North & South Wind.\n #### b. Wind is significantly stronger between April - August, which pretty much overlap with rainy season."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ggplot(data=subset(Total_result,X1=='sum_uwind')) + geom_bar(aes(x=month,y=value),stat='identity', fill=\"steelblue\")+ggtitle('North & South Wind')+ theme_minimal()\nggplot(data=subset(Total_result,X1=='sum_vwind')) + geom_bar(aes(x=month,y=value),stat='identity', fill=\"steelblue\")+ggtitle('East & West Wind')+ theme_minimal()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 5. Estimated electricity production (From Vincent Larmet's notebook, https://www.kaggle.com/vlarmet/an-r-notebook-for-no2-emission-factor#Exploratory-Data-Analysis)\n\nSince dataset from EIA doesn't offer the geo-location of each plant, I input manually the geo-location of each plant.\nHowever EIA offers monthly estimate of electricity generation and also people pointed out that there's some problem in GGPD dataset, it is needed.\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ggplot(elm3,aes(x=month,y=value,colour=factor(AER.Fuel.Type.Code)))+\n  geom_line(size=1)+\n  scale_x_continuous(breaks=1:12,label=labs)+\n  theme_bw()+\n  theme(axis.text.x = element_text(angle = 45, hjust = 1))+\n  labs(colour=\"Fuel type\",y=\"MWh\")+\n  ggtitle(\"Estimated production from 2018-07 to 2019-06\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3. Remarks on each database\n"},{"metadata":{},"cell_type":"markdown","source":"## a. NO2 emission and Nightlight (Year average)\n![](https://imgur.com/imKBWis.png)\n![](https://imgur.com/9SWV4Tj.png)\n![](https://imgur.com/l8xbM3j.png)\n\n### Quick impression: \n#### 1. NO2 emission looks very smoothed, distributed linear, layered. You can't see the island.\n#### 2. Nightlight is very sharp and specific, you can immediately distinguish landarea\n#### 3. There are parts overlapping, high NO2 emission and high nightlight (Camuy, Morovis), low NO2 emission and low nightlight (Utuado)\n#### 4. There are parts dicrepency, low NO2 emission and high nightlight (Juana diaz)\n#### 5. However, there's never any case high NO2 emission but low nightlight\n#### 6. Since NO2 emission is due to the burning fuel, nightlight is necessary but not all nightlight may cause NO2 emission\n\nData source on Nightlight:\nThe Earth Observations Group (EOG) at NOAA/NCEI is producing a version 1 suite of average radiance composite images using nighttime data from the Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB).\n\nhttps://eogdata.mines.edu/download_dnb_composites.html\n\n## b. Global Forecast System (GFS) dataset\n\n### Rain\n\nThe total amount of precipiation (kg/m^2), is aggregated by municipality level.\n\n## c. OpenStreetMap (OSM)\n\nThe detail explanation of what extracted are in following link. Dots on the map are powerplants from EIA\nhttps://tools.wmflabs.org/query2map/osmhack.php?lat=18.228047&lon=-66.389008&name=Puerto+Rico\n\n\n### Road\nOnly the main road is extracted. For the analysis, the total length (km) of the road is aggregated by municipality level\nFor road following query are used: \"highway\" IN ('motorway','primary','secondary','trunk')\n![](https://imgur.com/0R2sHkh.png)\n\n\n### Industrial zone surface\nIn OSM, they offer the land type, the main NO2 emission is also due to contruction and industrial area.\nFor Industrial zone follwoing query are used: \"landuse\" IN ('construction','industrial')\n![](https://imgur.com/yIbCXN4.png)\n"},{"metadata":{},"cell_type":"markdown","source":"# 4. Municipality table"},{"metadata":{},"cell_type":"markdown","source":"#### This aggregate all the data in municipality level by month.\n#### Each municipality, has NO2 emission level, nightlight level, Precipitation level, Industry_surface, Road_length. \n#### Note that Industry_surface and Road_length is not a time-variant variable in this dataset\n#### You can download easily on the output folder, Municipality.csv file."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"generation_2018$Month <- match(generation_2018$month,month.name)\ngeneration_2019$Month <- match(generation_2019$month,month.name)\nnames(generation_2018)[3] <- 'electricity_gen'\nnames(generation_2019)[3] <- 'electricity_gen'\ngeneration <- rbind(subset(generation_2018,Month%in%c(7,8,9,10,11,12)),subset(generation_2019,Month%in%c(1,2,3,4,5,6)))\ngeneration$month=NULL\ngeneration$electricity_gen[which(generation$electricity_gen<1)]=1\ngeneration$electricity_gen <- log10(generation$electricity_gen)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"monthly_table <- melted_table[which(melted_table$X1%in%c('Y','N','uwind','vwind','rain')),] #'road','industry'\n#monthly_table\nMonthly_result <-  dcast(monthly_table,NAME_1+X2~X1)\nMonthly_result$month <- substr(Monthly_result$X2,3,4)\nMonthly_result$wind <- abs(Monthly_result$uwind)+abs(Monthly_result$vwind)\nMonthly_result <- Monthly_result[,c('NAME_1','N','rain','wind','Y','month')]\n\n\n###### prepare road and industry to merge since it is not time-variant vairables\nosm <- melted_table[which(melted_table$variable%in%c('road','industry')),]\nosm$value[is.na(osm$value)]=0\nosm <- dcast(osm,NAME_1~X1)\n\nComplete <- merge(Monthly_result,osm)\nnames(Complete) <- c('Municipality','Nightlight','Rain','Wind','NO2','Month','Industry_surface','Road_length')\nComplete <- Complete[,c('Municipality','Month','NO2','Nightlight','Rain','Wind','Industry_surface','Road_length')]\nComplete$Month <- as.numeric(Complete$Month)\nComplete <- merge(Complete,generation,by=c('Month','Municipality'),all.x=TRUE)\ngeneration$electricity_gen <- 10^generation$electricity_gen\nComplete$electricity_gen[is.na(Complete$electricity_gen)]=0\n\nwrite.csv(Complete,'Municipality.csv')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 5. Regression"},{"metadata":{},"cell_type":"markdown","source":"#### In this section, based on municipality table, I run a linear regression\n#### The purpose of regression is separating NO2 emission coming from plant itself such as traffic,burning,human activities,wind,rain.\n#### The sequence of months can make a nice panel data too, but here it is pool based analysis \n\n#### The equation I run for linear regression is below\n### NO2 emission = NO2 emissions Powerplants + NO2 emissions rest(Nightlight + Precipitation + Road + Landuse)\n\n"},{"metadata":{},"cell_type":"markdown","source":"First, scaling the variables for the regression (for 0 to 100)\nAbove is the boxplot before the scaling \nand below is the boxplot after the scaling"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"boxplot(Complete[,c(3:9)])\ndat <- Complete[,c(3:9)]\nregression <- data.frame(lapply(dat, function(x) scale(x, center = FALSE, scale = max(x, na.rm = TRUE)/100)))\n#regression <- rescale(Complete[,c(3:9)],to = c(0,1))[,]\nregression <- data.frame(regression)\nboxplot(regression)\nnames(regression)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"reg_result <- lm(NO2~Nightlight+Rain+Industry_surface+Road_length+electricity_gen,data=regression)\nsummary(reg_result)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 6. Result"},{"metadata":{},"cell_type":"markdown","source":"Based on the above linear regression result, I just take the NO2 emission produced by electricity_gen\nThe method is following. \n1. since the variables used in the model are scaled, what we extract is just the ratio. (electricity_gen)/NO2 will give you how much NO2 of the municipality is due to electricity generation\n2. Use the ratio calculated from 1. from the NO2 of the municiaplity multiply actual (before scaled) NO2\n3. Divide the result from 2. to each powerplant for the case when there are more than one plant in the municipality proportionally with electricity generation\n4. Divde the result from 3. by generated electricity \n5. The output of 4. is the Simplified Emissions Factor of the powerplant"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#Step 1.\nNO2_from_electricity <- coefficients(reg_result)[6]*regression$electricity_gen\nRatio_NO2_electricity  <- NO2_from_electricity/regression$NO2\n#Step 2.\nRatio_NO2_electricity  <- NO2_from_electricity/regression$NO2\nAbs_NO2_electricity  <- Ratio_NO2_electricity*Complete$NO2\n#Step 3.\nto_merge <- cbind(Complete[,c('Month','Municipality','NO2')],Ratio_NO2_electricity,Abs_NO2_electricity)\nto_merge <- merge(joined,to_merge,by=c('Month','Municipality'))\nto_merge$electricity_reg <- 0\nto_merge$electricity_reg[which(to_merge$Month%in%c(1,2,3,4,5,6))] <- to_merge$y2019[which(to_merge$Month%in%c(1,2,3,4,5,6))]\nto_merge$electricity_reg[which(to_merge$Month%in%c(7,8,9,10,11,12))] <- to_merge$y2019[which(to_merge$Month%in%c(7,8,9,10,11,12))]\nto_merge <- merge(generation,to_merge,by=c('Month','Municipality'))\nto_merge$ratio_municipality <- to_merge$electricity_reg/to_merge$electricity_gen\nto_merge$Abs_NO2_elec_plant <- to_merge$Abs_NO2_electricity*to_merge$ratio_municipality\n#Step 4.\nto_merge$SEF <- to_merge$Abs_NO2_elec_plant/to_merge$electricity_reg\nto_merge <- to_merge[which(to_merge$electricity_reg>0),]\n##### Keep only rhe result \nresult <- to_merge[,c('Month','Municipality','Plant.Name','electricity_gen','NO2','Abs_NO2_elec_plant','SEF')]\n##### For easier read multiply 10^10 on SEF\nresult$SEF <- result$SEF*10^10","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### So here's final result of my notebook. \n#### Each plant has each Simplified Emission Factor each month. \n#### Here's explanation of each variable: \n#### Municipality : Name of Municipality where the plant belongs to\n#### Electricity_gen : Generated electericity for corresponding month (municipality level)\n#### NO2 : The amount of emission of NO2 (municipality level)\n#### Abs_NO2_elec_plant : Absolute amount of NO2 emission made by plant (plant level)\n#### SEF : Simplified Emission Factor, which is NO2 emission from the plant divided by amount of generated electricity. (note that it is 10^10 multiplied for easy reading)\n\nYou can download easily the csv file generated from below."},{"metadata":{"trusted":true},"cell_type":"code","source":"result\nwrite.csv(result,'result.csv')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":4}