{"cells":[{"metadata":{"_kg_hide-input":true,"trusted":true,"collapsed":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)\n#install.packages('geojsonio')\n#library(geojsonio)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Goal: Estimate the emission\n\nthe dataset suggested from the competition are following\n\nThe emission : NO2 emission from Sentinel-5p \n\nPollution : Powerplants in the country \n\nTarget period : 1 year (July 1st 2018 - June 29th 2019, 387 raster files for NO2 emissions)\n\n\n# Nitrogen Dioxide (NO2)\n\nAccording to United states Environmental protection agency (EPA), NO2 primarily gets in the air from the **burning of fuel**. NO2 forms from emissions from cars, trucks and buses, power plants, and off-road equipment.\n\nhttps://www.epa.gov/no2-pollution/basic-information-about-no2#What%20is%20NO2\n\n\nPowerplant may be one of big factor for the emission but the problem is that it is **not time variant**. \nOnce the powerplant is located somewhere there's no added information on the amount of emission, it wouldn't be useful to forecast the amount of emission unless future plan of establishing additional powerplant. \n\n\n# Powerplants and NO2 emission\n\n\nSince powerplants are given from the competition, I had a look how powerplant's location and NO2 emission distributed.\n\nEach tif file from Sentinel-5p has 1 week interval each file, and there are 387 files, overlapping many periods. \n\n\n\n\n\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"table(substr(list.files(path = '../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2'),9,14))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"You can see from the above result that each month have more or less 30 images. \nFirst of all, 388 raster images are even hard to see one by one. I would aggregate the data by month to explorer the data easier.\n\nEach tiff file, there are 12 bands. Other bands must be related with NO2 emission too, but first band looks most relavant for me.\n\nTotal vertical column of NO2 (ratio of the slant column density of NO2 and the total air mass factor)\nhttps://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S5P_OFFL_L3_NO2\n\nThe purpose now is the simplyfy the problem and data. Here 1. Aggregate data by month 2. Use only band 1. \n\n\n\nFirst, let's see whole year aggregate NO2 emission in the region\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"files = 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)\nrs2 = as.data.frame(rs1,xy=TRUE)\nggplot()+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_raster(aes(x=x,y=y,fill=layer,alpha=layer),data=rs2)+geom_sf()+theme_minimal()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"** And here's overlapping with powerplant, the size of circle is proportional to the capacitymw **"},{"metadata":{"trusted":true},"cell_type":"code","source":"\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\nlibrary(sf)\nplants_geo <- st_as_sf(plants, coords = c(\"lon\", \"lat\"),crs = 4326, agr = \"constant\")\n\nggplot() +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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since even though it is very accurate where is NO@ emission happening, it is hard to track and also hard to compare where to where.\nAggregate by admin unit always make this job easier"},{"metadata":{"trusted":true},"cell_type":"code","source":"library(sf)\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)\nggplot()+geom_sf(data=pri)\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pri$NO2 <- extract(rs1, pri,fun=mean)\nggplot()+geom_sf(aes(fill=NO2),data=pri)\nggplot()+geom_sf(aes(fill=NO2),data=pri)+geom_sf(data=plants_geo,aes(size=capacity_mw))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n\nIt will be continued..."},{"metadata":{},"cell_type":"markdown","source":"\nAlso I have a doubt that the plant itself. I don't know how other people felt, but the location of platns and Dense NO2 emission area didn't seem that match for me. \n\n\nFor time variant data,number of vehicles on the road will be important. \nFor this purpose, what important for calculation will be road, population, buildings. \nI suggest three datasets which sounds more logical to me for the amount of emission. \n\n1. Nightlight data \nActivities made during the night\n\n\n2. Openstreet map \n\na) Road \n\nb) gas stations? \n\nc) \n\n\n\n3. population grid \n"}],"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}