{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Project overview\n\nDevelop a methodology to calculate an average historical emissions factor of electricity generated for a sub-national region, using remote sensing data and techniques.\n\nThe Environmental Insights Explorer team at Google is keen to gather insights on ways to improve calculations of global emissions factors for sub-national regions. The ultimate goal of this challenge is to test if calculations of emissions factors using remote sensing techniques are possible and on par with calculations of emissions factors from current methodologies."},{"metadata":{},"cell_type":"markdown","source":"# Exploratory Data Analysis"},{"metadata":{},"cell_type":"markdown","source":"In this section, i will explore various Puerto Rico characteristics :\n* NO2, cloud and aerosol absorption index data from TROPOMI\n* electricity production\n* population density\n* traffic from commuting data\n* San Juan hotspot\n* Wind data"},{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"#Install useful packages if not installed\npkg<-c(\"raster\",\n       \"sf\",\n       \"gstat\",\n       \"data.table\",\n       \"pracma\",\n       \"tmap\",\n       \"cppRouting\",\n       \"stars\",\n      \"doParallel\",\n      \"repr\")\npkg<-pkg[!(pkg %in% installed.packages()[,\"Package\"])]\nif (length(pkg)>0) install.packages(pkg)\n\n#load packages\nlibrary(tidyverse) \nlibrary(raster) #raster manipulation\nlibrary(sf) #spatial object manipulation\nlibrary(gstat) #geostatistics (kriging, interpolation...etc)\nlibrary(data.table) #fast and convenient data manipulation\nlibrary(pracma) #I use it for erfc() function\nlibrary(tmap) #map library\nlibrary(cppRouting) #fast routing library\nlibrary(stars) # I use it for converting raster into polygons\nlibrary(doParallel) #multi threading\nlibrary(repr) # notebook output formating\noptions(warn=-1)\n\nlist.files(path = \"../input\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Data from sentinel TROPOMI\n\nThe TROPOspheric Monitoring Instrument (TROPOMI) is the satellite instrument on board of the Copernicus Sentinel-5 Precursor satellite. "},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"files<-list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2\")\nrs<-raster(paste0(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/\",files[1]),band=1)\n\n#read all rasters then put it in a table : one column per day\n\nread_raster_band<-function(file,band){\n  rs<-raster(paste0(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/\",file),band=band)\n  rs<-flip(rs,direction=\"y\")\n  result<-list(getValues(rs))\n  names(result)<-substr(file,9,16)\n  return(result)\n}\n#NO2\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files,fun=read_raster_band,band=1)\nstopCluster(cl)\nnam<-sapply(ok,names)\nno2<-lapply(ok,function(x) do.call(rbind,x))\nno2<-t(do.call(rbind,no2))\ncolnames(no2)<-nam\n#NaN to NA\nfor (i in 1:ncol(no2)) no2[is.nan(no2[,i]),i]<-NA\nno2<-as.data.frame(no2)\n            \n\n\n#cloud\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files,fun=read_raster_band,band=7)\nstopCluster(cl)\nnam<-sapply(ok,names)\ncloud<-lapply(ok,function(x) do.call(rbind,x))\ncloud<-t(do.call(rbind,cloud))\ncolnames(cloud)<-nam\ncloud<-cloud[,colnames(no2)]\n#NaN to NA\nfor (i in 1:ncol(cloud)) cloud[is.nan(cloud[,i]),i]<-NA\ncloud<-as.data.frame(cloud)\n            \n#Aerosol absorption index\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files,fun=read_raster_band,band=6)\nstopCluster(cl)\nnam<-sapply(ok,names)\naai<-lapply(ok,function(x) do.call(rbind,x))\naai<-t(do.call(rbind,aai))\ncolnames(aai)<-nam\naai<-aai[,colnames(no2)]\n#NaN to NA\nfor (i in 1:ncol(aai)) cloud[is.nan(aai[,i]),i]<-NA\naai<-as.data.frame(aai)\n            \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### NO2 annual mean and standard deviation"},{"metadata":{},"cell_type":"markdown","source":"TROPOMI data have to be flipped because the y axis is reversed. I have opened an issue [here](https://github.com/rspatial/raster/issues/95)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"options(repr.plot.width=15,repr.plot.height=8)\nrs2<-raster(matrix(rowMeans(no2,na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs2,main=\"NO2 annual mean in mol/m^2\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We clearly see that the biggest hotspot is at San Juan city."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rs2<-raster(matrix(apply(no2,1,sd,na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs2,main=\"NO2 annual standard deviation in mol/m^2\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Cloud fraction annual mean and standard deviation"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rs2<-raster(matrix(rowMeans(cloud,na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs2,main=\"Cloud fraction annual mean\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rs2<-raster(matrix(apply(cloud,1,sd,na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs2,main=\"Cloud fraction annual standard deviation\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Aerosol Absorption Index annual mean and standard deviation"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rs2<-raster(matrix(rowMeans(aai,na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs2,main=\"Aerosol Absorption Index annual mean\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rs2<-raster(matrix(apply(aai,1,sd,na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs2,main=\"Aerosol Absorption Index annual standard deviation\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### NO2 by month"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"yrmth<-unique(substr(colnames(no2),1,6))\npar(mfrow=c(1,1))\nfor (i in yrmth){\n    vars<-grepl(i,colnames(no2))\n    rs2<-raster(matrix(rowMeans(no2[,vars],na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\n    plot(rs2,main=paste(i))\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Electricity production\n### Power plant data provided by Kaggle."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"options(repr.plot.width=15,repr.plot.height=8)\ndf<-read.csv(\"../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv\",stringsAsFactors = F)\ndf$.geo<-gsub('\\\\{\"type\":\"Point\",\"coordinates\":\\\\[','',df$.geo)\ndf$.geo<-gsub('\\\\]}','',df$.geo)\ndf$lon<-sapply(strsplit(df$.geo,\",\"),function(x) as.numeric(x[1]))\ndf$lat<-sapply(strsplit(df$.geo,\",\"),function(x) as.numeric(x[2]))\ndf<-subset(df,select=-.geo)\npoint=st_as_sf(df,coords = c(\"lon\",\"lat\"),crs=4326)\n\n#Plot power plants \nrs2<-raster(matrix(rowMeans(no2,na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\n#plot(rs2,main=\"NO2 annual mean in mol/m^2\")\n#pointsp<-as(point,\"Spatial\")\n#plot(pointsp,add=TRUE,col=\"black\")\n\n               p<-tm_shape(rs2)+\n               tm_raster(style=\"cont\",palette=\"-RdYlBu\")+\n               tm_shape(point,xlim=c(-67.5,-65))+\n               tm_dots(col=\"primary_fuel\",size=.5)+\n               tm_layout(legend.outside=TRUE,legend.outside.position=\"right\")\n               p\n                             \n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since, NO2 emission come from combustion, we can visualize only fossil-fuel powered power plants. The dot size is varying with capacity in MW.  "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"p<-tm_shape(rs2)+\ntm_raster(style=\"cont\",title=\"NO2 mean\",palette=\"-RdYlBu\")+\ntm_shape(point[point$primary_fuel %in% c(\"Oil\",\"Coal\",\"Gas\"),],title=\"Fuel\")+\ntm_dots(col=\"primary_fuel\",size=\"capacity_mw\",title=\"Fuel\",scale=2)+\ntm_layout(legend.outside=TRUE,legend.outside.position=\"right\")\np\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Capacity by primary fuel"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"options(repr.plot.width=10,repr.plot.height=6)\ndf<-as.data.table(point)[,.(capacity=sum(capacity_mw)),.(primary_fuel)]\n               \np<-ggplot(df,aes(x=primary_fuel,y=capacity))+\n               geom_bar(stat=\"identity\")+\ntheme_bw()\np","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"According to Global Power Plant Database, oil is highly predominant in electricity production.  "},{"metadata":{},"cell_type":"markdown","source":"### Data from EIA.gov  \n\nLike other says, the generation seems to be wrong, so let's investigate [eia data](https://www.eia.gov/electricity/data/eia923/) to have a more detailed information on electricity production.  \nUnfortunately, we only have data for 2017 and 2018.  \nWe can plot electricity production by fuel type by month. "},{"metadata":{"trusted":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)]\nggplot(elm2,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\")+\n  ggtitle(\"2017 production by fuel type\")\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\nggplot(elm,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\")+\n  ggtitle(\"2018 production by fuel type\")\n        \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"According to EIA, natural gas is much more important than GPPD data. Maria storm is clearly visible on the 2017 plot.  \n  \n**Power plants fuel type :**                       \n* RFO : residual petroleum (= crude oil) \n* DFO : distilated petroleum (= diesel)                      \n* COL : coal  \n* HYC : hydroelectric  \n* NG : natural gas  \n* SUN : solar  \n* WND : wind  \n* OTH : other  \n\nIt seems that natural gas and residual petroleum are equally predominant in 2017 but RFO highly decreases in 2018, however there is some missing values especially in Aguirre power plant, which is the biggest RFO-fueled plant in Puerto Rico.  \n\nPlot the amount of energy produced by fossil-fuel powered power plant"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#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\nggplot(data=elm, aes(x=Plant.Name, y=value, fill=AER.Fuel.Type.Code)) +\n  geom_bar(stat=\"identity\")+\n  theme_bw()+\n  theme(axis.text.x = element_text(angle = 45, hjust = 1))+\n  labs(fill=\"Fuel type\",y=\"MWh\",x=\"\")+\n  ggtitle(\"2017 production by power plant\")\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\nggplot(data=elm, aes(x=Plant.Name, y=value, fill=AER.Fuel.Type.Code)) +\n  geom_bar(stat=\"identity\")+\n  theme_bw()+\n  theme(axis.text.x = element_text(angle = 45, hjust = 1))+\n  labs(fill=\"Fuel type\",y=\"MWh\",x=\"\")+\n  ggtitle(\"2018 production by power plant\")\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since emission factor has to be computed from 2018-07 to 2019-06, we have to estimate a new dataset for this time period.  For each 2019 month, we can average 2017 and 2018 production by each power plant. For 2018 missing values, we can impute 2017 ones.  \n  "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"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\nggplot(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":"## Puerto Rico characteristics\n\n### Population density\nThe biggest sources of NOx emissions come from traffic, industry and energy so let's investigate the population density of Puerto Rico.   \nPopulation density from 2020 come from [this paper](https://www.nature.com/articles/sdata201545#MOESM111).  \n  \nPlot population density and fossil-fuel powered power plants.  "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"options(repr.plot.width=15,repr.plot.height=8)\npop<-raster(\"../input/population-density-in-puerti-rico/PRI_ppp_v2b_2020.tif\")\n#cut borders of the raster                      \next<-extent(pop)\next[1]<--67.5\next[2]<--65.5\npop<-crop(pop,ext)\n                      \np<-tm_shape(pop)+\ntm_raster(style=\"cont\",title=\"Population\")+\ntm_shape(point[point$primary_fuel %in% c(\"Oil\",\"Coal\",\"Gas\"),],title=\"Fuel\")+\ntm_dots(col=\"primary_fuel\",size=\"capacity_mw\",title=\"Fuel\",scale=2,palette=c(\"blue\",\"green\",\"grey\"))\np ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Industry location\n                      \nAccording to [Wikipedia](https://en.wikipedia.org/wiki/Economy_of_Puerto_Rico), Puerto Rico's industry is dominated by pharmaceutical and petrochemical industries.  \nI have extracted industrial areas from OpenStreetMap (fclass=industrial).  \n  \nPlot industries (red dots) and fossil-fuel powered power plants."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"options(repr.plot.width=15,repr.plot.height=8)\nindustry<-read_sf(\"../input/various-data-about-puerto-rico/puerto_rico_osm_landuse.shp\")\nindustry<-st_centroid(industry)\n#read shapefile of county\nshp<-read_sf(\"../input/various-data-about-puerto-rico/puerto_rico_county.shp\")\nshp<-st_transform(shp,crs=4326)\nshp<-shp[shp$ALAND!=0,]\nshp$county<-paste0(shp$STATEFP,shp$COUNTYFP)\nshp<-st_crop(shp,ext)\nshp2<-summarise(group_by(shp,county))\n\n                      \n  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"p<-tm_shape(shp2)+\ntm_polygons(col=\"white\",title=\"Counties\")+\ntm_shape(industry[industry$fclass %in% c(\"industrial\"),])+\ntm_dots(col=\"red\",size=.3)+\ntm_shape(point[point$primary_fuel %in% c(\"Oil\",\"Coal\",\"Gas\"),],title=\"Fuel\")+\ntm_dots(col=\"primary_fuel\",size=\"capacity_mw\",title=\"Fuel\",scale=2,palette=c(\"blue\",\"green\",\"grey\"))+\ntm_layout(legend.outside=FALSE)\np ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Estimate traffic from commuting data  \n \nCommuting data coming from [2015 census](https://www.census.gov/data/tables/2015/demo/metro-micro/commuting-flows-2015.html).\nWe have county pairs (residence county - workplace county) and the flow of commuters for each pair.  We can use Openstreetmap data to construct a routable network and look at the most frequented counties. A graph is composed of edges (= road) connected by nodes (= road intersection)\n  \nHead of commuting data :  \n                  "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#import data and cleaning\ncom<-read.csv(\"../input/various-data-about-puerto-rico/commuting_puerto_rico_2015.csv\",colClasses=c(\"character\",\"character\",\"NULL\",\"NULL\",\"character\",\"character\",\"NULL\",\"NULL\",\"numeric\",\"numeric\",\"NULL\",\"NULL\",\"NULL\",\"NULL\"))\ncom<-com[!is.na(com$County.FIPS.Code) & !is.na(com$County.FIPS.Code.1),]\ncom$County.FIPS.Code<-sapply(com$County.FIPS.Code,function(x){\n  nzero<-nchar(x)\n  return(paste0(paste(rep(0,3-nzero),collapse = \"\"),x))\n})\n\ncom$County.FIPS.Code.1<-sapply(com$County.FIPS.Code.1,function(x){\n  nzero<-nchar(x)\n  return(paste0(paste(rep(0,3-nzero),collapse = \"\"),x))\n})\n\nhead(com)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Road shapefile come from [OSM](http://download.geofabrik.de/north-america/us.html). I have extracted data with [osm2po](https://osm2po.de/) tool."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"cent<-st_centroid(st_transform(shp2,3991))\ngraph<-read_sf(\"../input/various-data-about-puerto-rico/osm2po_puerto_rico.shp\")\n#Plot road network\nsuppressMessages(graph<-st_crop(graph,ext))\nplot(st_geometry(graph),main=\"Openstreetmap roads\")\n\n\n#Extract nodes\ncoo<-data.frame(id=c(graph$osm_source,graph$osm_target),\n                X=c(graph$x1,graph$x2),\n                Y=c(graph$y1,graph$y2))\ncoo<-coo[!duplicated(coo$id),]\n\n#Spatial join network nodes to county centroids buffer to avoid creating huge distance matrix\ncoo<-st_as_sf(coo,coords = c(\"X\",\"Y\"),crs=4326)\ncoo<-st_transform(coo,3991)\ncoo2<-st_join(coo,st_buffer(cent,4000))\ncoo2<-coo2[!is.na(coo2$county),]\ncent<-cent[!cent$county %in% setdiff(shp2$county,coo2$county),]\nd<-pointDistance(st_coordinates(coo2),st_coordinates(cent),allpairs = TRUE,lonlat = FALSE)\nd<-apply(t(d),1,which.min)\n\ncounty_nodes<-data.frame(county=cent$county,\n                         node=as.character(coo2$id[d]),stringsAsFactors = FALSE)\n#Graph preparation\nst_geometry(graph)<-NULL\ngraph$km<-as.numeric(graph$km)\ngraph$cost<-as.numeric(graph$cost)\ngraph$reverse_co<-as.numeric(graph$reverse_co)\ngraph<-data.frame(graph)\ndat<-graph[,c(\"osm_source\",\"osm_target\",\"km\",\"cost\")]\ndat2<-graph[graph$reverse_co<100000,c(\"osm_source\",\"osm_target\",\"km\",\"reverse_co\")]\ncolnames(dat2)<-c(\"osm_target\",\"osm_source\",\"km\",\"cost\")\ndat<-rbind(dat,dat2)\ndat<-dat[dat$osm_source!=dat$osm_target,]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Head of routable graph (each row is an edge between two nodes)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"head(dat)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Compute shortest paths for each county pairs and plot the flow of people/commuters crossing each county \n#### Quantile scale"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gr<-makegraph(dat[,c(1,2,4)])\n\n\n#Merge network nodes to counties\ncom$from<-paste0(com[,1],com[,2])\ncom$to<-paste0(com[,3],com[,4])\ncom<-merge(com,county_nodes,by.x=\"from\",by.y=\"county\")\ncolnames(com)[9]<-\"from_node\"\ncom<-merge(com,county_nodes,by.x=\"to\",by.y=\"county\")\ncolnames(com)[10]<-\"to_node\"\n\n#contract cppRouting graph for faster computing\ngr2<-cpp_contract(gr,silent = TRUE)\n#shortest paths\nt1<-get_path_pair(gr2,com$from_node,com$to_node,long = T)\n#Spatial join of network nodes to county polygons\nshp<-st_transform(shp,3991)\ncoun<-st_join(shp,coo)\nst_geometry(coun)<-NULL\ncoun$id<-as.numeric(coun$id)\n#Merging and agregating commuters flow by county\nt1<-merge(t1,coun[,c(\"id\",\"COUSUBFP\")],by.x=\"node\",by.y=\"id\",all.x=TRUE)\nt1<-t1[,-1]\nt1<-t1[!duplicated(t1),]\ncom$Workers.in.Commuting.Flow<-as.numeric(gsub(\" \",\"\",com$Workers.in.Commuting.Flow))\nt1<-merge(t1,com[,c(\"from_node\",\"to_node\",\"Workers.in.Commuting.Flow\")],by.x=c(\"from\",\"to\"),by.y=c(\"from_node\",\"to_node\"),all.x=TRUE)\nflow<-setDT(t1)[,.(flow=sum(Workers.in.Commuting.Flow,na.rm=T)),.(COUSUBFP)]\n\nshp<-merge(shp,flow,by=\"COUSUBFP\",all.x=T)\nshp$flow[is.na(shp$flow)]<-0\n#plot\np<-tm_shape(shp)+\n  tm_fill(col=\"flow\",style=\"quantile\",palette=\"Reds\",n=8,title=\"Commuting flow\")+\n  tm_shape(point[point$primary_fuel %in% c(\"Oil\",\"Coal\",\"Gas\"),],title=\"Fuel\")+\n  tm_dots(col=\"primary_fuel\",size=\"capacity_mw\",title=\"Fuel\",scale=2,palette=c(\"blue\",\"green\",\"grey\"))+\n  tm_layout(legend.outside=FALSE,legend.position=c(\"left\",\"bottom\"))\np ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Continuous scale"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"p<-tm_shape(shp)+\n  tm_fill(col=\"flow\",style=\"cont\",palette=\"Reds\",title=\"Commuting flow\")+\n  tm_shape(point[point$primary_fuel %in% c(\"Oil\",\"Coal\",\"Gas\"),],title=\"Fuel\")+\n  tm_dots(col=\"primary_fuel\",size=\"capacity_mw\",title=\"Fuel\",scale=2,palette=c(\"blue\",\"green\",\"grey\"))+\n  tm_layout(legend.outside=FALSE,legend.position=c(\"left\",\"bottom\"))\np ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### San Juan\n\nWe can compare NO2 variations in San Juan and the rest of Puerto Rico.  \nThe hotspot pixels are selected by taking the 99th percentile of averaged NO2 distribution. "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rs2<-raster(matrix(rowMeans(no2,na.rm=T),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nlim<-quantile(c(rowMeans(no2,na.rm=T)),\n             probs=seq(0.01,1,0.01))   \nlim<-lim[99]                      \nrs3<-rs2*0\nrs3[rs2>lim]<-1\nplot(rs3,main=\"Hotspot\")\nhs_pix<-which(getValues(rs3)==1)\n\ntemp<-as.data.table(no2)\ntemp[,pixel:=1:ncell(rs),]\ntemp<-melt(temp,id.vars = \"pixel\",variable.factor = TRUE)\ntemp2<-temp[pixel %in% hs_pix,.(value=mean(value,na.rm=T)),.(variable)]\ntemp2$class<-\"hotspot\"\ntemp3<-temp[,.(value=mean(value,na.rm=T)),.(variable)]\ntemp3$class<-\"all\"\ntemp4<-temp[!pixel %in% hs_pix,.(value=mean(value,na.rm=T)),.(variable)]\ntemp4$class<-\"without hotspot\"\ntemp2<-rbind(temp2,temp3,temp4)\n\n\nlabels<-unique(as.character(temp2$variable))\nlabels<-gsub(\"_\",\"\",labels)\nbreaks<-which(substr(labels,7,8)==\"01\")\n\noptions(repr.plot.width=10,repr.plot.height=6)\np<-ggplot(data=temp2,aes(x=as.numeric(variable),y=value,colour=class))+\n  geom_line()+\n  scale_x_continuous(breaks=breaks,labels=labels[breaks])+\n  theme_bw()+\n  theme(axis.text.x = element_text(angle = 45, hjust = 1))+\n  labs(colour=\"\",x=\"\",y=\"NO2 density\")+\n  ggtitle(\"Daily averaged NO2 density\")\n\np     ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Week-end effect \n\nComparison between weekdays and week-end (saturday and sunday).  "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"semaine<-as.data.table(no2)\nsemaine[,pixel:=1:ncell(rs),]                      \nsemaine<-melt(semaine,id.vars = c(\"pixel\"),variable.factor = FALSE)\nsemaine$month<-substr(semaine$variable,5,6)\nsemaine$variable<-gsub(\"_\",\"\",semaine$variable)\n#day corresponding to the date\nweek_days<-data.frame(date=unique(semaine$variable),\n                      day=weekdays(as.Date(unique(semaine$variable),format=\"%Y%m%d\")))\n\nsemaine$weekday<-week_days$day[match(semaine$variable,week_days$date)]\nsemaine$weekend<-ifelse(semaine$weekday %in% c(\"Saturday\",\"Sunday\"),1,0)\nsemaine$hotspot<-ifelse(semaine$pixel %in% hs_pix,1,0)\nsemaine2<-semaine[,.(no2=mean(value,na.rm=T)),.(month,weekend,hotspot)]\nindex<-data.frame(month=unique(semaine2$month),ind=1:12)\nsemaine2$index<-index$ind[match(semaine2$month,index$month)]\n\np<-ggplot(data=semaine2,aes(x=index,y=no2,colour=factor(weekend),linetype=factor(hotspot)))+\n  geom_line()+\n  scale_x_continuous(breaks=1:12,labels=index$month)+\n  theme_bw()+\n  labs(linetype=\"Hotspot\",colour=\"Week-end\",x=\"Months\",y=\"NO2 density\")+\n  ggtitle(\"Week-end effect on NO2 density\")\n\np","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Wind data  \n                      \nI choose to use [ERA interim](https://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=ml/) data instead of GFS because of the much higher resolution, plus ERA interim have datasets at different altitude. It seems better to have a higher altitude than 10m because we have a more representative view of the tropospheric behavior. ERA interim data are daily data and consist of u and v components at about 500m of altitude (level 53 of the dataset). I have selected 18 UTC to be close to NO2 data (between 15h and 18h UTC). \nWind speed is simply sqrt(u²+v²), expressed in m/s. \nWind direction is expressed in degrees from North (e.g 0° = from south to north, -90° = from East to West).  \nWind data have been resampled by the nearest neighbor method to match TROPOMI extent and resolution.  \n\n"},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"\n#First, wind data need to be resampled to NO2 to match the extent and resolution\nfiles2<-list.files(\"../input/era-interim-wind-data-in-puerto-rico/era_interim\")\n#create output folder\ndir.create(\"../working/era_resampled\")\n\n\nresample_era<-function(file,model){\n  era<-raster(paste0(\"../input/era-interim-wind-data-in-puerto-rico/era_interim/\",file))\n  n<-nbands(era)\n  for (j in c(\"u\",\"v\")){\n    for (k in 1:n){\n      rs<-raster(paste0(\"../input/era-interim-wind-data-in-puerto-rico/era_interim/\",file),band=k,varname=paste(j))\n      day<-as.character(k)\n      if (nchar(day)==1) day<-paste0(\"0\",day)\n      date<-paste0(substr(file,14,19),day)\n      \n      rs<-resample(rs,model,method=\"ngb\")\n      \n      \n      writeRaster(rs,paste0(\"../working/era_resampled/\",paste0(j,\"_\",date,\".tif\")),overwrite=TRUE)\n    }\n    \n    \n  }\n}\n\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files2,fun=resample_era,model=rs2)\nstopCluster(cl)\n\n##Read data\n#u component\nfiles2<-list.files(\"../working/era_resampled\")\nfiles2<-files2[grepl(\"u_\",files2)]\nread_raster<-function(file){\n  rs<-raster(paste0(\"../working/era_resampled/\",file),band=1)\n  result<-list(getValues(rs))\n  names(result)<-substr(file,3,10)\n  return(result)\n}\n\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files2,fun=read_raster)\nstopCluster(cl)\nnam<-sapply(ok,names)\ndfu<-lapply(ok,function(x) do.call(rbind,x))\ndfu<-t(do.call(rbind,dfu))\ncolnames(dfu)<-nam\ndfu<-dfu[,colnames(no2)]#keep the same column order between all tables            \n            \n#read v component\nfiles2<-list.files(\"../working/era_resampled\")\nfiles2<-files2[grepl(\"v_\",files2)]\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files2,fun=read_raster)\nstopCluster(cl)\nnam<-sapply(ok,names)\ndfv<-lapply(ok,function(x) do.call(rbind,x))\ndfv<-t(do.call(rbind,dfv))\ncolnames(dfv)<-nam\ndfv<-dfv[,colnames(dfu)]\n            \n#Wind speed\ndfw<-sqrt(dfu^2+dfv^2)\ncolnames(dfw)<-colnames(dfu)\n            \n#mean direction\nwind_dir<-matrix(atan2(c(dfu/dfw),c(dfv/dfw))*180/pi,ncol=ncol(dfu),nrow=nrow(dfu))\ncolnames(wind_dir)<-colnames(dfu)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Maps of averaged wind speed and averaged direction"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"options(repr.plot.width=15,repr.plot.height=8)\nrs5<-raster(matrix(rowMeans(wind_dir,na.rm=T),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs5,main=\"Averaged wind direction\")\nrs5<-raster(matrix(apply(wind_dir,1,sd,na.rm=T),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs5,main=\"Wind direction annual standard deviation\")\n\nrs4<-raster(matrix(rowMeans(dfw,na.rm=T),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs4,main=\"Averaged wind speed\")\nrs4<-raster(matrix(apply(dfw,1,sd,na.rm=T),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs4,main=\"Wind speed annual standard deviation\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Wind is particularly constant in speed and direction (about -70°)\n"},{"metadata":{},"cell_type":"markdown","source":"# Methodology\n\nI choose to focus exclusively on Aguirre power plant (RFO fueled in the south) because :\n* it is far from large cities\n* sparsely populated compared to other PP\n* not many industries\n* not many traffic\n* visible from TROPOMI data (I guess)\n\nWe can try to fit an Exponentially Modified Gaussian (EMG) function to its plume for estimating the NO2 quantity emission. For the complete methodology, please read this following papers :\n* [Enhanced Capabilities of TROPOMI NO2: Estimating NOX from North\nAmerican Cities and Power Plants](https://pubs.acs.org/doi/pdf/10.1021/acs.est.9b04488)\n* [Direct observation of changing NOx lifetime in North American cities](https://authors.library.caltech.edu/99738/2/aax6832-Laughner-SM.pdf)\n\nBasically, the method can be separated into three steps :\n\n1. rotate NO2 data depending the daily wind direction, then take the average of all the rotated days. The NO2 hotspot is called a plume.\n2. the plume is spatially integrated along the x axis (i.e wind direction). See [page 17 of this paper](https://kfolkertboersma.files.wordpress.com/2019/05/quantifying_nox_emissions_from_paris.pdf) for detailed methodology.\n3. fit an EMG function to the line density previously created\n\nThe NO2 emission (in mol/hour) is derived from the estimated parameters of the EMG function. The selected NO2 data for the first step has to be the less cloudy as possible.\n\nZoom on Aguirre power plant.\n"},{"metadata":{"_kg_hide-output":false,"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rs2<-raster(matrix(rowMeans(no2,na.rm=TRUE),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nrs_crop<-crop(rs2,extent(c(-66.5,-66.15,17.9,18.1)))\np<-tm_shape(rs_crop)+\n   tm_raster(style=\"cont\",palette=\"-RdYlBu\")+\n   tm_shape(point[point$gppd_idnr==\"WRI1026809\",])+\n   tm_dots(col=\"primary_fuel\",size=2)+\n   tm_layout(legend.outside=TRUE,legend.outside.position=\"right\")\n   p\n                  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"#Function to crop and rotate rasters\nfun2<-function(pp,#power plant id\n               no2,#no2 data frame\n               cloud_lim,#cloud fraction limit\n               cloud_cover,#rate of the area with cloud fraction > cloud_lim\n               radius_stat, #buffer around the power plant, in which cloud cover, averaged wind speed and averaged wind direction are computed\n               radius_rot, #rotated area radius\n               cloud,# cloud data frame\n               wind_dir,#wind direction data frame\n               windspeed,#wind speed data frame\n               point2,#shapefile of power plant location\n               raster_model){ #a template raster for the output\n  ##buffer\n  #reproject in a plane CRS for meter unit\n  point_2154<-st_transform(point2,2154)\n  #radius rotation buffer\n  point2<-st_buffer(point_2154,dist=radius_rot)\n  point2<-st_transform(point2,4326)\n  #radius stat buffer\n  point3<-st_buffer(point_2154,dist=radius_stat)\n  point3<-st_transform(point3,4326)\n  \n  #Get pixel IDs of the two buffers\n  pix<-rasterize(point2[point2$gppd_idnr==pp,],raster_model)\n  pix<-getValues(pix)\n  pix<-which(!is.na(pix))\n  pix2<-rasterize(point3[point3$gppd_idnr==pp,],raster_model)\n  pix2<-getValues(pix2)\n  pix2<-which(!is.na(pix2))\n  \n  #select the day where cloud condition is respected\n  vars<-which((colSums(cloud[pix2,] > cloud_lim,na.rm=F)/length(pix2)) < cloud_cover)\n  #averaged direction\n  meandir<-colMeans(wind_dir[pix2,vars])[!is.na(colMeans(no2[pix,vars]))]\n  #averaged speed\n  meanspeed<-colMeans(windspeed[pix2,vars])[!is.na(colMeans(no2[pix,vars]))]\n  \n  #output raster \n  grid<-raster(extent(point2[point2$gppd_idnr==pp,]),resolution=res(raster_model))\n  result=list()\n  \n  for (i in vars){\n    #raster from no2 data.frame\n    rs<-raster(matrix(no2[,i],ncol=ncol(raster_model),nrow=nrow(raster_model),byrow = TRUE),template=raster_model)\n    df<-no2[pix,i]\n    #if missing value, next\n    if (any(is.na(df))) next\n    #averaged wind direction\n    direction<-mean(wind_dir[pix,i])\n    #crop no2 raster\n    cropped<-crop(rs,point2[point2$gppd_idnr==pp,])\n    #raster to polygons\n    pol<-st_as_sf(st_as_stars(cropped))\n    ge<-st_coordinates(st_centroid(pol)) \n    ax<-st_coordinates(st_centroid(point2[point2$gppd_idnr==pp,]))\n    rotang=90-direction\n    #rotang=45\n    rot = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)\n    #Rotation to W-E direction\n    x=(ge[,1]-ax[,1])*cos(-rotang*pi/180)-(ge[,2]-ax[,2])*sin(-rotang*pi/180)+ax[,1]\n    y=(ge[,1]-ax[,1])*sin(-rotang*pi/180)+(ge[,2]-ax[,2])*cos(-rotang*pi/180)+ax[,2]\n    pol2<-st_as_sf(data.frame(x,y),coords=c(\"x\",\"y\"),crs=4326)\n    pol2$value<-pol[[1]]  \n    \n    rst<-rasterize(pol2,grid,field=\"value\",na.rm=T)\n    rst<-resample(rst,rst)\n    result[[length(result)+1]]<-rst\n  }\n  \n  final<-list()\n  final[[1]]<-result\n  final[[2]]<-meandir\n  final[[3]]<-meanspeed\n  \n  return(final)\n  \n}\n\n#Function to agregate rotated rasters\nag_rasters<-function(obj,#resulting object of the previous function\n                     minspeed=-Inf,\n                     maxspeed=Inf, #wind speed criterions\n                     minangle=-360,\n                     maxangle=360, #wind direction criterions\n                     na=FALSE){\n  mod<-obj[[1]][[1]]\n  df<-do.call(cbind,lapply(obj[[1]],getValues))\n  vars<-which(obj[[2]] > minangle & obj[[2]] < maxangle & obj[[3]] < maxspeed & obj[[3]] > minspeed)\n  res<-raster(matrix(rowMeans(as.matrix(df[,vars]),na.rm=na),nrow=nrow(mod),ncol=ncol(mod),byrow = T),template=mod)\n  return(res)\n}\n\n#Function to compute line density\n\nline_dens<-function(raster_result,#resulting object of the previous function\n                    largeur, #size of the rectangle for calculating densities\n                    fill=NULL, #NA values to fill with (in the borders)\n                    pts, #power plant shapefile\n                    pp#power plant ID\n                    \n){\n  if (!is.null(fill)) raster_result[is.na(raster_result)]<-fill\n  crs(raster_result)<-crs(pts)\n  #reproject raster to a plane projection\n  rs2154<-projectRaster(raster_result,crs=\"+proj=lcc +lat_1=18.43333333333333 +lat_2=18.03333333333333 +lat_0=17.83333333333333 +lon_0=-66.43333333333334 +x_0=200000 +y_0=200000 +ellps=GRS80 +units=m +no_defs\")\n  #raster to polygons\n  pol<-rasterToPolygons(rs2154,dissolve = FALSE,na.rm=TRUE)\n  pol<-st_as_sf(pol)\n  pol<-st_cast(pol,\"POLYGON\")\n  ext<-extent(pol)\n    \n  #Create vertical grid shapefile \n  x0<-st_coordinates(st_transform(pts[pts$gppd_idnr==pp,],crs=crs(pol)))\n  mat<-data.frame(x=rep(seq(ext[1],ext[2],res(rs2154)[1]),each=2),\n                  y=rep(c((x0[2]-largeur),(x0[2]+largeur)),length(seq(ext[1],ext[2],res(rs2154)[1]))))\n  mat$x<-mat$x+(res(rs2154)[1]/2)\n  mat2<-split(mat,f=rep(1:(nrow(mat)/2),each=2))\n  mat2<-lapply(mat2,function(x) st_linestring(as.matrix(x),dim=\"XY\"))\n  mat2<-st_sfc(mat2)\n  mat2<-st_sf('id'=unique(mat$x),'geometry'=mat2,crs=crs(pol))\n               \n  #crop the rotated plume with the rectangle for visual inspection             \n  cropped<-crop(rs2154,mat2)\n  \n  #intersection between vertical lines and the plume\n  inter<-st_intersection(mat2,pol)\n  inter$len<-as.numeric(st_length(inter))\n  #Remove lines in the borders\n  naid<-setDT(inter)[,.(len=sum(len)),.(id)]\n  #Calculate density\n  inter<-setDT(inter)[!id %in% naid$id[naid$len!=2*largeur],.(density=sum(len*layer)),.(id)] \n  inter<-inter[order(inter$id),]\n\n  return(list(rs=cropped,\n             inter=inter))\n}\n               \n#Function to plot line density\nplot_ld<-function(line_density,fit=NULL){\n    p<-ggplot(line_density,aes(x=id,y=density))+\n    geom_point()+\n    labs(x=\"meters\",y=\"mol/m\")+\n    theme_bw()\n    \n    if (!is.null(fit)) p<-p+geom_line(data=data.frame(id=line_density$id,density=fit),aes(x=id,y=density),colour=\"red\")\n    return(p)\n}\n               \n#Function to compute NO2 emissions after fitting\ncompute_no2<-function(parameters,\n                      fun2_output,\n                     minangle=-Inf,\n                     maxangle=Inf,\n                     minspeed=-360,\n                     maxspeed=360){\ntau<-parameters[2]/mean(fun2_output[[3]][fun2_output[[2]] > minangle & fun2_output[[2]] < maxangle & fun2_output[[3]] < maxspeed & fun2_output[[3]] > minspeed])\nno2_em=parameters[1]/tau\n#convert to ton/h\nreturn(no2_em*46.055*3600/1000000)\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Aguirre NO2 emission\nLet's try in Aguirre power plant."},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"test<-fun2(pp=\"WRI1026809\",\n          no2 = data.frame(no2),\n          cloud_lim = 0.5,\n          cloud_cover = 0.01,\n          radius_stat = 20000,\n          radius_rot=45000,\n          cloud=data.frame(cloud),\n          wind_dir=data.frame(wind_dir),\n          windspeed = data.frame(dfw),\n          point2=point,\n          raster_model = rs)\n               \nrstest<-ag_rasters(test)\nplot(rstest,main=paste(length(test[[1]]),\"rotated days\"))\nplot(point[point$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")   ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Infortunately, Aguirre is near the border of the tif so with rotations we lose informations.\nLet's try with angle restriction"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rstest<-ag_rasters(test,minangle=-90)\nplot(rstest,main=paste(length(test[[1]][test[[2]] > -90]),\"rotated days - wind direction > -90\"))\nplot(point[point$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")\n               \nrstest<-ag_rasters(test,minangle=-75)\nplot(rstest,main=paste(length(test[[1]][test[[2]] > -75]),\"rotated days - wind direction > -75\"))\nplot(point[point$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")       ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Well, the plume seems to look good but we cannot use all of the data. Plus, a big part of NO2 background (i.e around the plume) is missing.\nThat's why I have downloaded the same data from Earth Engine in an extended area.\nLet's see it !\n"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"#Read extended NO2 data\nfiles<-list.files(\"../input/s5p-no2-tropomi-in-puerto-rico\")\nold_no2_files<-list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2\")\nfiles=files[files %in% substr(old_no2_files,9,nchar(old_no2_files))]\nprint(length(files))\n#My template, if for some unknown reasons all the files are not aligned\nrs_model<-raster(paste0(\"../input/s5p-no2-tropomi-in-puerto-rico/\",files[1]))\n\nread_raster_no2<-function(file,model){\n  rs<-raster(paste0(\"../input/s5p-no2-tropomi-in-puerto-rico/\",file),band=1)\n  rs<-flip(rs,direction=\"y\")\n  if (ncell(rs) != ncell(model)) rs<-resample(rs,model,method=\"ngb\")\n  result<-list(getValues(rs))\n  names(result)<-substr(file,1,8)\n  return(result)\n}\n\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files,fun=read_raster_no2,model=rs_model)\nstopCluster(cl)\nnam<-sapply(ok,names)\nno2<-lapply(ok,function(x) do.call(rbind,x))\nno2<-t(do.call(rbind,no2))\ncolnames(no2)<-nam\n\n\n\n#Read cloud data\nread_raster_cloud<-function(file,model){\n  rs<-raster(paste0(\"../input/s5p-no2-tropomi-in-puerto-rico/\",file),band=7)\n  rs<-flip(rs,direction=\"y\")\n  if (ncell(rs) != ncell(model)) rs<-resample(rs,model,method=\"ngb\")\n  result<-list(getValues(rs))\n  names(result)<-substr(file,1,8)\n  return(result)\n}\n            \ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files,fun=read_raster_cloud,model=rs_model)\nstopCluster(cl)\nnam<-sapply(ok,names)\ncloud<-lapply(ok,function(x) do.call(rbind,x))\ncloud<-t(do.call(rbind,cloud))\ncolnames(cloud)<-nam\ncloud<-cloud[,colnames(no2)]\n              \n\n              \n#First, wind data need to be resampled to NO2 to match the extent and resolution\nfiles2<-list.files(\"../input/era-interim-wind-data-in-puerto-rico/era_interim\")\n#create output folder\ndir.create(\"../working/ext_era_resampled\")\n\n\nresample_era<-function(file,model){\n  era<-raster(paste0(\"../input/era-interim-wind-data-in-puerto-rico/era_interim/\",file))\n  n<-nbands(era)\n  for (j in c(\"u\",\"v\")){\n    for (k in 1:n){\n      rs<-raster(paste0(\"../input/era-interim-wind-data-in-puerto-rico/era_interim/\",file),band=k,varname=paste(j))\n      day<-as.character(k)\n      if (nchar(day)==1) day<-paste0(\"0\",day)\n      date<-paste0(substr(file,14,19),day)\n      \n      rs<-resample(rs,model,method=\"ngb\")\n      \n      \n      writeRaster(rs,paste0(\"../working/ext_era_resampled/\",paste0(j,\"_\",date,\".tif\")),overwrite=TRUE)\n    }\n    \n    \n  }\n}\n\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files2,fun=resample_era,model=rs_model)\nstopCluster(cl)\n\n##Read data\n#u component\nfiles2<-list.files(\"../working/ext_era_resampled\")\nfiles2<-files2[grepl(\"u_\",files2)]\nread_raster<-function(file){\n  rs<-raster(paste0(\"../working/ext_era_resampled/\",file),band=1)\n  result<-list(getValues(rs))\n  names(result)<-substr(file,3,10)\n  return(result)\n}\n\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files2,fun=read_raster)\nstopCluster(cl)\nnam<-sapply(ok,names)\ndfu<-lapply(ok,function(x) do.call(rbind,x))\ndfu<-t(do.call(rbind,dfu))\ncolnames(dfu)<-nam\ndfu<-dfu[,colnames(no2)]            \n            \n#read v component\nfiles2<-list.files(\"../working/ext_era_resampled\")\nfiles2<-files2[grepl(\"v_\",files2)]\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files2,fun=read_raster)\nstopCluster(cl)\nnam<-sapply(ok,names)\ndfv<-lapply(ok,function(x) do.call(rbind,x))\ndfv<-t(do.call(rbind,dfv))\ncolnames(dfv)<-nam\ndfv<-dfv[,colnames(no2)]\n            \n#Wind speed\ndfw<-sqrt(dfu^2+dfv^2)\ncolnames(dfw)<-colnames(dfu)\n            \n#mean direction\nwind_dir<-matrix(atan2(c(dfu/dfw),c(dfv/dfw))*180/pi,ncol=ncol(dfu),nrow=nrow(dfu))\ncolnames(wind_dir)<-colnames(dfu)  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#Plot\nrs<-raster(matrix(rowMeans(no2,na.rm=TRUE),ncol=ncol(rs_model),nrow=nrow(rs_model),byrow=TRUE),template=rs_model)\nplot(rs,main=\"NO2 annual mean extended version\")\n\nrs<-raster(matrix(rowMeans(cloud,na.rm=TRUE),ncol=ncol(rs_model),nrow=nrow(rs_model),byrow=TRUE),template=rs_model)\nplot(rs,main=\"Cloud fraction annual mean extended version\")  \n\nrs4<-raster(matrix(rowMeans(dfw,na.rm=T),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs4,main=\"Averaged wind speed - extended\")\n\nrs5<-raster(matrix(rowMeans(wind_dir,na.rm=T),ncol=ncol(rs),nrow=nrow(rs),byrow = T),template=rs)\nplot(rs5,main=\"Averaged mean direction - extended\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The function for raster rotation take many arguments :\n* considered power plant\n* max cloud fraction allowed\n* rate of the considered area exceeding max cloud fraction\n* radius of rotation area\n* radius of statistic area (where cloud fraction criterion is applied)\n\nHere are the averaged rotated days centered on Aguirre power plant with the following parameters : \n* cloud fraction criterion = 0.5\n* rate of statistic area exceeding cloud fraction criterion : 0.05\n* radius of rotation area : 45km\n* radius of statistic area : 20km"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"test<-fun2(pp=\"WRI1026809\",\n          no2 = data.frame(no2),\n          cloud_lim = 0.5,\n          cloud_cover = 0.05,\n          radius_stat = 20000,\n          radius_rot=45000,\n          cloud=data.frame(cloud),\n          wind_dir=data.frame(wind_dir),\n          windspeed = data.frame(dfw),\n          point2=point,\n          raster_model = rs)\nrstest<-ag_rasters(test)\nplot(rstest,main=paste(length(test[[1]]),\"rotated days\"))\nplot(point[point$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's try different angle restriction to have the best neat plume possible"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"par(mfrow=c(2,3))\nfor (i in seq(-90,-40,10)) {\n  rstest<-ag_rasters(test,minangle = i)\n  plot(rstest,main=paste(i))\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Fitting EMG with wind direction > -70 seems to be the best choice.  \nWe can spatially integrate the plume within a rectangle of 30km by ~90km (2*rotation radius)"},{"metadata":{"trusted":true,"collapsed":true,"_kg_hide-input":true},"cell_type":"code","source":"rstest<-ag_rasters(test,minangle=-70)\nlide<-line_dens(raster_result =  rstest,\n                largeur = 15000, #2 * 15km\n                fill=NULL,\n                pts=point,\n                pp=\"WRI1026809\")\n\n#Plot the plume\nplot(lide$rs,main=\"NO2 plume in mol/m^2\")\npoint2<-st_transform(point,crs=crs(lide$rs))\nplot(point2[point2$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")\n\n\n#cost function (RMSE)\nfit_gaussian<-function(x,dat,res){\n  test<-x[5]+(x[1]/(2*x[2]))*exp((x[3]/x[2])+((x[4]^2)/(2*(x[2]^2)))-(dat/x[2]))*erfc((-1/sqrt(2))*(((dat-x[3])/x[4])-(x[4]/x[2])))\n  return(sum((test-res)^2))\n}\n#predict function\npred_gaussian<-function(x,dat){\n  test<-x[5]+(x[1]/(2*x[2]))*exp((x[3]/x[2])+((x[4]^2)/(2*(x[2]^2)))-(dat/x[2]))*erfc((-1/sqrt(2))*(((dat-x[3])/x[4])-(x[4]/x[2])))\n  return(test)\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And finally plot line density"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#Plot line density\np<-ggplot(lide$inter,aes(x=id,y=density))+\ngeom_point()+\nlabs(x=\"meters\",y=\"mol/m\")+\ntheme_bw()\np\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Fitting parameters are very sensitive here, we have 5 parameters to estimate :\n\n* no2 quantity in the plume : the integral within the plume in mol\n* the e-folding distance : distance with which the no2 density decay by e(1)\n* the location of the plume center (x)\n* the standard deviation of the Gaussian function\n* no2 background in mol/m\n\nIf there are NO2 contamination before the plume, NO2 background will be higher and NO2 emission from the source will be under-estimated.\nSo there must be a visual inspection before fitting for possibly removing potential contamination in the curve.  \nSo, let's fit with stats::optim function.\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#fitting EMG with stats::optim function\nlide2<-lide$inter\n#lide2<-lide2[lide2$id < 255000,]\n\n\npar<-optim(par=c(7000,20000,235000,.8,1.55),\n           fn=fit_gaussian,\n           dat=lide2$id,\n           res=lide2$density,\n           control = list(maxit=10000)\n)\n\n#predict\nYfit<-pred_gaussian(par$par,lide$inter$id)\n#plot line densities and fitted EMG\nplot_ld(lide$inter,Yfit)\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#NO2 emission\nprint(paste(\"Annual NO2 emission for Aguirre plant is\",round(compute_no2(parameters=par$par,fun2_output=test,minangle= -70),4),\"ton/h\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's try for March where the plume is the cleanest in the dataset. (see EDA part)\n\n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"vars<-substr(colnames(no2),1,6) %in% c(\"201903\")\n  test_03<-fun2(pp=\"WRI1026809\",\n             no2 = no2[,vars],\n             cloud_lim = 0.5,\n             cloud_cover = 0.05,\n             radius_stat = 20000,\n             radius_rot=65000,\n             cloud=cloud[,vars],\n             wind_dir=wind_dir[,vars],\n             windspeed = dfw[,vars],\n             point2=point,\n             raster_model = rs)\nrstest<-ag_rasters(test_03)\nplot(rstest,main=paste(length(test_03[[1]]),\"rotated days\"))\nplot(point[point$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The NO2 plume looks good enough.  \nSpatial integration + line density :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rstest<-ag_rasters(test_03)\nlide<-line_dens(raster_result =  rstest,\n                largeur = 15000, #2 * 10km\n                fill=NULL,\n                pts=point,\n                pp=\"WRI1026809\")\n\n#Plot the plume\nplot(lide$rs)\nplot(point2[point2$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")\n#Plot line density\nplot_ld(lide$inter)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can easily avoid the NO2 perturbation appearing near 270000km to fitting only the power plant plume."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#fitting EMG with stats::optim function\nlide2<-lide$inter\nlide2<-lide2[lide2$id < 270000,]\n\npar<-optim(par=c(3000,20000,235000,.8,1.5),\n           fn=fit_gaussian,\n           dat=lide2$id,\n           res=lide2$density,\n           control = list(maxit=10000)\n)\n\n#predict\nYfit<-pred_gaussian(par$par,lide$inter$id)\n#plot line densities and fitted EMG\nplot_ld(lide$inter,Yfit)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(paste(\"2019 March NO2 emission for Aguirre plant is\",round(compute_no2(parameters=par$par,fun2_output=test_03),4),\"ton/h\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since we have electricity production for that power plant from September to December and according the averaged NO2 map by month in EDA part, we must focus on November and December where the plume seems to be the most neat.  \n  \nRotation :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"vars<-substr(colnames(no2),1,6) %in% c(\"201811\",\"201812\")\n  test_1112<-fun2(pp=\"WRI1026809\",\n             no2 = no2[,vars],\n             cloud_lim = 0.5,\n             cloud_cover = 0.05,\n             radius_stat = 20000,\n             radius_rot=65000,\n             cloud=cloud[,vars],\n             wind_dir=wind_dir[,vars],\n             windspeed = dfw[,vars],\n             point2=point,\n             raster_model = rs)\nrstest<-ag_rasters(test_1112)\nplot(rstest,main=paste(length(test_1112[[1]]),\"rotated days\"))\nplot(point[point$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's try with  different angle restriction"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"par(mfrow=c(2,3))\nfor (i in seq(-90,-40,10)) {\n  rstest<-ag_rasters(test_1112,minangle = i)\n  plot(rstest,main=paste(i))\n}\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Wind direction > -80° seems to be the better option."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"rstest<-ag_rasters(test_1112,minangle= -80)\nlide<-line_dens(raster_result =  rstest,\n                largeur = 10000, #2 * 10km\n                fill=NULL,\n                pts=point,\n                pp=\"WRI1026809\")\n\n#Plot the plume\nplot(lide$rs)\nplot(point2[point2$gppd_idnr==\"WRI1026809\",],add=TRUE,col=\"black\")\nplot_ld(lide$inter)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can avoid the NO2 contamination at about 255000km "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#fitting EMG with stats::optim function\nlide2<-lide$inter\nlide2<-lide2[lide2$id < 250000,]\n\npar<-optim(par=c(3000,20000,235000,.8,0.85),\n           fn=fit_gaussian,\n           dat=lide2$id,\n           res=lide2$density,\n           control = list(maxit=10000)\n)\n\n#predict\nYfit<-pred_gaussian(par$par,lide$inter$id)\n#plot line densities and fitted EMG\nplot_ld(lide$inter,Yfit)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(paste(\"2018 November-December NO2 emission for Aguirre plant is\",round(compute_no2(parameters=par$par,fun2_output=test_1112,minangle=-80),4),\"ton/h\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Aguirre emission factor\n\nNow we have NO2 emission for Aguirre Plant and the real electricity production for 2018 November-December timestamp, we can compute emission factor for that particular plant with that simple formula : EF = NO2 emission / electricity production with NO2 emission in mass unit/hour, electricity production in MWh/hour"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"no2_emission=compute_no2(parameters=par$par,fun2_output=test_1112,minangle=-80)\naguirre_elec=sum(eltot$y2019[eltot$Plant.Name==\"Aguirre Plant\" & eltot$month %in% c(\"November\",\"December\")])\nprint(paste(\"Aguirre electricity generation in November-December 2018 :\",aguirre_elec,\"MWh\"))\nprint(paste(\"Averaged hourly production :\",aguirre_elec/(62*24),\"MWh/h\"))\n\naguirre_ef=1000*no2_emission/(aguirre_elec/(62*24))\nprint(paste(\"Aguirre NO2 emission factor is\",aguirre_ef,\"kg/MWh\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Extrapolate to Puerto Rico\n\nI have found on [this document](https://www.edf.org/sites/default/files/10071_EDF_BottomBarrel_Ch3.pdf) (page 9) some informations relative to RFO/DFO/Natural gas NOx emissions, we can use these data to extrapolate NO2 emissions to the whole island. These data have been compilated from EPA.gov.  \nRFO fuel seems to be much more \"dirty\" than DFO and natural gas. Here are the factors of NO2 emissions relative to RFO :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ef_rate<-data.frame(fuel=c(\"RFO\",\"DFO\",\"NG\",\"COL\",\"SUN\",\"HYC\",\"OTH\",\"WND\"),\n                    factor=c(1,58/166,42/166,NA,0,0,0,0))\nef_rate","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Aguirre plant is fired by a mix of RFO-DFO. Here is the average mix for November-December :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"RFO_rate=sum(eltot$y2019[eltot$Plant.Name==\"Aguirre Plant\" & eltot$month %in% c(\"November\",\"December\") & eltot$AER.Fuel.Type.Code==\"RFO\"])/sum(eltot$y2019[eltot$Plant.Name==\"Aguirre Plant\" & eltot$month %in% c(\"November\",\"December\")])\nDFO_rate=sum(eltot$y2019[eltot$Plant.Name==\"Aguirre Plant\" & eltot$month %in% c(\"November\",\"December\") & eltot$AER.Fuel.Type.Code==\"DFO\"])/sum(eltot$y2019[eltot$Plant.Name==\"Aguirre Plant\" & eltot$month %in% c(\"November\",\"December\")])\n\nprint(paste(\"RFO :\",RFO_rate))\nprint(paste(\"DFO :\",DFO_rate))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"According NO2 emission factors from EPA and NO2 emission factor previously computed, we can separate emission factors for RFO and DFO :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"RFO_ef=aguirre_ef/(RFO_rate+DFO_rate*ef_rate[ef_rate$fuel==\"DFO\",\"factor\"])\nDFO_ef=RFO_ef*ef_rate[ef_rate$fuel==\"DFO\",\"factor\"]\n\nprint(paste(\"Aguirre RFO emission factor :\",RFO_ef))\nprint(paste(\"Aguirre DFO emission factor :\",DFO_ef))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have just a missing factor for coal-fired power plants. AES Puerto Rico plant is pretty young, according to this [site](https://www.modernpowersystems.com/features/featurepuerto-rico-s-first-coal-plant-will-be-among-the-world-s-cleanest/) the plant is equipped by a Selective Non-catalytic Reduction system (SNCR) which can dramatically reduce Nox emissions.  \nI have used [Coal-fired Characteristics and Controls](https://www.epa.gov/airmarkets/power-plant-data-highlights) from EPA.gov to estimate NOx emissions of AES Puerto Rico coal-fired power plant. I have chosen only plants equipped with SNCR system.\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"coal<-fread(\"../input/electricity-production-in-puerto-rico/coalunitcharacteristics2019.csv\",sep=\";\")\ncoal$`2019 NOX (tons)`<-as.numeric(gsub(\" \",\"\",coal$`2019 NOX (tons)`))\ncoal$`2019 Generation (MWh)`<-as.numeric(gsub(\" \",\"\",coal$`2019 Generation (MWh)`))\ncoal$ef<-coal$`2019 NOX (tons)`/coal$`2019 Generation (MWh)`\ncoal_ef=mean(coal$ef[coal$`NOX Controls`==\"Selective Non-catalytic Reduction\"],na.rm=T)*1000000\nprint(paste(\"Averaged Nox emissions factor for coal-fired power plant equipped with SNCR system :\",coal_ef,\"g/MWh\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In papers listed in Methodology part, NOx/NO2 ratio is commonly admitted to be about 1.33 so :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"coal_ef=coal_ef/1.33\nprint(paste(\"Averaged NO2 emissions factor for coal-fired power plant with SNCR :\",coal_ef,\"g/MWh\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, we have all the factors relative to RFO fuel so we can estimate monthly emission factors.\nThe hypothesis here is that individual emission factor (per fuel) is constant through time and power plants.  \nSo, variation come from fuel proportions variation.  \nComplete factors by fuel :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ef_rate[ef_rate$fuel==\"COL\",\"factor\"]<-coal_ef/(RFO_ef*1000)#in gram\nef_rate","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, here is NO2 emission factor for the whole island (don't forget that electricity generation has been estimated for 2019):"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"elm3<-merge(elm3,ef_rate,by.x=\"AER.Fuel.Type.Code\",by.y=\"fuel\")\nelm3$no2<-elm3$value*elm3$factor*RFO_ef\n\nem_factor<-elm3[,.(value=sum(no2)/sum(value)),.(month)]\nggplot(em_factor,aes(x=month,y=value))+\n  geom_point(size=1)+\n  geom_line()+\n  scale_x_continuous(breaks=1:12,label=labs)+\n  theme_bw()+\n  theme(axis.text.x = element_text(angle = 45, hjust = 1))+\n  ylim(c(0,.5))+\n  labs(y=\"kg/MWh\")+\n  geom_text(aes(label=round(value,3)),hjust=.5, vjust=-1)+\n  ggtitle(\"NO2 emission factor\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can also do it for the whole time period :"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"elm4<-elm3[,.(value=sum(value)),.(AER.Fuel.Type.Code)]\nelm4<-merge(elm4,ef_rate,by.x=\"AER.Fuel.Type.Code\",by.y=\"fuel\")\nelm4$no2<-elm4$value*elm4$factor*RFO_ef\nprint(paste(\"Annual averaged NO2 emission factor is :\",sum(elm4$no2)/sum(elm4$value),\"kg/MWh\"))              ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusion\n\nI think this methodology could be scalable depending several conditions :\n* one or multiple PPs should be quite isolated from other NO2 sources (traffic, cities, industries)\n* few clouds area \n* accurate data on electricity production\n\nHere, I have focused on one specific power plant because the others don't meet first and second condition, but the more we can estimate the better.  \n\n## What I have tried\nObviously, I have tried this method on other power plants, without success because of NO2 contamination. I have been surprised by the coal power plant which is completely unvisible.\n\nI've been highly impressed by this recent [paper](https://advances.sciencemag.org/content/5/11/eaax9800.full) which can separate important sources from diffused pollution using TROPOMI data. \nUnfortunately, this method is not appropriated for Puerto Rico area because of the (almost) constant wind direction. It consist of computing daily NO2 *Fluxes*  with the following formula :  \n**F = Vw**  \nwith F the flux, V the NO2 column density and w the wind speed. We got two fluxes, Fu and Fv for u and v components. Then, the localised NO2 sources are revealed by computing the *divergence* (i.e spatial derivative) of the averaged fluxes. \n\n### Temporal mean of U and V fluxes\n\nFirst, we have to spatially interpolate wind data to match NO2 resolution. I have chosen the inverse distance weighting method. Computing time for interpolation is pretty long so for this notebook, I randomly choose 100 days."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"##Read data\n#NO2 data\nfiles<-list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2\")\nrs<-raster(paste0(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/\",files[1]),band=1)\n\nread_raster_band<-function(file,band){\n  rs<-raster(paste0(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/\",file),band=band)\n  rs<-flip(rs,direction=\"y\")\n  result<-list(getValues(rs))\n  names(result)<-substr(file,9,16)\n  return(result)\n}\n\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files,fun=read_raster_band,band=1)\nstopCluster(cl)\nnam<-sapply(ok,names)\nno2<-lapply(ok,function(x) do.call(rbind,x))\nno2<-t(do.call(rbind,no2))\ncolnames(no2)<-nam\n#NaN to NA\nfor (i in 1:ncol(no2)) no2[is.nan(no2[,i]),i]<-NA\nno2<-as.data.frame(no2)\n                      \n#sample 100 days with less missing values as possible\nset.seed(123)\nna_distr<-quantile(colSums(is.na(no2)))            \nno2<-no2[,colSums(is.na(no2)) <= na_distr[\"50%\"]]            \nno2_files<-sample(colnames(no2),100)    \n\nfiles2<-list.files(\"../input/era-interim-wind-data-in-puerto-rico/era_interim\")\n#create output folder\ndir.create(\"../working/era_resampled_idw\")\n#no2 raster model\nfiles<-list.files(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2\")\nrs_model<-raster(paste0(\"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/\",files[1]),band=1)\n\n#create centroids of wind data for interpolation\nera<-raster(paste0(\"../input/era-interim-wind-data-in-puerto-rico/era_interim/\",files2[1]),varname=\"u\",band=1)\npol<-rasterToPolygons(era,dissolve = T)\npts<-st_centroid(st_as_sf(pol))\n\n#interpolation function\nresample_era_idw<-function(no2_file,file,model,pts){\n\n  file<-file[substr(file,14,19)==substr(no2_file,1,6)]  #pick in the right netcdf file\n  era<-raster(paste0(\"../input/era-interim-wind-data-in-puerto-rico/era_interim/\",file))\n  n<-nbands(era)\n  for (j in c(\"u\",\"v\")){\n    for (k in 1:n){\n      \n      \n      rs<-raster(paste0(\"../input/era-interim-wind-data-in-puerto-rico/era_interim/\",file),band=k,varname=paste(j))\n      day<-as.character(k)\n      if (nchar(day)==1) day<-paste0(\"0\",day)\n      date<-paste0(substr(file,14,19),day)\n        \n      if (date != no2_file)  next\n            \n      rasValue=extract(rs, as(pts,\"Spatial\"))\n      pts2<-pts\n      pts2$value<-rasValue\n      pts2<-as(pts2,\"Spatial\")\n      crs(pts2)<-crs(model)\n      mg <- gstat(formula=value~1, locations=pts2, nmax=8)\n      nn <- interpolate(model, mg)\n      \n      \n      \n      \n      writeRaster(nn,paste0(\"../working/era_resampled_idw/\",paste0(j,\"_\",date,\".tif\")),overwrite=TRUE)\n    }\n    \n    \n  }\n}\n\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, {library(\"raster\");library(\"gstat\");library(\"sf\")})\nok<-parLapply(cl,no2_files,fun=resample_era_idw,file=files2,model=rs_model,pts=pts)\nstopCluster(cl)\n\n            \n#u component\nfiles2<-list.files(\"../working/era_resampled_idw\")\nfiles2<-files2[grepl(\"u_\",files2)]\nread_raster<-function(file){\n  rs<-raster(paste0(\"../working/era_resampled_idw/\",file),band=1)\n  result<-list(getValues(rs))\n  names(result)<-substr(file,3,10)\n  return(result)\n}\n\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files2,fun=read_raster)\nstopCluster(cl)\nnam<-sapply(ok,names)\ndfu<-lapply(ok,function(x) do.call(rbind,x))\ndfu<-t(do.call(rbind,dfu))\ncolnames(dfu)<-nam\ndfu<-dfu[,no2_files]#keep the same column order between all tables            \n            \n#read v component\nfiles2<-list.files(\"../working/era_resampled_idw\")\nfiles2<-files2[grepl(\"v_\",files2)]\ncl<-makeCluster(detectCores())\nclusterEvalQ(cl, library(\"raster\"))\nok<-parLapply(cl,files2,fun=read_raster)\nstopCluster(cl)\nnam<-sapply(ok,names)\ndfv<-lapply(ok,function(x) do.call(rbind,x))\ndfv<-t(do.call(rbind,dfv))\ncolnames(dfv)<-nam\ndfv<-dfv[,no2_files]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"options(repr.plot.width=15,repr.plot.height=8)\nrs<-raster(matrix(rowMeans(dfu,na.rm=T),ncol=ncol(rs_model),nrow=nrow(rs_model),byrow=T),template=rs_model)\nplot(rs,main=\"Wind u component spatially interpolated\")\nrs<-raster(matrix(rowMeans(dfv,na.rm=T),ncol=ncol(rs_model),nrow=nrow(rs_model),byrow=T),template=rs_model)\nplot(rs,main=\"Wind v component spatially interpolated\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Average fluxes in mol.m^-1.s^-1"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#flux u\nFu=no2[,no2_files]*dfu\nrsu<-raster(matrix(rowMeans(Fu,na.rm=T),ncol=ncol(rs_model),nrow=nrow(rs_model),byrow = T),template=rs_model)\np<-tm_shape(rsu)+\n  tm_raster(style=\"cont\",palette=\"-RdYlBu\",title=\"Flux u\")+\n  tm_layout(legend.outside=TRUE)\np\n#flux v\nFv=no2[,no2_files]*dfv\nrsv<-raster(matrix(rowMeans(Fv,na.rm=T),ncol=ncol(rs_model),nrow=nrow(rs_model),byrow = T),template=rs_model)\np<-tm_shape(rsv)+\n  tm_raster(style=\"cont\",palette=\"-RdYlBu\",title=\"Flux v\")+\n  tm_layout(legend.outside=TRUE)\np","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Divergence map\n\nWith NO2 point sources in red. (mol.m^-2.s^-1)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#reproject rasters to a plane projection, in meter\nrsu<-projectRaster(rsu,\n                  crs=\"+proj=lcc +lat_1=18.43333333333333 +lat_2=18.03333333333333 +lat_0=17.83333333333333 +lon_0=-66.43333333333334 +x_0=200000 +y_0=200000 +ellps=GRS80 +units=m +no_defs\")\nrsv<-projectRaster(rsv,\n                  crs=\"+proj=lcc +lat_1=18.43333333333333 +lat_2=18.03333333333333 +lat_0=17.83333333333333 +lon_0=-66.43333333333334 +x_0=200000 +y_0=200000 +ellps=GRS80 +units=m +no_defs\")\n#resolution\nresol<-res(rsu)\n#fourth central-finite difference dv/dy\nmat<-matrix((1/(12*resol[2]))*c(0,0,-1,0,0,\n                                0,0,8,0,0,\n                                0,0,0,0,0,\n                                0,0,-8,0,0,\n                                0,0,1,0,0),byrow = T,ncol=5,nrow=5)\n\ntest2<-focal(rsv,w=mat,fun=sum)\n\n#du/dx\nmat<-matrix((1/(12*resol[1]))*c(0,0,0,0,0,\n                                0,0,0,0,0,\n                                1,-8,0,8,-1,\n                                0,0,0,0,0,\n                                0,0,0,0,0),byrow = T,ncol=5,nrow=5)\ntest<-focal(rsu,w=mat,fun=sum)\n\npoint2<-st_transform(point,\n                     crs=\"+proj=lcc +lat_1=18.43333333333333 +lat_2=18.03333333333333 +lat_0=17.83333333333333 +lon_0=-66.43333333333334 +x_0=200000 +y_0=200000 +ellps=GRS80 +units=m +no_defs\")\n#Plot divergence map\np<-tm_shape(test+test2)+\n  tm_raster(style=\"cont\",palette=\"-RdYlBu\",title=\"Divergence\",midpoint=0)+\n  tm_shape(point2[point2$primary_fuel %in% c(\"Oil\",\"Gas\",\"Coal\"),])+\n  tm_dots(size=.3,title=\"Power plants\")+\n  tm_layout(legend.outside=TRUE)\np\n","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}