{"cells":[{"metadata":{},"cell_type":"markdown","source":"### **DS4G: Environmental Insights Explorer**\n#### Exploring alternatives for emissions factor calculations  \n\nAllison Smith"},{"metadata":{},"cell_type":"markdown","source":"**Objective**: Develop a methodology to calculate an average annual historical emissions factor for the sub-national region.  \n\n**Bonus Objective 1**: Smaller time slices of the average historical emissions factors, such as one per month for the 12-month period\n\n**Bonus Objective 2**: Develop methodologies for calculating marginal emissions factors for the sub-national region."},{"metadata":{},"cell_type":"markdown","source":"### Abstract\n\nThe code in the submission will produce an annual emissions factor, monthly emissions factors, and marginal emissions factors for nitrogen dioxide from power plants in Puerto Rico using data from satellite based sensors.  The methodology combines geochemistry, grid analysis, and predictive modeling.  There are also sections discussing sources of nitrogen dioxide, assumptions of the methods, validation that is needed, and extensions to additional sub-national regions."},{"metadata":{},"cell_type":"markdown","source":"### Table of Contents\n\nSources of $NO_2$ emissions\n\nGrid Cell Methodology\n\nSimplified Emissions Factor\n\nConvert .tif to NetCDF (.nc) format\n\nRegrid Data\n\nCalculate the Annual Emissions Factor\n\nCalculate the Monthly Emissions Factors\n\nCalculate the Marginal Emission Factors\n\nConclusions\n\nValidation \n\nExtension to Other Sub-national Regions\n\nAdditional Thoughts\n"},{"metadata":{},"cell_type":"markdown","source":"### Sources of $NO_2$ emissions\n\nPrimary sources of $NO_{2}$ emissions are from power plants, motor vehicles, and silage<sup>1</sup>.  The $NO_{2}$ data from the sensors on the Sentinel 5P OFFL satellite are a combination of emissions from all of these sources. In mathematical terms:  \n\n$NO_{2_{Sentinel}} = NO_{2_{powerplants}} + NO_{2_{vehicles}} + NO_{2_{silage}}$ \n\nThe objective is to calculate an emissions factor for power plants.  In order to achieve this objective, the fraction of emissions measured by the sensor must be attributed to power plants. Power plant emissions are difficult to attribute because of the variability in the characteristics of power plants, which use different types of fuels, have different pollution scrubbers, and have fluctuating emissions depending on demand. \n\n$NO_{2}$ emissions from motor vehicles are less impacted by the emissions from any individual vehicle. Therefore, the collective emissions across many vehicles can be generalized as long as there are no contributions from power plants or silage. \n\n$NO_{2_{vehicles}} = NO_{2_{Sentinel}}$ \n\nwhere $NO_{2_{powerplants}} = 0$ and $NO_{2_{silage}} = 0$\n\nSilage describes the process of fermenting green foliage crops to be used for animal feed.  $NO_{2}$ emissions from silage occur in agricultural areas and have a seasonal component depending on the harvesting patterns of a geographic area.  For the purposes of this analysis, I will make the assumption that grid cells with power plants do not have large agricultural operations.    \n\nTherefore, $NO_2$ from powerplants would be calculated as: \n\n$NO_{2_{powerplants}} = NO_{2_{Sentinel}} - NO_{2_{vehicles}}$\n\nwhere $NO_{2_{silage}} = 0$\n\nReference:   \n<sup>1</sup>https://toxtown.nlm.nih.gov/chemicals-and-contaminants/nitrogen-oxides"},{"metadata":{},"cell_type":"markdown","source":"### Grid cell methodology\n\nSatellite sensors and data assimilation products provide data on grids.  I will regrid the data to 0.01 degree latitude by 0.01 degree longitude grid cells for Puerto Rico which is approximately 1 km x 1 km.  Grid cells will be divided based on primary sources of $NO_2$.  Grid cells with agriculture will be excluded from analyses in order exclude the potential effects of $NO_2$ from silage which is outside the scope of the objectives.  Grid cells without power plants or agriculture will be used to build a predictive model for $NO_2$ emissions from motor vehicles in a grid cell.  Then the predictive model will be used to predict the $NO_2$ emissions from motor vehicles in grid cells with power plants.  The $NO_2$ emissions from motor vehicles can then be subtracted from the total $NO_2$ emissions measured by the satellite to get the $NO_2$ emissions produced by the power plant.  These $NO_2$ emissions can then be used to calculate the emissions factor.  \n\nThere are many assumptions with this method.  \n\nThere is an assumption that all $NO_2$ emissions in grid cells without power plants and agriculture comes from motor vehicles. \n\nThere is an assumption that there is little advection of $NO_2$ from the grid cells near power plants.  I filter the data for the days with the lowest wind speeds but there is still advection. \n\nThere is an assumption that the residence time for $NO_2$ is approximately 1 day. However, the residence time of $NO_2$ varies over time depending on local environmental conditions.  If the residence time is less than 1 day, the emissions will be under estimated. Whereas, if the residence time is more than 1 day , the emissions will be over estimated. "},{"metadata":{},"cell_type":"markdown","source":"### Simplified Emissions Factor\n\nI will use the simplified emissions factor that was presented in the overview for the Kaggle challenge because I do not have data for emissions reduction technology for the power plants in Puerto Rico. \n\nEmissions = Activity_rate * Emissions_factor\n\nRearranging: \n\nEmissions_factor = Emissions/Activity_rate\n\n$NO_2$ Emissions Factor:\n\nEmissions_factor = $NO_{2_{powerplants}}$/estimated_generation_gwh\n\nwhere the $NO_{2_{powerplants}}$ is calculated using the grid cell methodology and the estimated_generation_gwh is from the Global Database of Power Plants. "},{"metadata":{},"cell_type":"markdown","source":"### Convert .tif Files to NetCDF (.nc)\nShell script for converting the s5p_no2 data in the starter pack from multiple .tif files into a single .nc file. I converted all of the data sets because I find the .nc format to be more intuitive and easier to use for analyses.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"##Shell Script\n##requires gdal, nco, cdo to run \n##change directory\n# cd ../data/starter_pack/s5p_no2/\n\n# filenamelist=`ls tif/*.tif`\n# for tif_file in $filenamelist\n# do\n\n#     echo $tif_file\n#     #tif_file=\"tif/s5p_no2_20180701T161259_20180707T175356.tif\"\n\n#     #set names\n#     filename=${tif_file%.tif}\n#     nc_file=nc/${filename#tif/}.nc\n\n#     #convert tif to netCDF\n#     gdal_translate -of NETCDF -co \"FORMAT=NC4\" ${tif_file} convert.nc\n\n#     #rename th variables in the netcdf files\n#     ncrename -v Band1,NO2_column_number_density convert.nc\n#     ncrename -v Band2,tropospheric_NO2_column_number_density convert.nc\n#     ncrename -v Band3,stratospheric_NO2_column_number_density convert.nc\n#     ncrename -v Band4,NO2_slant_column_number_density convert.nc\n#     ncrename -v Band5,tropopause_pressure convert.nc\n#     ncrename -v Band6,absorbing_aerosol_index convert.nc\n#     ncrename -v Band7,cloud_fraction convert.nc\n#     ncrename -v Band8,sensor_altitude convert.nc\n#     ncrename -v Band9,sensor_azimuth_angle convert.nc\n#     ncrename -v Band10,sensor_zenith_angle convert.nc\n#     ncrename -v Band11,solar_azimuth_angle convert.nc\n#     ncrename -v Band12,solar_zenith_angle convert.nc\n\n#     #extract time from the file name and add as a dimension to the netCDF file\n#     timestrT=\"$(cut -d'_' -f3 <<<${filename})\"\n#     timestr=\"${timestrT//T}\"\n#     unix=\"$(date -j -u -f \"%Y%m%d%H%M%S\" \"${timestr}\" +\"%s\")\"\n#     # printf \"/***.nco ncap2 script***/\n#     # defdim(\\\"time\\\",1);\n#     # time[time]=\"$unix\";\n#     # time@long_name=\\\"Time\\\";\n#     # time@units=\\\"seconds since 1970-01-01 00:00:00\\\";\n#     # time@standard_name=\\\"time\\\";\n#     # time@axis=\\\"T\\\";\n#     # time@coordinate_defines=\\\"point\\\";\n#     # time@calendar=\\\"standard\\\";\n#     # /***********/\" > time.nco\n\n#     ncap2 -Oh -s 'defdim(\"time\",1);time[time]=1561830483;time@long_name=\"Time\";time@units=\"seconds since 1970-01-01 00:00:00\";time@standard_name=\"time\";time@axis=\"T\";time@coordinate_defines=\"point\";time@calendar=\"standard\";' convert.nc convert1.nc\n#     date_string=\"$(date -u -r ${unix} +'%Y-%m-%d,%H:%M:%S')\"\n#     cdo settaxis,\"${date_string}\" convert1.nc convert2.nc\n\n#     #save final version of netcdf file \n#     mv convert2.nc ${nc_file}\n\n#     #clean up files\n#     rm convert.nc\n#     rm convert1.nc\n#     #rm time.nco \n\n# done\n\n# cdo mergetime nc/*.nc no2_1year.nc","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Python packages"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport xarray as xr\n#import xesmf as xe\nimport json\nimport matplotlib.pyplot as plt\nimport cartopy.crs as ccrs\nimport seaborn as sns\nimport json\nimport datetime as dt\nfrom sklearn.model_selection import train_test_split\nfrom xgboost import XGBRegressor\nfrom sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Regrid data\n\n**The files that result from running the regridding code in this section were uploaded to Kaggle to be used to calculate the emissions factor.**  They are annual_ds4g_emissions.nc and monthly_ds4g_emissions.nc.  \n\n\nThe code in this section requires xesmf and additional package dependencies to be installed.  The NetCDF files that were generated by the shell script were also too large to upload.   "},{"metadata":{},"cell_type":"markdown","source":"The data sets have different latitude, longitude, and time intervals. The data need to be converted to the same intervals, which I will refer to as grids, for analysis. I regridded the data using the xesmf package in python which makes regridding much easier."},{"metadata":{"trusted":true},"cell_type":"code","source":"#Reusable function for doing the regridding:\n# def make_regridder(ds, ds_base, variable, algorithm='bilinear'):   \n#     if 'latitude' in ds[variable].dims:\n#        ds = ds.rename({'latitude': 'lat', 'longitude': 'lon'}).set_coords(['lon', 'lat'])\n#     ds_regrid = xr.Dataset({'lat': (['lat'], np.arange(np.floor(ds_base['lat'].min().values*10)/10, np.ceil(ds_base['lat'].max().values*10)/10, 0.01)),\n#                      'lon': (['lon'], np.arange(np.floor(ds_base['lon'].min().values*10)/10, np.ceil(ds_base['lon'].max().values*10)/10, 0.01)),\n#                     }\n#                    )\n\n#     regridder = xe.Regridder(ds, ds_regrid, algorithm)\n#     regridder.clean_weight_file()\n#     return regridder","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Create a base grid"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Base grid is based on the s5p_no2 data, i.e. all other data sets will be regrid to the no2 grid\n# ds_s5p = xr.open_dataset('../data/starter_pack/s5p_no2/no2_1year.nc')\n# ds_no2_clouds = ds_s5p[['NO2_column_number_density', 'cloud_fraction']]\n# no2_regridder = make_regridder(ds_no2_clouds, ds_no2_clouds, 'NO2_column_number_density')\n# ds_base_regrid = no2_regridder(ds_no2_clouds)\n# ds_base_regrid = ds_base_regrid.where(ds_base_regrid['NO2_column_number_density']!=0.)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Plot the NO_2 data for a single time point\n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_base_regrid.NO2_column_number_density.isel(time=0).plot(ax=ax, transform=ccrs.PlateCarree());\n# ax.coastlines()\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Add land mask to the base grid\nIn the plot above, there are $NO_2$ measurements both over land and sea.  The sources of $NO_2$ most critical for these analyses are based on land so a high resolution land mask will be useful for doing analyses and creating plots. I typically use a very high resolution SST file (1 km) to create land masks. A freely available data set can be accessed via the link below. Google earth engine only has the GHRSST 4 km data. The data need to be downloaded from NOAA directly to get the 1 km data."},{"metadata":{"trusted":true},"cell_type":"code","source":"#Download Super High Resolution GHRSST SST file (0.01 degree grid)\n#https://coastwatch.pfeg.noaa.gov/erddap/griddap/jplG1SST.nc?SST[(2017-09-13T00:00:00Z):1:(2017-09-13T00:00:00Z)][(17.005):1:(19.005)][(-69.995):1:(-64.005)],mask[(2017-09-13T00:00:00Z):1:(2017-09-13T00:00:00Z)][(17.005):1:(19.005)][(-69.995):1:(-64.005)],analysis_error[(2017-09-13T00:00:00Z):1:(2017-09-13T00:00:00Z)][(17.005):1:(19.005)][(-69.995):1:(-64.005)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Read in SST data, the file is in the Kaggle folder input/ds4g-emissions-nc\n# ds_sea = xr.open_dataset('../input/ds4g-emissions-nc/jplG1SST_e435_8209_9395.nc')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Plot SST to view the grid\n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_sea.SST.isel(time=0).plot(ax=ax, transform=ccrs.PlateCarree());\n# ax.coastlines()\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Regrid SST data to the same grid as the NO2 data\n# sea_regridder = make_regridder(ds_sea, ds_base_regrid, 'SST')\n# ds_sea_regrid = sea_regridder(ds_sea)\n# ds_sea_regrid = ds_sea_regrid.where(ds_sea_regrid['SST']!=0.)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Plot the regridded data\n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_sea_regrid.SST.isel(time=0).plot(ax=ax, transform=ccrs.PlateCarree());\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Create a land mask and add it to the base data set\n# land_ones = ds_sea_regrid.SST.isel(time=0).fillna(1)\n# land_mask = land_ones.where(land_ones ==1.)\n# land_mask = land_mask.where(land_mask.lat<18.5)\n# land_mask = land_mask.drop('time')\n# ds_base_regrid.coords['land_mask'] = land_mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Plot NO2 for just the land\n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_base_regrid['NO2_column_number_density'].isel(time=0).where(ds_base_regrid.land_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Compute daily averages for s5p_no2 data\nThe s5p_no2 data are available on nearly daily intervals. Occasionally, multiple data points are recorded on the same day.  The next step is to standardize the data to a daily time interval. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_base = ds_base_regrid.resample(time='1D').mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Note the changes to the timestamp when the data are plotted\n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_base['NO2_column_number_density'].isel(time=0).where(ds_base.land_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Add weather data to the base grid\nThere were two data sets in the starter pack that included weather data - the gfs data set and the gldas data set.  The gfs data set is based on forecasts and the gldas data set is based on recorded data. There is almost complete overlap in the variables available in each of the data sets. I decided to mostly use the gldas data set because recorded data is better than forecasts for historical assessments.  The gldas data set also has more variables available. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_gldas = xr.open_dataset('../data/starter_pack/gldas/gldas_1year.nc')\n# ds_gldas = ds_gldas.drop('crs')\n# gldas_regridder = make_regridder(ds_gldas, ds_base_regrid, 'Tair_f_inst',  'nearest_s2d')\n# ds_gldas_regrid = gldas_regridder(ds_gldas)\n# ds_gldas_regrid = ds_gldas_regrid.where(ds_gldas_regrid['Tair_f_inst']!=0.)\n# ds_gldas_regrid.coords['land_mask'] = land_mask","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The gldas data are only available for land and the original grid size is fairly large and does not cover many of the coastal sites where power plants are located in Puerto Rico. I filled the coastal areas with the nearest neighbor pixel.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_gldas_regrid_fill = ds_gldas_regrid.ffill(dim='lat')\n# ds_gldas_regrid_fill = ds_gldas_regrid_fill.bfill(dim='lat')\n# ds_gldas_regrid_fill = ds_gldas_regrid_fill.ffill(dim='lon')\n# ds_gldas_regrid_fill = ds_gldas_regrid_fill.bfill(dim='lon')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# #The gldas data without the coastal pixels filled\n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_gldas_regrid['Tair_f_inst'].isel(time=0).where(land_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# #The gldas data with the coastal pixels filled\n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_gldas_regrid_fill['Tair_f_inst'].isel(time=0).where(land_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The gldas data are available at 3 hour intervals. Daily mean, max, and min have the potential to be useful depending on the weather variable being measured. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_gldas_daily_mean = ds_gldas_regrid_fill.resample(time='1D').mean()\n# ds_gldas_daily_max = ds_gldas_regrid_fill.resample(time='1D').max()\n# ds_gldas_daily_min = ds_gldas_regrid_fill.resample(time='1D').min()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_base['gldas_wind_mean'] = ds_gldas_daily_mean['Wind_f_inst']\n# ds_base['gldas_airT_mean'] = ds_gldas_daily_mean['Tair_f_inst']\n# ds_base['gldas_airT_max'] = ds_gldas_daily_max['Tair_f_inst']\n# ds_base['gldas_airT_min'] = ds_gldas_daily_min['Tair_f_inst']\n# ds_base['gldas_lwdown_mean'] = ds_gldas_daily_mean['LWdown_f_tavg']\n# ds_base['gldas_pres_mean'] = ds_gldas_daily_mean['Psurf_f_inst']\n# ds_base['gldas_humidity_mean'] = ds_gldas_daily_mean['Qair_f_inst']\n# ds_base['gldas_heatflux_mean'] = ds_gldas_daily_mean['Qg_tavg']\n# ds_base['gldas_rain_max'] = ds_gldas_daily_max['Rainf_f_tavg']\n# ds_base['gldas_SWdown_max'] = ds_gldas_daily_max['SWdown_f_tavg']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Add GFS wind data to the base grid\nI used GFS for wind speed data because GFS had coastal coverage without filling with the nearest pixel. Wind is particularly important for determining emissions.  Advection of $NO_2$ away from point source grid cells occurs on windy days.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_gfs = xr.open_dataset('../data/starter_pack/gfs/gfs_1year.nc')\n# ds_gfs = ds_gfs.drop('crs')\n# gfs_regridder = make_regridder(ds_gfs, ds_base_regrid, 'temperature_2m_above_ground')\n# ds_gfs_regrid = gfs_regridder(ds_gfs)\n# ds_gfs_regrid = ds_gfs_regrid.where(ds_gfs_regrid['temperature_2m_above_ground']!=0.)\n# ds_gfs_regrid.coords['land_mask'] = land_mask\n# ds_gfs_regrid['wind_speed'] = np.sqrt(np.square(ds_gfs_regrid.u_component_of_wind_10m_above_ground) + np.square(ds_gfs_regrid.v_component_of_wind_10m_above_ground))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_gfs_regrid['wind_speed'].isel(time=6).where(ds_gfs_regrid.land_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_gfs_daily_mean = ds_gfs_regrid.resample(time='1D').mean()\n# ds_base['gfs_wind_mean'] = ds_gfs_daily_mean['wind_speed']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Add night time lights data to the base grid\nAdditional data are needed to calculate the $NO_2$ emissions from motor vehicles.  Ideally, a gridded data set with the number of gas stations per grid cell would be available in Google Earth Engine. This data set is not available (yet). The options that I did find in Google Earth Engine that could be proxies for motor vehicles are population density and nighttime lights. The \npopulation density data product available on Google Earth Engine is the GPWv411: Population Density (Gridded Population of the World Version 4.11). The most recent version of night time lights is the VIIRS Nighttime Day/Night Band Composites Version 1.\n\nI decided to use the night time lights because the population density are modeled estimates based on censuses and the night time lights are data regularly captured with a sensor on a satellite. Data captured regularly via a satellite platform tend to be more reliable and up-to-date.   \n\nhttps://developers.google.com/earth-engine/datasets/catalog/NOAA_VIIRS_DNB_MONTHLY_V1_VCMCFG\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"##Data were downloaded using the Google Earth Engine Interface.  \n# // Create a geometry representing an export region.\n# var geometry = ee.Geometry.Rectangle([-67.4, 17.9, -65.1, 18.6])\n\n# // Load an image\n# var getimagecollection = ee.ImageCollection(\"NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG\")\n#       .filterDate('2018-07-01', '2019-07-01')\n#       .filterBounds(geometry);\n\n# var imageselect = getimagecollection.select('avg_rad').first();\n\n# // var getimage = ee.Image(imageselect).first();\n\n\n# // Export the image, specifying scale and region.\n# Export.image.toDrive({\n# image: imageselect,\n# description: 'ds4g_nighttime_lights2',\n# scale: 1000,\n# region: geometry,\n# });","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_nightlights = xr.open_dataset('../input/ds4g-emissions-nc/VIIRS_nighttime_lights.nc')\n# ds_nightlights2 = ds_nightlights.drop('crs')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# nl_regridder = make_regridder(ds_nightlights2, ds_base_regrid, 'avg_rad')\n\n# ds_nl_regrid = nl_regridder(ds_nightlights2)\n# ds_nl_regrid = ds_nl_regrid.where(ds_nl_regrid['avg_rad']!=0.)\n# ds_nl_regrid.coords['land_mask'] = land_mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_nl_regrid['avg_rad'].where(ds_nl_regrid.land_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_base['night_avg_rad'] = ds_nl_regrid['avg_rad']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Add land cover data to the base grid\nAdditional data are needed to exclude the $NO_2$ emissions from silage.  There are some land cover data sets available on Google Earth Engine. I used the GFSAD1000: Cropland Extent 1km Multi-Study Crop Mask, Global Food-Support Analysis Data becuase it was global and had a high resolution. A global data set is key to applying this method to other sub-national regions.  \n\nhttps://developers.google.com/earth-engine/datasets/catalog/USGS_GFSAD1000_V1"},{"metadata":{"trusted":true},"cell_type":"code","source":"##Data were downloaded using the Google Earth Engine Interface.  \n# // Create a geometry representing an export region.\n# var geometry = ee.Geometry.Rectangle([-67.4, 17.9, -65.1, 18.6])\n\n# // Load an image\n# var getimage = ee.Image(\"USGS/GFSAD1000_V1\")\n\n# var imageselect = getimage.select('landcover');\n\n# // var getimage = ee.Image(imageselect).first();\n\n\n# // Export the image, specifying scale and region.\n# Export.image.toDrive({\n# image: imageselect,\n# description: 'ds4g_landcover',\n# scale: 1000,\n# region: geometry,\n# });","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_landcover = xr.open_dataset('../input/ds4g-emissions-nc/GFSAD1000_landcover.nc')\n# ds_landcover = ds_landcover.drop('crs')\n# land_regridder = make_regridder(ds_landcover, ds_base_regrid, 'landcover_category',  'nearest_s2d')\n# ds_land_regrid = land_regridder(ds_landcover)\n# ds_land_regrid.coords['land_mask'] = land_mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_land_regrid['landcover_category'].where(ds_land_regrid.land_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_base['landcover_category'] = ds_land_regrid['landcover_category']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Add masks for power plants to the base grid"},{"metadata":{},"cell_type":"markdown","source":"The grid cells that contain power plants need to be selected with a mask to analyze the emissions from power plants.  I also created a mask to exclude the grid cells with power plants from the analysis of the emissions from motor vehicles. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#Read in the data \n# plants = pd.read_csv('../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\n# plants = plants[['capacity_mw', 'estimated_generation_gwh', 'primary_fuel', '.geo']]\n\n#Extract the latitude and longitude coordinates from a json string into separate columns\n# coordinates = pd.json_normalize(plants['.geo'].apply(json.loads))['coordinates']\n# plants[['longitude', 'latitude']] = pd.DataFrame(coordinates.values.tolist(), index= coordinates.index)\n# plants.drop('.geo', axis=1, inplace=True)\n\n\n# pd.json_normalize was missing from Kaggle, I uploaded the data with latitude and longitude as separate columns\nplants = pd.read_csv('../input/ds4g-emissions-nc/gppd_120_pr_lat_lon.csv')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Filter for plants that burn fossil fuels and generate NO2\nplants_fossil = plants[plants['primary_fuel'].isin(['Oil', 'Gas', 'Coal'])].copy()\nplants_fossil.reset_index(drop=True, inplace=True)\n# plants_fossil['grid_lon'] = np.nan\n# plants_fossil['position_lon'] = np.ones\n# plants_fossil['grid_lat'] = np.nan\n# plants_fossil['position_lat'] = np.ones\n\n#Map the fossil fuel plants to the nearest grid cells in the base grid \n# lons = ds_base.lon.values\n# a=0 \n# for lon in plants_fossil.longitude:\n#     lon_diff = abs(lon-lons) \n#     plants_fossil.at[a,'grid_lon'] = lons[np.argmin(lon_diff)]\n#     plants_fossil.at[a,'position_lon'] = np.argmin(lon_diff)\n#     a=a+1\n\n# lats = ds_base.lat.values\n# a=0 \n# for lat in plants_fossil.latitude:\n#     lat_diff = abs(lat-lats) \n#     plants_fossil.at[a,'grid_lat'] = lats[np.argmin(lat_diff)]\n#     plants_fossil.at[a,'position_lat'] = np.argmin(lat_diff)\n#     a=a+1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Calculate the number of plants in each grid cell\nplants_fossil['num_plants'] = 1\n#plants_fossil_grid = plants_fossil[['grid_lon', 'grid_lat', 'position_lat', 'position_lon', 'num_plants']].groupby(['grid_lon', 'grid_lat', 'position_lat', 'position_lon'], as_index=False).sum()\nplants_fossil_grid = pd.read_csv('../input/ds4g-emissions-nc/plants_fossil_grid.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Save data for future use\n#plants_fossil.to_csv('plants_fossil.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Convert data frame into a grid that can be added as a mask in the base grid\n# Also creating a mask with the positions of the latitude and longitudes in the grid to use for the marginal emissions analysis\n# plants_mask = 0 * np.ones((ds_base.dims['lat'], ds_base.dims['lon'])) * np.isnan(ds_base.NO2_column_number_density.isel(time=0)) \n# position_lat_id = 0 * np.ones((ds_base.dims['lat'], ds_base.dims['lon'])) * np.isnan(ds_base.NO2_column_number_density.isel(time=0))\n# position_lon_id = 0 * np.ones((ds_base.dims['lat'], ds_base.dims['lon'])) * np.isnan(ds_base.NO2_column_number_density.isel(time=0))\n# plants_mask = plants_mask.drop('time')\n\n# # Create masks for the fossil fuel power plants, \n# # The mask also includes the grid cells immediately surround the grid cell with the power plant\n# for x in plants_fossil_grid.index:\n#     plants_mask[(plants_fossil_grid.at[x,'position_lat']-2):(plants_fossil_grid.at[x,'position_lat']+2),(plants_fossil_grid.at[x,'position_lon']-2):(plants_fossil_grid.at[x,'position_lon']+2)]=1\n#     position_lat_id[(plants_fossil_grid.at[x,'position_lat']-2):(plants_fossil_grid.at[x,'position_lat']+2),(plants_fossil_grid.at[x,'position_lon']-2):(plants_fossil_grid.at[x,'position_lon']+2)]=plants_fossil_grid.at[x,'position_lat']\n#     position_lon_id[(plants_fossil_grid.at[x,'position_lat']-2):(plants_fossil_grid.at[x,'position_lat']+2),(plants_fossil_grid.at[x,'position_lon']-2):(plants_fossil_grid.at[x,'position_lon']+2)]=plants_fossil_grid.at[x,'position_lon']\n\n# # Add the masks to the base grid array    \n# plants_mask = plants_mask.where(plants_mask == 1.)\n# position_lat_id = position_lat_id.where(position_lat_id >= 1.)\n# position_lon_id = position_lon_id.where(position_lon_id >= 1.)\n# ds_base.coords['plants_mask'] = (('lat', 'lon'), plants_mask)\n# ds_base.coords['no_plants_mask'] = ds_base.plants_mask.fillna(0).where((ds_base.plants_mask != 1) & (ds_base.land_mask == 1))\n# ds_base.coords['no_plants_mask']  = ds_base.no_plants_mask + 1\n# ds_base.coords['position_lat_id'] = (('lat', 'lon'), position_lat_id)\n# ds_base.coords['position_lat_id'] = ds_base.position_lat_id.where(ds_base.position_lat_id >= 1)\n# ds_base.coords['position_lon_id'] = (('lat', 'lon'), position_lon_id)\n# ds_base.coords['position_lon_id'] = ds_base.position_lon_id.where(ds_base.position_lon_id >= 1)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot the grid cell areas with power plants \n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_base['NO2_column_number_density'].isel(time=0).where((land_mask==1) & (plants_mask==1)).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot the grid cell ares without power plants\n# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_base['NO2_column_number_density'].isel(time=0).where(ds_base.no_plants_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Calculate daily average for the year"},{"metadata":{},"cell_type":"markdown","source":"$NO_2$ emission concentrations are variable because $NO_2$ can be advected by wind. I filtered the data for the days with the lowest wind speeds and then calculated an annual average for the geographic area.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_base_annual = ds_base.where((ds_base.gfs_wind_mean <= 2)).mean(dim=['time'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_base_annual['NO2_column_number_density'].where((ds_base_annual.land_mask == 1) & (ds_base_annual.no_plants_mask ==1)).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Base grid file is finished. \n#Save a copy to a NetCDF file to be used to calculate the emissions factor.\n#ds_base_annual.to_netcdf('annual_ds4g_emissions.nc')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Calculate daily averages for each month"},{"metadata":{},"cell_type":"markdown","source":"I also calculated monthly averages but the threshold filter for wind speed had to be considerably higher to get adequate geographic coverage for a month. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# ds_base_monthly = ds_base.where((ds_base.gfs_wind_mean <= 5)).resample(time='1M').mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\n# ds_base_monthly['NO2_column_number_density'].isel(time=0).where(ds_base_monthly.land_mask == 1).plot(ax=ax, transform=ccrs.PlateCarree())\n# ax.set_extent([-67.5, -65, 17.5, 19])\n# ax.set_aspect(\"equal\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Save a copy to a NetCDF file to be used to calculate the emissions factor for each month.\n#ds_base_monthly.to_netcdf('monthly_ds4g_emissions.nc')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Calculate the annual emissions factor"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Read in the data set from the saved file created using the code in the regridding section above.\nds = xr.open_dataset('../input/ds4g-emissions-nc/annual_ds4g_emissions.nc')\n\n#Use data set created above.\n#ds = ds_base_annual","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(ds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Create a plot of the NO2 emissions from the data set\nax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\nds.NO2_column_number_density.where(ds.land_mask==1).plot(ax=ax, transform=ccrs.PlateCarree());","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Create a plot of the gfs wind speed from the data set\nax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\nds.wind_speed_mean.where(ds.land_mask==1).plot(ax=ax, transform=ccrs.PlateCarree());","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Create a plot of the NO2 emissions for just the power plant grid cells\nax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\nds.NO2_column_number_density.where(ds.plants_mask==1).plot(ax=ax, transform=ccrs.PlateCarree());","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Build model to predict motor vehicle $NO_2$ emissions for each grid cell\nUse masks and landcover data to filter for grid cells that are located on land, have no power plants nearby, and do not have agricultural operations. "},{"metadata":{"trusted":true},"cell_type":"code","source":"ds['landcover_category'] = ds['landcover_category'].fillna(0)\nds_land = ds.where((ds.land_mask == 1) & \n                   (ds.no_plants_mask == 1) & \n                   (ds.landcover_category == 0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=-65, central_latitude=18))\nds_land.NO2_column_number_density.plot(ax=ax, transform=ccrs.PlateCarree());","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ds_vehicles = ds_land\nds_vehicles = ds_vehicles.drop(['no_plants_mask','plants_mask','land_mask', 'position_lat_id', 'position_lon_id'])\ndf_vehicles = ds_vehicles.to_dataframe().dropna()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Explore features\n# df_vehicles.describe()\n# sns.pairplot(df_vehicles)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Extract features that are important for predicting $NO_2$ emissions in grid cells where the main source of $NO_2$ emissions is motor vehicles. "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_vehicles = df_vehicles[['NO2_column_number_density', 'cloud_fraction', 'night_avg_rad', 'gldas_wind_mean', 'gldas_airT_max', 'gldas_pres_mean', 'gldas_lwdown_mean']].copy()\ndf_vehicles['arcsin_cloud_fraction'] = np.arcsin(np.sqrt(df_vehicles['cloud_fraction']))\ndf_vehicles['log_night_avg_rad'] = np.log(df_vehicles['night_avg_rad'])\ndf_vehicles.drop(['cloud_fraction', 'night_avg_rad'], axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Generate a plot to make sure that none of the variables are highly correlated. \n\ndf_vehicles2 = df_vehicles.copy()\ndf_vehicles2['NO2_column_number_density'] = df_vehicles['NO2_column_number_density']*(10**5) #Kaggle does not always plot very small numbers so rescaling.\nsns.pairplot(df_vehicles2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Build an XGBoost model to predict the NO2_column_number_density for grid cells in which the primary source of $NO_2$ emissions is from motor vehicles. "},{"metadata":{"trusted":true},"cell_type":"code","source":"X = np.array(df_vehicles.drop(['NO2_column_number_density'], axis=1))\ny = np.array(df_vehicles['NO2_column_number_density'])\ny = y*(10**5) #The numbers are very small so rescaling for model training and prediction\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=122)\nmodel = XGBRegressor(learning_rate = 0.01, objective='reg:squarederror', n_estimators = 500)\nmodel.fit(X_train, y_train)\ny_train_pred = model.predict(X_train)\ny_pred = model.predict(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(y_test, y_pred) #Kaggle has limits on small scales, so these should be 10**-5","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Check the performance of the model. The metrics indicate that the data are not overfit during training."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(r2_score(y_train, y_train_pred))\nprint(r2_score(y_test, y_pred))\nprint(mean_squared_error(y_test*10**-5, y_pred*10**-5))\nprint(mean_absolute_error(y_test*10**-5, y_pred*10**-5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Calculate vehicle emissions for the grid cells with power plants"},{"metadata":{},"cell_type":"markdown","source":"I built a model that predicts motor vehicle emissions for a grid cell. Next, I used the model to predict the motor vehicle emissions in power plant grid cells."},{"metadata":{"trusted":true},"cell_type":"code","source":"#Create a data frame with only the power plant grid cells\nds_powerplants = ds.where(ds.plants_mask == 1)\nds_powerplants = ds_powerplants.drop(['no_plants_mask','plants_mask','land_mask'])\ndf_powerplants = ds_powerplants.to_dataframe().dropna()\ndf_powerplants = df_powerplants[['NO2_column_number_density', 'cloud_fraction', 'night_avg_rad', 'gldas_wind_mean', 'gldas_airT_max', 'gldas_pres_mean', 'gldas_lwdown_mean', 'position_lat_id', 'position_lon_id']].copy()\ndf_powerplants['arcsin_cloud_fraction'] = np.arcsin(np.sqrt(df_powerplants['cloud_fraction']))\ndf_powerplants['log_night_avg_rad'] = np.log(df_powerplants['night_avg_rad'])\ndf_powerplants.drop(['cloud_fraction', 'night_avg_rad'], axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_powerplants = np.array(df_powerplants.drop(['NO2_column_number_density', 'position_lat_id', 'position_lon_id'], axis=1))\ny_powerplants = np.array(df_powerplants['NO2_column_number_density'])\ny_powerplants = y_powerplants*(10**5) #The numbers are very small so rescaling for model training and prediction\ndf_powerplants['predict_vehicles'] = model.predict(X_powerplants)*(10**-5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The $NO_2$ from motor vehicles is subtracted from the total $NO_2$ in grid cells with power plants.  The excess $NO_2$ is the $NO_2$ produced by power plants. "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_powerplants['NO2_excess'] = df_powerplants['NO2_column_number_density'] - df_powerplants['predict_vehicles']\n\n# In some cases, the predicted vehicle emissions were higher than the measured emissions in a grid cell. \n# The values in these cells were converted to zero in order to quantify the emissions factor\ndf_powerplants.loc[df_powerplants['NO2_excess']<0, 'NO2_excess'] = 0","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Annual emissions factor for Puerto Rico\nThe satellite sensor passes over Puerto Rico once per day. Therefore, I make the assumption the residence time for $NO_2$ in a grid cell is approximately 1 day and calculate the annual emissions factor. If the residence time is more or less than 1 day, then a more complicated dynamic model would need to be developed.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"#mol/m^2 per day converted to mol/year, a grid cell is approximately 1 km^2\ntotal_emissions = round(df_powerplants['NO2_excess'].sum()*(1000*1000)*365, 2) \nprint(total_emissions)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#gwh for a year using all sources of power including wind, hydro, and solar\ntotal_generated = round(plants['estimated_generation_gwh'].sum())\nprint(total_generated)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#(mol NO2/year)/(mwh/year) = mol NO2/mwh of fossil fuel generated\nemissions_factor = total_emissions/total_generated #(mol NO2/year)/(mwh/year) = mol NO2/mwh of fossil fuel generated\nprint(round(emissions_factor, 3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**The estimated annual emissions factor for Puerto Rico is 0.103 mol $NO_2$ per megawatt hour of power generated.**"},{"metadata":{},"cell_type":"markdown","source":"During data exploration, I calculated the maxium potential generation using the capacity_mw.  In some cases, the estimated_generation_gwh for a power plant in the dataset exceeded the maximum potential annual generation, which is not possible. \n\nThere is an additional reason that I do not think that the estimated_generation_gwh values are correct for individual power plants.  Some of the gas power plants have the exact same estimated_generation_gwh to the 10th decimal place. This is very unlikely. Also, the values could not be that precise.   \n"},{"metadata":{},"cell_type":"markdown","source":"### Calculate the monthly emissions factor\nOne of the objectives is to calculate an emissions factor for each month.  I used the same method that I used to calculate the annual emissions factor.  The wind speed threshold was much higher for the monthly data to get geographic coverage.  Therefore, the emissions facor calculations are likely to be higher. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Read in data from the saved file created using the code in the regridding section above.\nds_monthly = xr.open_dataset('../input/ds4g-emissions-nc/monthly_ds4g_emissions.nc')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Build model to predict motor vehicle $NO_2$ emissions for each grid cell"},{"metadata":{"trusted":true},"cell_type":"code","source":"ds_monthly['landcover_category'] = ds_monthly['landcover_category'].fillna(0)\nds_monthly_vehicles = ds_monthly.where((ds_monthly.land_mask == 1) & \n                   (ds_monthly.no_plants_mask == 1) & \n                   (ds_monthly.landcover_category == 0))\nds_monthly_vehicles = ds_monthly_vehicles.drop(['no_plants_mask','plants_mask','land_mask', 'position_lat_id', 'position_lon_id'])\ndf_monthly_vehicles = ds_monthly_vehicles.to_dataframe().dropna()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Extract important features for predicting $NO_2$ emissions."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_monthly_vehicles = df_monthly_vehicles[['NO2_column_number_density', 'cloud_fraction', 'night_avg_rad', 'gldas_wind_mean', 'gldas_airT_max', 'gldas_pres_mean', 'gldas_lwdown_mean']].copy()\ndf_monthly_vehicles['arcsin_cloud_fraction'] = np.arcsin(np.sqrt(df_monthly_vehicles['cloud_fraction']))\ndf_monthly_vehicles['log_night_avg_rad'] = np.log(df_monthly_vehicles['night_avg_rad'])\ndf_monthly_vehicles.drop(['cloud_fraction', 'night_avg_rad'], axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_monthly_vehicles2 = df_monthly_vehicles.copy()\ndf_monthly_vehicles2['NO2_column_number_density'] = df_monthly_vehicles['NO2_column_number_density']*(10**5)\nsns.pairplot(df_monthly_vehicles2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Build a XGBoost model for predicting monthly $NO_2$."},{"metadata":{"trusted":true},"cell_type":"code","source":"X_monthly_vehicles = np.array(df_monthly_vehicles.drop(['NO2_column_number_density'], axis=1))\ny_monthly_vehicles = np.array(df_monthly_vehicles['NO2_column_number_density'])\ny_monthly_vehicles = y_monthly_vehicles*(10**5) #The numbers are very small so rescaling for model training and prediction\nX_monthly_train, X_monthly_test, y_monthly_train, y_monthly_test = train_test_split(X_monthly_vehicles, y_monthly_vehicles, test_size=0.33, random_state=122)\nmodel_monthly = XGBRegressor(learning_rate = 0.01, objective='reg:squarederror', n_estimators = 2000)\nmodel_monthly.fit(X_monthly_train, y_monthly_train)\ny_monthly_train_pred = model_monthly.predict(X_monthly_train)\ny_monthly_pred = model_monthly.predict(X_monthly_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(y_monthly_test, y_monthly_pred)  #The numbers are very small, multiply numbers in axes by 10**-5 for actual values","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Evaluate the performance of the model.  The model does not overfit the data. "},{"metadata":{"trusted":true},"cell_type":"code","source":"print(r2_score(y_monthly_train, y_monthly_train_pred))\nprint(r2_score(y_monthly_test, y_monthly_pred))\nprint(mean_squared_error(y_monthly_test*10**-5, y_monthly_pred*10**-5))\nprint(mean_absolute_error(y_monthly_test*10**-5, y_monthly_pred*10**-5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Calculate vehicle emissions for the grid cells with power plants\nI built a model that predicts motor vehicle emissions for a grid cell. Next, I used the model to predict the motor vehicle emissions in power plant grid cells."},{"metadata":{"trusted":true},"cell_type":"code","source":"ds_monthly_powerplants = ds_monthly.where(ds_monthly.plants_mask == 1) \nds_monthly_powerplants = ds_monthly_powerplants.drop(['no_plants_mask','plants_mask','land_mask'])\ndf_monthly_powerplants = ds_monthly_powerplants.to_dataframe().dropna()\ndf_monthly_powerplants = df_monthly_powerplants[['NO2_column_number_density', 'cloud_fraction', 'night_avg_rad', 'gldas_wind_mean', 'gldas_airT_max', 'gldas_pres_mean', 'gldas_lwdown_mean']].copy()\ndf_monthly_powerplants['arcsin_cloud_fraction'] = np.arcsin(np.sqrt(df_monthly_powerplants['cloud_fraction']))\ndf_monthly_powerplants['log_night_avg_rad'] = np.log(df_monthly_powerplants['night_avg_rad'])\ndf_monthly_powerplants.drop(['cloud_fraction', 'night_avg_rad'], axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_monthly_powerplants = np.array(df_monthly_powerplants.drop(['NO2_column_number_density'], axis=1))\ny_monthly_powerplants = np.array(df_monthly_powerplants['NO2_column_number_density'])\ny_monthly_powerplants = y_monthly_powerplants*(10**5) #The numbers are very small so rescaling for model training and prediction\ndf_monthly_powerplants['predict_vehicles'] = model_monthly.predict(X_monthly_powerplants)*(10**-5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The $NO_2$ from motor vehicles is subtracted from the total $NO_2$ in grid cells with power plants.  The excess $NO_2$ is the $NO_2$ produced by power plants. "},{"metadata":{"trusted":true},"cell_type":"code","source":"df_monthly_powerplants['NO2_excess'] = df_monthly_powerplants['NO2_column_number_density'] - df_monthly_powerplants['predict_vehicles']\ndf_monthly_powerplants[df_monthly_powerplants['NO2_excess']<0] = 0","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The number of days in each month is extracted from the dates in order ot be used to convert the daily emissions to monthly emissions. "},{"metadata":{"trusted":true},"cell_type":"code","source":"monthly_EF = pd.DataFrame(df_monthly_powerplants['NO2_excess'].groupby('time').sum())\nmonthly_EF['date'] = monthly_EF.index\nmonthly_EF['days_in_month'] = monthly_EF['date'].dt.day\nmonthly_EF.drop('date', axis=1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Total emissions were calculated for each month. I divided the estimated_generation_gwh by 12 months in order to get the generation per month. However, the estimated_generation likely fluctuates on a monthly basis so monthly power generation data would improve calculations of monthly emissions factors.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"#mol/m^2 per day converted to mol/month\nmonthly_EF['total_emissions'] = round(monthly_EF['NO2_excess']*(1000*1000)*monthly_EF['days_in_month'], 2) \nmonthly_EF['total_generated'] = round(plants['estimated_generation_gwh'].sum()/12)\nmonthly_EF['EF'] = monthly_EF['total_emissions']/monthly_EF['total_generated']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Monthly emissions factor for Puerto Rico"},{"metadata":{"trusted":true},"cell_type":"code","source":"monthly_EF","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt = monthly_EF.EF.plot()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**The monthly emissions factors for Puerto Rico ranged from 0.06 to 0.14 from July, 2018 to June, 2019.**"},{"metadata":{},"cell_type":"markdown","source":"### Calculate the marginal emissions factor\nThe emissions factor for each type of fossil fuel can be calculated using grid cell locations grouped by primary fuel type.  Some grid cells had both a gas power plant and an oil power plant so these grid cells are labeled as 'Mixed'."},{"metadata":{"trusted":true},"cell_type":"code","source":"plants_fossil = pd.read_csv('../input/ds4g-emissions-nc/plants_fossil.csv')\nplants_fossil_sum = plants_fossil[['position_lat', 'position_lon', 'estimated_generation_gwh', 'num_plants']].groupby(['position_lat', 'position_lon'], as_index=False).sum()\nplants_fossil_1plant = pd.merge(plants_fossil, plants_fossil_sum[plants_fossil_sum['num_plants']==1])\nplants_fossil_1plant = plants_fossil_1plant[['position_lat', 'position_lon', 'primary_fuel']]\ndf_powerplants_em = df_powerplants[['position_lat_id', 'position_lon_id', 'predict_vehicles', 'NO2_excess']].copy()\ndf_powerplants_loc = df_powerplants_em.groupby(['position_lat_id', 'position_lon_id'], as_index=False).sum()\ndf_powerplants_loc[['position_lat_id','position_lon_id']] = df_powerplants_loc[['position_lat_id','position_lon_id']].applymap(np.int64)\ndf_powerplants_loc.rename({'position_lat_id':'position_lat', 'position_lon_id':'position_lon'}, axis=1, inplace=True)\nplants_1type = pd.merge(plants_fossil_sum, df_powerplants_loc, how='left', on=['position_lat','position_lon'])\nplants_alltypes = pd.merge(plants_1type, plants_fossil_1plant, how='left', on=['position_lat','position_lon'])\nplants_alltypes['primary_fuel'] = plants_alltypes['primary_fuel'].fillna('Mixed')\nplants_primaryfuel = plants_alltypes[['primary_fuel','NO2_excess', 'estimated_generation_gwh', 'num_plants']].groupby(['primary_fuel'], as_index=False).sum()\nplants_primaryfuel['NO2_excess_annual'] = plants_primaryfuel['NO2_excess']*(1000*1000)*365\nplants_primaryfuel['EF'] = plants_primaryfuel['NO2_excess_annual']/plants_primaryfuel['estimated_generation_gwh']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plants_primaryfuel","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**It is possible to calculate marginal emissions factors. There is a lot of variability for the different fuels.** \n\nIn particular, $NO_2$ excess for the coal power plant is very small, and the estimated_generation_gwh is very large.  I calculated the percentage of the estimated_generation_gwh for the coal power plant relative to the total generation for Puerto Rico and found that coal is estimated to generate 92% of the power for Puerto Rico. Therefore, I would expect $NO_2$ emissions for coal to be much higher. I also estimated the fraction of coal compared to the total capacity and found it to be 7% of the total capacity. These calculations do not seem to be in agreement with each other so this may be related to the method that was used to calculate estimated_generation_gwh.  "},{"metadata":{},"cell_type":"markdown","source":"### Conclusions\nEmissions factors were calculated for the sub-national region of Puerto Rico for monthly and yearly time scales.  As more data are collected by the satellite sensor over time, additional analyses can be done to analyze seasonal patterns and trends in $NO_2$ emissions over longer periods of time.  Marginal emissions factors were also calculated but were highly variable and need further investigation.\n\nThe benefit of using a satellite sensor to calculate emissions factors is that data are collected over large geographic areas making it possible to calculate emissions factors over spatial and temporal scales that otherwise would not be possible.  The drawback is that there are many sources of $NO_2$ and attributing the emissions to power plants will have measurement errors that vary with local environmental conditions.   \n\n**All these calculations of emissions factors need validation.** \n"},{"metadata":{},"cell_type":"markdown","source":"### Validation\nThese results have **not** been validated with *in situ* measurements for power plants. If data exist, a sub-national region where $NO_2$ emissions are explicitly measured for power plants could be used to validate the results.    \n"},{"metadata":{},"cell_type":"markdown","source":"### Extension to other sub-national areas\nThis method for calculating emissions factors could be applied to other sub-national areas.  All the data sets used were global data sets.  In particular, it would be interesting calculate emissions factors for areas with known power generation and known emissions factors.  The Global Power Plant Database includes power generation for invdividual power plants when it is available.  \n\nPuerto Rico is an island so there are some unique challenges with measuring emissions in coastal areas where many power plants are located along the land-sea margin.  For example, gldas data are only available for land and did not cover all the areas with power plants in Puerto Rico without extending the data set to surrounding grid cells. Coastal areas also tend to experience daily wind events which advect $NO_2$ emissions away from source grid cells. Therefore, calculating emissions for a region that is not influenced by coastal weather patterns would also be beneficial for evaluating the methodology. "},{"metadata":{},"cell_type":"markdown","source":"### Additional Thoughts\nWind advection transports pollutants away from point sources such as a power plants, which makes it difficult to attribute pollutants measured by satellites. I filtered the measurements for the lowest wind speed days in order to calculate the emissions factor for the power plants.  However, this approach excludes a large proportion of the available data.  A pollution transport simulation model which can trace the movement of pollutants back to original sources would be useful for quantifying the total pollutants emitted by a power plant.  This type of model could also be used to evaluate the residence time of $NO_2$ in a grid cell.  This would be particularly useful for calculating higher temporal resolution, months or weeks, emissions factors. I used a much higher wind speed threshold to calculate the monthly emissions factors because the coverage area was too low using a low wind speed threshold. \n\n\n$NO_{2_{Sentinel}} = NO_{2_{powerplants}} + NO_{2_{vehicles}} + NO_{2_{silage}}$ - $NO_{2_{wind}}$\n\n\nOzone is generated when $NO_2$ reacts with VOCs on hot sunny days.  Depending on the rate of conversion $NO_2$ may be underestimated when these conditions occur.  \n\n$NO_{2_{Sentinel}} = NO_{2_{powerplants}} + NO_{2_{vehicles}} + NO_{2_{silage}}$ - $NO_{2_{wind}}$ - $NO_{2_{ozone}}$\n\nThere are sources of $NO_2$ emissions in addition to power plants, motor vehicles, and silage.  These could also be taken into account. "},{"metadata":{},"cell_type":"markdown","source":"### Acknowledgements\n\nThis analysis utilized data from the following sources:\n\nGlobal Power Plant database by WRI  \nSentinel 5P OFFL NO2 by EU/ESA/Copernicus  \nGlobal Forecast System 384-Hour Predicted Atmosphere Data by NOAA/NCEP/EMC  \nGlobal Land Data Assimilation System by NASA  \nGFSAD1000: Cropland Extent 1km Multi-Study Crop Mask, Global Food-Support Analysis Data by USGS  \nVIIRS Nighttime Day/Night Band Composites Version 1 by NOAA  \nGHRSST Global 1-km Sea Surface Temperature (G1SST) by NOAA  \nGPWv411: Population Density (Gridded Population of the World Version 4.11) by CIESIN  "},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":4}