{"cells":[{"metadata":{"_kg_hide-input":true,"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"library(tidyverse) # metapackage with lots of helpful functions\nlibrary(rgdal)\nlibrary(rnaturalearth)\nlibrary(rnaturalearthdata)\nlibrary(ggplot2)\nlibrary(raster)\nlibrary(viridis)\nlibrary(stars)\nlibrary(data.table)\nlibrary(sf)\nlibrary(RColorBrewer)\n#devtools::install_github(\"ropensci/geojsonio\")\n#library(\"geojsonio\")\n\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 <- calc(rs, sum,na.rm = TRUE)\nrs1 <- sum(rs,na.rm = TRUE)\nrs1 <- flip(rs1, direction = \"y\")\nrs2 = as.data.frame(rs1,xy=TRUE)\n#ggplot()+geom_raster(aes(x=x,y=y,fill=layer,alpha=layer),data=rs2)\n\nworld <- ne_countries(scale = \"medium\", returnclass = \"sf\")\nworld=world[which(world$adm0_a3=='PRI'),]\nggplot(data = world) + geom_sf()+geom_raster(aes(x=x,y=y,fill=layer,alpha=layer),data=rs2)\n\n\n\nlist.files(path = '../input/ds4g-environmental-insights-explorer/eie_data/')\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#print(list.files('..input/puerto'))\npri <- st_read('../input/puerto/gadm36_PRI_shp/gadm36_PRI_1.shp')\n#shape <- read_sf(dsn = '../input/puerto/gadm36_PRI_shp/gadm36_PRI_1.shp')#, layer = \"SHAPEFILE\")\n#readRDS('../input/puerto/gadm36_PRI_1_sp.rds')\n#shape <- as.data.frame(shape,xy=TRUE)\n#str(shape)\n#ggplot()+geom_sf(data=pri)\n\npri$NO2 <- extract(rs1, pri,fun=mean)\n#ggplot()+geom_sf(aes(fill=NO2),data=pri)+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")\n#ggplot()+geom_sf(aes(fill=NO2),data=pri)+geom_sf(data=plants_geo,aes(size=capacity_mw))+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")\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\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 <- calc(rs, sum,na.rm = TRUE)    \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#ggplot()+geom_raster(aes(x=x,y=y,fill=layer,alpha=layer),data=rs2)\n\n\n\na <- c()\nb <- c()\nfor(i in 1:12)\n{    \na <- cbind(a,rank(-as.data.frame(pri[,paste0('Y',names(select)[i])])[,1]))\nb <- cbind(b,sum(as.data.frame(pri[,paste0('SY',names(select)[i])])[,1]))\n    }\n#a\n#b\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n# Key figures\n\nToo much information is often poisonous, here are very simplified key figures to start.\nAlso my strategy is very simple too. There will be more detail come afterward, but for me NO2 is following.\n\n## NO2 = Powerplants + Cars\n\n## NO2 Distribution \n\n## Whole year average, NO2 emission longitude corrected (Thank you Rasyid!)\n\n![](https://imgur.com/9YJD4HG.png)\n\n## Which month in Puerto Rico had the most NO2 emission?\n\n#### You can see a clear pattern that between March-July it is significantly higher than August-February.\n#### According to following website, it seems like it is due to the precipitation. Puerto Rico has constant temparature through the year, but the rainy season is more less overlapped with low polluted season. It can be a key clue for finding the factors\nhttps://www.climatestotravel.com/climate/puerto-rico\n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"b <- data.frame(formatC(b, digits = 2, format = \"f\"))\nb <- transpose(b)\nb <- cbind(c(rep(2018,6),rep(2019,6)),c(month.abb[c(7:12,1:6)]),b)\nnames(b) <- c('Year','Month','Total_NO2')\nb$Month <- factor(b$Month, levels = b$Month)\nprint(b)\nggplot(data=b) + geom_bar(aes(x=Month,y=Total_NO2),stat='identity', fill=\"steelblue\")+labs(y='Total NO2 emission each month in Puerto Rico')+ theme_minimal()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n## Which municipality in Puerto Rico had the most NO2 emission? \n\n#### Raster data is very precise, but always very hard to read and summarize. Therefore, I simplify so we can catch pattern from it. I made an average of NO2 emission for each municipality. Puerto Rico has 78 municipalities and it looks like below map. Below table, Total NO2 emission is summing up all the NO2 emission and then divided by the territory size multiply 10^6 to read it easier, and I show rank of each month which municipality emit the most NO2. \n\n#### My conclusion is there's very strong **seasonal effect** in NO2 emission in Puerto Rico. Even the highest emit region is keep changing depends on month. You can see municipalitiy Quebradillas on average emit the most NO2 through the year, but it ranks top only in May and June, and September October it is rather pretty less emitting zone.\n\n#### Toa Baja and Catano seem more constant in a sense of ranking, and the very interesting fact is Toa Baja and Catano are near San Juan (capital) area, and Quebradillas is east coastal area. \n\n![](https://imgur.com/l8xbM3j.png)\n![](https://imgur.com/hPgxXaP.png)\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"a <- data.frame(a)\nresult <- data.frame(cbind(as.character(pri$NAME_1),round(pri$NO2*100000),rank(-pri$NO2),a)[order(pri$NO2,decreasing = TRUE),])\nnames(result) <- c('Municipality','Total_NO2emission','Total_Rank',month.abb[c(7:12,1:6)])\nrownames(result) <- c(seq(1,78))\nresult","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n\n\n\n### What's the relationship with the nightlight and NO2 emission?\n#### Here's monthly average of nightlight for each municipalities\n\n![](https://imgur.com/nscgZBf.png)\n\n![](https://imgur.com/GcBDdy3.png)\n\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"pri_road <- st_read('../input/puerto-road/road_pri.shp')\nnames(pri_road)\nggplot()+geom_sf(aes(col='red'),data=pri_road)+scale_fill_gradient(low = \"#56B1F7\", high = \"#132B43\")+ggtitle(paste0('Y',names(select)[i]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Nightlight \n\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\nhttps://ngdc.noaa.gov/eog/viirs/download_dnb_composites.html\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(list.files(path = '../input/nightlight-pri'))\ninstall.packages('reposTools')\nlibrary(reposTools)\nnew <- untgz(list.files(path = '../input/ds4g-environmental-insights-explorer/nightlight-pri/')[1])","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}