{"cells":[{"metadata":{},"cell_type":"markdown","source":"

Remote sensing and Emissions Factors

\nThis Notebook shows the work I have done with the S5P and GFS data (I haven't used GLDAS data yet). Takes information also from other Notebooks I made:\n- gppd with additional information (monthly generation and emissions estimate per Power Plant, and the corresponding estimation of global Emissions Factor for 2018) for the second semester of 2018 (first part of the rolling period). See https://www.kaggle.com/ajulian/eia-923-input-nox-emissions-and-ef-reference\n- estimate about the mobile fuel combustion annual emissions. See https://www.kaggle.com/ajulian/activities-and-ghg-precursor-gases\n\n

Some details about \"plot_ee_data_on_map\"

\nFirst I'll describe the function plot_ee_data_on_map, which is extensively used in the notebook.\n\n

Main parameters

\n\"products\" is a list of S5P products; products admitted:\n- \"NO2\" (default): the column will be \"tropospheric_NO2_column_number_density\". EE provides S5P NO2 images since 2018-06-28\n- \"CO\": the column will be \"CO_column_number_density\". EE provides S5P CO images since 2018-06-28\n- \"SO2\": the column will be \"SO2_column_number_density\". EE provides S5P SO2 images since 2018-12-05 (almost 6 months later than NO2 and CO!!)\n- additionally, for \"NO2\" and \"SO2\" the band \"cloud_fraction\" is (can be) displayed, but hidden by default\n\n\"zoom_country\" has two posible values:\n- True (default): Zooms at the Country level; shows images from the product(s) selected and wind arrows around Puerto Rico (unless (lat, long) are provided)\n- False: Zooms at the Power Plant level; shows the product(s) selected around a Power Plant identified by PPindex. The wind arrows direction are given by the u, v wind GFS bands and are shown around the selected location; those u and v bands are used to rotate the products rectangle, thus trying to face the emissions plume.\n\nlat, long: \n- if zoom_country is True and (lat, long) are not provided, by default they are set to the center of Puerto Rico\n- if zoom_country is False, (lat, long) are obtained through PPindex \n\nPPindex is a number [0, 7] identifying one of the eight largest Power Plants\n\n\"fit2country\" is a flag (False by default) to fit the images to the country shape. Makes the images look nice, but the wind pushes many emissions from the North to the sea and with the flag they are not shown, so I almost do not use it.\n\n\"gppd_gdf\" is the Power Plants dataset originally provided with extra information:\n- real power generation and estimated NO2 emissions at the month and Power Plant levels\n (I worked this in another notebook: https://www.kaggle.com/ajulian/eia-923-input-nox-emissions-and-ef-reference)\n- geojson info, taken from https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson\n\n

Other displayable data: Temperature and wind

\n- Temperature: the GFS band \"temperature_2m_above_ground\" is displayed as image, but hidden by default\n- Wind: the GFS \"u\" and \"v\" wind bands are displayed as arrows, and visible by default \n\n

Visualization

\nRegarding color visualization: min (green) and max (red) are computed from the image to be displayed (not so easy in EE as in numpy), and taken from https://www.kaggle.com/ianakoto/emmision-factor\nThis means a \"red\" at the Power Plant level may be an \"orange\" or less when displayed at the Country level. Since this may be sometimes confusing, they are also printed (a colormap bar would be even better)\n\nRemember: there is an icon at the right top of the maps allowing switching on and off each layer."},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport math\nimport folium\nimport matplotlib.pyplot as plt\nimport mplleaflet # matplotlib to leaflet\n\n# Connect to Earth Engine\nimport ee\nfrom kaggle_secrets import UserSecretsClient\nfrom google.oauth2.credentials import Credentials\n\n# Trigger the authentication flow.\n#ee.Authenticate()\n\n# Retrieve your refresh token.\n#!cat ~/.config/earthengine/credentials\n\nuser_secret = \"AJR_EIE_test\" # Your user secret, defined in the add-on menu of the notebook editor\nrefresh_token = UserSecretsClient().get_secret(user_secret)\ncredentials = Credentials(\n None,\n refresh_token=refresh_token,\n token_uri=ee.oauth.TOKEN_URI,\n client_id=ee.oauth.CLIENT_ID,\n client_secret=ee.oauth.CLIENT_SECRET,\n scopes=ee.oauth.SCOPES)\n\n# Initialize GEE\nee.Initialize(credentials=credentials)\ns5p_NOx_clean = None\ns5p_NOx_clean2 = None\ngppd_gdf = None","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"def add_ee_layer(self, ee_image_object, vis_params, name, show=True):\n map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)\n folium.raster_layers.TileLayer(\n tiles = map_id_dict['tile_fetcher'].url_format,\n attr = \"Map Data © Google Earth Engine\",\n name = name,\n overlay = True,\n control = True,\n show = show\n ).add_to(self)\n\ndef rotate_around_point(point, radians, origin=(0, 0)):\n \"\"\"Rotate a point around a given point.\"\"\"\n x, y = point\n ox, oy = origin\n\n qx = ox + math.cos(radians) * (x - ox) + math.sin(radians) * (y - oy)\n qy = oy + -math.sin(radians) * (x - ox) + math.cos(radians) * (y - oy)\n\n return qx, qy\n\n#def rect_rotated(long, lat, side, u_pp, v_pp):\ndef rect_rotated(long, lat, rectangle, u_pp, v_pp):\n # clockwise\n #long1 = long - side/2; lat1 = lat # left bottom\n #long2 = long + side/2; lat2 = lat + side # right top\n long1 = rectangle.bounds().getInfo()['coordinates'][0][0][0]\n lat1 = rectangle.bounds().getInfo()['coordinates'][0][0][1]\n long2 = rectangle.bounds().getInfo()['coordinates'][0][1][0]\n lat2 = rectangle.bounds().getInfo()['coordinates'][0][2][1]\n\n alpha = np.arctan2(u_pp, v_pp)\n \n nlong1, nlat1 = rotate_around_point((long1, lat1), alpha, origin=(long, lat))\n nlong2, nlat2 = rotate_around_point((long1, lat2), alpha, origin=(long, lat))\n nlong3, nlat3 = rotate_around_point((long2, lat2), alpha, origin=(long, lat))\n nlong4, nlat4 = rotate_around_point((long2, lat1), alpha, origin=(long, lat))\n \n rect = ee.Geometry.Polygon([[nlong1, nlat1], [nlong2, nlat2], [nlong3, nlat3], [nlong4, nlat4]])\n return rect\n\ndef add_gfs_wind_arrows_layer(gfs_image_col, band_u, band_v, begin_date, end_date, \n lat, long, zoom_country, rectangle, showMode):\n wind_speed_scale_factor = 50.0\n # dates are transformed to \"previous\" hours before begin_date\n previous = 6 # hours\n begin_date = ee.Date(begin_date).advance(-previous, \"hour\")\n end_date = ee.Date(end_date).advance(-previous, \"hour\")\n #print(begin_date.format().getInfo(), end_date.format().getInfo())\n image_uv = (gfs_image_col\n .filterDate(begin_date, end_date)\n .first())\n \n img = image_uv.addBands(ee.Image.pixelLonLat()) # generates bands \"latitude\" and \"longitude\"\n\n if zoom_country:\n # at the country level, rectangle is given\n scale = 10000\n imgList = img.reduceRegion(reducer=ee.Reducer.toList(),\\\n geometry=rectangle,\\\n maxPixels=1e13,\\\n scale=scale); # WARNING: scale=1000 blocks the Puerto Rico map\n # TODO compute for Puerto Rico\n u_pp_mean = 0\n v_pp_mean = 0\n else:\n # at the PP level, rectangle is computed here:\n # - PP is at the middle of a side\n # - rotation is averaged from few (u, v) \"pixels\"\n # - rotated \"rectangle\" (in fact is a polygon) returns for the emissions image\n scale = 8000\n wind_side = 0.2 # GFS resolution is 0.25 arc??\n lat1 = lat-wind_side/2; long1 = long-wind_side/2\n lat2 = lat+wind_side/2; long2 = long+wind_side/2\n wind_rect = ee.Geometry.Rectangle([long1, lat1, long2, lat2])\n img_min_uvList = image_uv.reduceRegion(reducer=ee.Reducer.toList(),\\\n geometry=wind_rect,\\\n scale=scale);\n u_pp = img_min_uvList.get(band_u).getInfo()\n v_pp = img_min_uvList.get(band_v).getInfo()\n # print(\"El viento es\", u_pp, v_pp)\n imgList = img.reduceRegion(reducer=ee.Reducer.toList(),\\\n geometry=wind_rect,\\\n maxPixels=1e13,\\\n scale=scale);\n \n u_pp_mean = np.mean(u_pp); v_pp_mean = np.mean(v_pp) # useful\n rectangle = rect_rotated(long, lat, rectangle, u_pp_mean, v_pp_mean)\n \n y = imgList.get(\"latitude\").getInfo() # list\n x = imgList.get(\"longitude\").getInfo()\n u_orig = np.array((ee.Array(imgList.get(band_u)).getInfo()))\n v_orig = np.array((ee.Array(imgList.get(band_v)).getInfo()))\n u = u_orig / wind_speed_scale_factor; v = v_orig / wind_speed_scale_factor\n\n x_mesh, y_mesh = np.meshgrid(x, y, sparse=True)\n # print(\"x después\", x_mesh); print(\"y después\", y_mesh)\n U = u.T\n V = v.T\n # print(\"U\", U); print(\"V\", V)\n fig, ax = plt.subplots()\n kw = dict(color='black', alpha=0.8, scale=1)\n q = ax.quiver(x_mesh, y_mesh, U, V, **kw)\n # fig has no data before plotted (ax.quiver) in matplotlib\n gj = mplleaflet.fig_to_geojson(fig=fig)\n\n # feature group allows to have all wind arrows as a layer\n feature_group = folium.map.FeatureGroup(name=\"Wind arrows\")\n for feature in gj['features']:\n if feature['geometry']['type'] == 'Point':\n x_long, y_lat = feature['geometry']['coordinates']\n div = feature['properties']['html']\n\n icon_anchor = (feature['properties']['anchor_x'],\n feature['properties']['anchor_y'])\n\n icon = folium.features.DivIcon(div, icon_anchor=icon_anchor)\n marker = folium.Marker(location=(y_lat, x_long), icon=icon)\n feature_group.add_child(marker)\n # folium.Marker(location=(y_lat, x_long), icon=icon).add_to(Map)\n else:\n msg = \"Unexpected geometry {}\".format\n raise ValueError(msg(feature['geometry']))\n \n return feature_group, rectangle, u_pp_mean, v_pp_mean\n\ndef add_gfs_layers(Map, begin_date, end_date, lat, long, zoom_country, rectangle, showMode):\n dataset = \"NOAA/GFS0P25\"\n band_temp = 'temperature_2m_above_ground'\n band_u = \"u_component_of_wind_10m_above_ground\"\n band_v = \"v_component_of_wind_10m_above_ground\"\n gfs_image_col = (ee.ImageCollection(dataset)\n .select(band_temp, band_u, band_v)\n # in GFS, every 6 hours in a day (00, 06, 12, 18) 384 files are generated; \n # F000 is the first, contains no forecasting but real-time measures \n .filterMetadata(\"system:index\", \"contains\", \"F000\")\n )\n\n feature_group1, rectangle, u_pp_mean, v_pp_mean = \\\n add_gfs_wind_arrows_layer(\n gfs_image_col, band_u, band_v,\n begin_date, end_date, lat, long, zoom_country, rectangle, showMode)\n \n if showMode==True: \n feature_group1.add_to(Map)\n\n vis_temp_params = {\n 'min': -10,\n 'max': 40,\n 'opacity': 0.5,\n 'palette': ['blue', 'purple', 'cyan', 'green', 'yellow', 'red']}\n\n gfs_image = gfs_image_col.filterDate(begin_date, end_date).first().clip(rectangle).select(band_temp)\n Map.add_ee_layer(gfs_image, vis_temp_params, name=\"Temp 2m\", show=False)\n return rectangle, u_pp_mean, v_pp_mean\n\n# Products is a list; the products admitted are:\n# - \"NO2\" (default): the column will be \"tropospheric_NO2_column_number_density\"\n# - \"SO2\": the column will be \"SO2_column_number_density\"\n# - \"CO\": the column will be \"CO_column_number_density\"\n# - additionally, for \"NO2\" and \"SO2\" the band \"cloud_fraction\" is (can be) displayed, but hidden by default\n# \n# zoom_country has two posible values:\n# - True (default): shows the product(s) selected and wind arrows around Puerto Rico\n# - False: shows the product(s) selected around a Power Plant identified by PPidx.\n# the image rectangle follows the wind direction given by the u, v GFS bands\n#\n# lat, long: if zoom_country is True and (lat, long) are not provided, by default they are set to the center of Puerto Rico\n# if zoom_country is False, (lat, long) are obtained through PPindex \n#\n# PPindex is a number [0, 7] identifying one of the eight largest Power Plants\n#\n# gppd_gdf is the Power Plants dataset originally provided with extra information:\n# - real power generation and estimated NO2 emissions at the month and Power Plant levels\n# (I worked this in another notebook)\n# - geojson info, taken from https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson\n#\n# modified from https://www.kaggle.com/paultimothymooney/how-to-get-started-with-the-earth-engine-data\ndef plot_ee_data_on_map(begin_date, end_date, products=[\"NO2\"], \n zoom_country=True, lat=18.232527, long=-66.257565, \n PPindex=0, fit2country=False, opacity=1.0, height=500, side=0.2, showMode=True):\n\n if zoom_country: # zoom at the country level; if (lat, long) are not provided, gets those from Puerto Rico\n zoom_start = 9\n half_x = 1.0667494 # from original kaggle bounding box\n half_y = 0.3323767 # from original kaggle bounding box\n lat1 = lat - 1.5*half_y; long1 = long - 1.5*half_x # more at the West, due to San Juan, and South\n lat2 = lat + 1.5*half_y; long2 = long + half_x # more at the North, due to San Juan \n else: # zoom at a Power Plant level\n print(\"Zoom on\", gppd_gdf.iloc[PPindex][\"name\"])\n lat = gppd_gdf.iloc[PPindex][\"latitude\"]\n long = gppd_gdf.iloc[PPindex][\"longitude\"]\n if showMode==True:\n search_cols = gppd_gdf.columns\n search_cols_gen = [search_cols[i] for i in range(len(search_cols)) \\\n if search_cols[i][:3]==\"Net\"][6:] # second semester\n search_cols_emi = [search_cols[i] for i in range(len(search_cols)) \\\n if search_cols[i][:3]==\"Emi\"][6:-1] # second semester\n\n col_slice_gen = gppd_gdf.iloc[PPindex][search_cols_gen]\n col_slice_emi = gppd_gdf.iloc[PPindex][search_cols_emi]\n #print(col_slice_gen, col_slice_emi)\n print(\"The max gen value for\", gppd_gdf.iloc[PPindex][\"name\"], \"was\", max(col_slice_gen), \n \"MWh and happened in\", search_cols_gen[col_slice_gen.values.argmax()].split()[1], \"2018\")\n print(\"The max emissions value for\", gppd_gdf.iloc[PPindex][\"name\"], \"was\", round(max(col_slice_emi)), \n \"ton and happened in\", search_cols_emi[col_slice_emi.values.argmax()].split()[2], \"2018\")\n zoom_start = 10\n long1 = long - side/2; lat1 = lat # left bottom\n long2 = long + side/2; lat2 = lat + side # right top\n\n\n Map = folium.Map(location=[lat, long], zoom_start=zoom_start, height=height)\n # Map.add_ee_layer = add_ee_layer # does not work\n folium.Map.add_ee_layer = add_ee_layer\n \n rectangle = ee.Geometry.Rectangle([long1, lat1, long2, lat2]) # for PP level will be rotated\n rectangle, u_pp_mean, v_pp_mean = add_gfs_layers(Map, begin_date, end_date, \n lat, long, zoom_country, rectangle, showMode)\n \n for product in products:\n if product in (\"NO2\", \"CO\", \"SO2\"):\n region_scale = 3000 # A nominal scale in meters of the projection to work in.\n #Sentinel-5P Nitrogen Dioxide, Carbon Monoxide or Sulfur Dioxide\n dataset = \"COPERNICUS/S5P/OFFL/L3_\" + product\n # in case product=='NO2', the beginning of the column name is 'tropospheric_', otherwise ''\n column = (product=='NO2')* 'tropospheric_'\n # the end of the column is '_column_number_density' for the three S5p products\n column += product + '_column_number_density'\n if product==\"CO\":\n band_cloud_height = \"cloud_height\"\n else:\n # for NOx and SO2 there are AAI and cloud bands\n band_aai = \"absorbing_aerosol_index\" # not used (yet?)\n band_cloud_fraction = \"cloud_fraction\"\n \n if showMode==True: \n if product==\"NO2\" and s5p_NOx_clean is not None:\n cloud_image = (s5p_NOx_clean\n .select(band_cloud_fraction)\n .filterDate(begin_date, end_date)\n .mean()\n )\n else:\n cloud_image = (ee.ImageCollection(dataset)\n .select(band_cloud_fraction)\n .filterDate(begin_date, end_date)\n .mean()\n )\n\n vis_cloud_params = {\n 'min': 0, # no clouds, black pixel\n 'max': 1, # cloudy, white pixel\n 'opacity': 1,\n 'palette': ['black', 'white']}\n\n Map.add_ee_layer(cloud_image.clip(rectangle), vis_cloud_params, name=\"Clouds_\"+product, show=False)\n\n if product==\"NO2\" and s5p_NOx_clean is not None:\n sat_image = (s5p_NOx_clean\n .select(column)\n .filterDate(begin_date, end_date)\n .mean()\n )\n else:\n sat_image = (ee.ImageCollection(dataset)\n .select(column)\n .filterDate(begin_date, end_date)\n .mean()\n )\n \n if showMode==True:\n # min and max taken from https://www.kaggle.com/ianakoto/emmision-factor\n min_value = sat_image.reduceRegion(ee.Reducer.min(), rectangle, region_scale).getInfo()[column] \n max_value = sat_image.reduceRegion(ee.Reducer.max(), rectangle, region_scale).getInfo()[column]\n if type(min_value) is not float or type(max_value) is not float:\n print(\"Error: there may be no\", product,\"images between\", begin_date, \"and\", end_date, \"or the quality is bad\")\n return\n\n print(product + \" min is\", round(min_value*1e6, 2), \"*10^-6 mol/m^2; max is\", round(max_value*1e6, 2), \"*10^-6 mol/m^2\")\n\n vis_params = {\n 'min': min_value,\n 'max': max_value,\n 'opacity': opacity,\n 'palette': ['green', 'blue', 'yellow', 'red']}\n\n # fit2country clips ee image to the country borders; only makes sense with zoom at the country level\n if zoom_country and fit2country:\n countries = ee.FeatureCollection(\"USDOS/LSIB_SIMPLE/2017\")\n Map.add_ee_layer(sat_image.clip(rectangle).clip(countries), vis_params, product)\n else:\n Map.add_ee_layer(sat_image.clip(rectangle), vis_params, product)\n # Map.add_child(folium.map.LayerControl())\n \n elif product == \"GLDAS\":\n dataset = \"NASA/GLDAS/V021/NOAH/G025/T3H\" # not used (yet?)\n \n if gppd_gdf is not None and showMode==True:\n # taken from https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson\n for plant_lat, plant_long in zip(gppd_gdf.latitude, gppd_gdf.longitude):\n folium.Marker(location = (plant_lat, plant_long), icon = folium.Icon(color='blue')).add_to(Map)\n \n Map.add_child(folium.map.LayerControl()) # needs to be at the end\n \n if showMode==True:\n display(Map)\n else: # computing mode\n return sat_image, rectangle, u_pp_mean, v_pp_mean","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#####

Zoom at the Puerto Rico level

\nFirst we will zoom at the country level, just to see two maps:\n- the first shows an East-West wind pattern typical in Puerto Rico \n- the second is one month later, around the days 14th and 15th, sept. 2017, when the hurricane Maria crossed Puerto Rico and the East-West wind pattern changed completely. \n\nNo S5P products are selected, and in fact can't be because S5P was launched in 2018."},{"metadata":{"trusted":true},"cell_type":"code","source":"# lat = 40.416775; long = -3.703790 # Madrid\n# lat = 25.761681; long=-80.191788 # Miami\n\nbegin_date = '2017-08-14'; end_date = '2017-08-15'\nplot_ee_data_on_map(begin_date, end_date, products=[])\nbegin_date = '2017-09-14'; end_date = '2017-09-15'\nplot_ee_data_on_map(begin_date, end_date, products=[])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

Power Plant level

\nTo work at the Power Plant level we must first load the gppd file with the geo info. \"gppd_gdf\" is the Power Plants dataset originally provided with extra information:\n- real power generation and estimated NO2 emissions at the month and Power Plant levels\n (I worked this in another notebook: https://www.kaggle.com/ajulian/eia-923-input-nox-emissions-and-ef-reference)\n- geojson info, taken from https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson\n\nThe Power Plants have been filtered by power generation, so just the eight largest are included."},{"metadata":{"trusted":true},"cell_type":"code","source":"# geo stuff taken from https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson\nimport pandas as pd\nimport geopandas as gpd\nfrom shapely.geometry import Point\n\ngppd_df = pd.read_csv('/kaggle/input/eia-923-validation-and-nox-emissions-reference/gppd_120_pr_ef.csv') # just the eight bigger plants\n# We can create a new column containing this information\ngppd_df['longitude'] = [float(gppd_df['.geo'][point].split(\"[\")[1][:-2].split(\",\")[0]) for point in range(gppd_df.shape[0])]\ngppd_df['latitude'] = [float(gppd_df['.geo'][point].split(\"[\")[1][:-2].split(\",\")[1]) for point in range(gppd_df.shape[0])]\n\ngeometry_power_plants = [Point(x,y) for x,y in zip(gppd_df.longitude, gppd_df.latitude)]\ngppd_gdf = gpd.GeoDataFrame(gppd_df, crs = {'init': 'epsg: 4326'}, geometry = geometry_power_plants)\n\n# Saving the geodataframe for easy use later\n# gppd_gdf.to_file('Geolocated_gppd_120_pr.geojson', driver='GeoJSON')\nprint(\"These are the eight largest Power Plants considered:\")\ngppd_gdf[\"name\"]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To zoom a Power Plant environment we must set zoom_country=False and specify a PP index (0-7). A text shows the months with more electricity generation and more emissions"},{"metadata":{"trusted":true},"cell_type":"code","source":"zoom_country = False\nPPindex = 1; # 0-7, the 8 largest Power plants, sorted alphabetically\nbegin_date = '2019-04-15'; end_date = '2019-04-16' # max 50\nbegin_date = '2018-07-18'; end_date = '2018-07-20' # max 40\nplot_ee_data_on_map(begin_date, end_date, zoom_country=zoom_country, \n PPindex=PPindex)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The wind information from the GFS u and v bands (represented jointly by those arrows around the Power Plant icon) usually allows positioning the NO2 emissions square quite accurately. It is difficult, however, to separate emissions from San Juan and Palo Seco since they are quite close to each other.\n\n"},{"metadata":{},"cell_type":"markdown","source":"

Emissions Factors

\n

Cleaning images in SP5 NO2 product

\nI have tried to work mostly in EE rather than with the TIFF images provided for several reasons, the main one being the ability to use other products from s5p (SO2, CO) or other satellites data (gfs), and also to work in areas other than Puerto Rico. But first some image filtering needs to be performed to keep only images showing Puerto Rico; since I haven´t been able to do the filtering in EE, I took the image names from the TIFF images provided.\n\nAlso, the number of images is 387, which is bigger than the number of days in a year. The reason is there are pairs of images with the same date; since usually they are partial images which do not cover the whole Puerto Rico, I keep the largest one in each pair."},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport pandas as pd\n\ns5p_files_path = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/'\ns5p_files = os.listdir(s5p_files_path)\nprint(\"Number of TIFF images provided:\", len(s5p_files)) # 387\n\nblack_df = pd.read_csv(\"/kaggle/input/black-list/black_pd.csv\")\nblack_list = list(black_df[\"img_name\"])\n\ndataset = \"COPERNICUS/S5P/OFFL/L3_NO2\"\nband = \"tropospheric_NO2_column_number_density\"\nlat=18.232527; long=-66.257565 # in case lat, long are provided, they are ignored at the country level\nhalf_x = 1.0667494 # from original kaggle bounding box\nhalf_y = 0.3323767 # from original kaggle bounding box\nlat1 = lat - 1.5*half_y; long1 = long - 1.5*half_x # more at the West, due to San Juan, and South\nlat2 = lat + 1.5*half_y; long2 = long + half_x # more at the North, due to San Juan \nrectangle = ee.Geometry.Rectangle([long1, lat1, long2, lat2])\n\ns5p = ee.ImageCollection(dataset)\ns5p_files.sort() # ascending in time\ndate_prev=\"\"\nrepeated = 0\ns5p_images_repeated_list = []\ns5p_images_list = []\nfor i in range(len(s5p_files)-1, -1, -1):\n if s5p_files[i] not in black_list:\n product_id = s5p_files[i].split(\".\")[0].split(\"s5p_no2_\")[1][:-12]\n # print(product_id)\n date = product_id[:8]\n if (date == date_prev):\n s5p_images_repeated_list.append(i) # points to position in list\n repeated += 1\n\n img = s5p.filterMetadata(\"PRODUCT_ID\", \"contains\", product_id).first()\n s5p_images_list.append(img)\n\n date_prev = date\n product_id_prev = product_id\n else:\n s5p_files.pop(i)\nprint(\"Number of ee.Image in the original list:\", len(s5p_images_list))\nprint(\"Number of image pairs with same date:\", repeated)\n\n# we are going to pop one image for each pair in repeated date, \n# so we must reverse the list\ns5p_images_repeated_list.sort(reverse=True)\n# for images in same date, take the largest !! \nprint(\"Removing the smallest image when there are two in same date. SLOW!! (getInfo() inside a loop)\")\n# TODO Avoid loops calling getInfo(), which I think is a communication from server to client\nblack_list = []\nfor i in range(len(s5p_images_repeated_list)):\n s5p_images_list_id = s5p_images_repeated_list[i]\n img = s5p_images_list[s5p_images_list_id]\n img = img.reduceRegion(reducer=ee.Reducer.toList(),geometry=rectangle,\\\n maxPixels=1e13,scale=1000)\n no2 = np.array(ee.Array(img.get(band)).getInfo())\n\n img_prev = s5p_images_list[s5p_images_list_id-1]\n img_prev = img_prev.reduceRegion(reducer=ee.Reducer.toList(),geometry=rectangle,\\\n maxPixels=1e13,scale=1000)\n no2_prev = np.array(ee.Array(img_prev.get(band)).getInfo())\n print(\"Image 1 size:\", no2.shape, \"Image 2 size:\", no2_prev.shape, \n \" (\", i+1, \"/\", len(s5p_images_repeated_list), \")\")\n if (no2_prev.shape > no2.shape):\n s5p_images_list.pop(s5p_images_list_id)\n bad_image_name = s5p_files.pop(s5p_images_list_id)\n else:\n s5p_images_list.pop(s5p_images_list_id-1)\n bad_image_name = s5p_files.pop(s5p_images_list_id-1)\n black_list.append(bad_image_name)\nif len(black_list) > 0:\n black_pd = pd.DataFrame(black_list, columns=[\"img_name\"])\n black_pd.to_csv(\"/kaggle/working/black_pd.csv\", index=False)\n\nprint(black_list)\nprint(\"Number of ee.Image after removing same date images from the original list:\", len(s5p_images_list))\ns5p_NOx_clean = ee.ImageCollection(s5p_images_list)\nprint(\"Created clean collection with\", s5p_NOx_clean.size().getInfo(), \"images\")\nprint(\"Created filename collection with\", len(s5p_files), \"filenames\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have 344 images, which means in the 365 days of the rolling year there are 21 days without an image."},{"metadata":{},"cell_type":"markdown","source":"

NO2 plot for Puerto Rico

\nNow we can plot the mean NO2 emissions (as seen once a day from the s5p satellite) for the rolling year. This is not directly the NOx emissions in the Power Plants or any other emission source (such as vehicles), but it's related."},{"metadata":{"trusted":true},"cell_type":"code","source":"# toBands concatenates in a single image all the bands (in our case just the selected one) \n# from all images in an image collection \nregion_scale = 2000\n#no2_list = s5p_NOx_clean.select(band).toBands().reduceRegion(ee.Reducer.sum(), rectangle, region_scale).toArray().getInfo()\nno2_list = s5p_NOx_clean.select(band).toBands().reduceRegion(ee.Reducer.mean(), rectangle, region_scale).toArray().getInfo()\nprint(\"One point per day:\", len(no2_list))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We can take dates from the imageCollection, which again is slow since uses getInfo() inside a loop.\n# Two ways but both slow\n# datetime_list = s5p_NOx_clean.aggregate_array(\"system:time_start\").getInfo()\n# datetime_list = s5p_NOx_clean.reduceColumns(ee.Reducer.toList(), [\"system:time_start\"]).get('list').getInfo()\n# datetime_list = [ee.Date(millisec).format(\"YYYY-MM-dd HH:mm:ss\").getInfo() for millisec in datetime_list]\ndatetime_list = [f.split(\".\")[0].split(\"s5p_no2_\")[1].split(\"_\")[0].replace(\"T\", \" \") for f in s5p_files]\ndatetime_list = [dt[:4] + \"-\" + dt[4:6] + \"-\" + dt[6:8] + \" \" + dt[9:11] + \":\" + dt[11:13] + \":\" + dt[13:15] for dt in datetime_list]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import plotly.graph_objects as go\nfig = go.Figure()\ndate_list = [datetime.split(\" \")[0] for datetime in datetime_list]\nfig.add_trace(go.Scatter(x=date_list, y=no2_list, name=\"NO2 emissions\"))\nfig.update_layout(title=\"NO2 mean emissions in tropospheric vertical column\",\n yaxis_title=\"NO2 emissions (mol/m^2)\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

Cloud influence

\nClouds make NO2 measures useless; we can take advantage of cloud information in S5p: NO2 and SO2 products have a band named \"cloud_fraction\" where pixels have values between 0 (no cloud) and 1 (cloudy). So let's apply a mask for cloudy days and check the plot again:"},{"metadata":{"trusted":true},"cell_type":"code","source":"cloud_threshold = 0.5 # pixels with value of cloud_fraction > 0.5 (50%) will be masked\ncloud_band = \"cloud_fraction\"\ndef clean_band(image):\n cloud = image.select(cloud_band)\n mask = cloud.lte(cloud_threshold)\n return image.updateMask(mask)\n\ns5p_NOx_clean2 = s5p_NOx_clean.map(clean_band)\nno2_clean_list = s5p_NOx_clean2.select(band).toBands().reduceRegion(ee.Reducer.mean(), rectangle, region_scale).toArray().getInfo()\n#s5p_NOx_clean = s5p_NOx_clean.map(clean_band)\n#no2_clean_list = s5p_NOx_clean.select(band).toBands().reduceRegion(ee.Reducer.mean(), rectangle, region_scale).toArray().getInfo()\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(x=date_list, y=no2_list, name=\"NO2 emissions (as provided)\"))\nfig.add_trace(go.Scatter(x=date_list, y=no2_clean_list, name=\"NO2 emissions (filtering cloud)\"))\nfig.update_layout(title=\"NO2 mean emissions in tropospheric vertical column\",\n yaxis_title=\"NO2 emissions (mol/m^2)\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There has been not much change. Let's check when that maximum happened and see the map."},{"metadata":{"trusted":true},"cell_type":"code","source":"max_date = date_list[no2_clean_list.index(max(no2_clean_list))] # gets the day with the max value\nprint(\"The day with a maximum of average NO2 was\", max_date)\nbegin_date = '2019-04-15'; end_date = '2019-04-17'\nplot_ee_data_on_map(begin_date, end_date)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = go.Figure()\ntime_list = [datetime.split(\" \")[1] for datetime in datetime_list]\ndf = pd.DataFrame({\"no2\": no2_clean_list, \"time\": time_list})\ndf.sort_values(\"time\", inplace=True)\nfig.add_trace(go.Scatter(x=df.time, y=df.no2, mode=\"markers\"))\nfig.update_layout(title=\"NO2 emissions vs. time in tropospheric vertical column\",\n yaxis_title=\"NO2 emissions (mol/m^2)\")\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I cannot see any variation of NO2 measurements with time, which would reflect a diference in lifetime. In fact, this should be easier to check at a Power Plant level but would require daily measures of some activity: electricity generation or emissions."},{"metadata":{},"cell_type":"markdown","source":"

Emissions Factors at the Power Plant and month level

\nIn order to compute the Emissions Factors, we need to estimate the emissions in the Power Plant.\nWe are going to apply a basic box model to convert from NO2 detected by the satellite (in mol/m2) to NO emitted in the Power Plant (in kg/h). This is a classic model which comes from \"Introduction to atmospheric chemistry\" (D.J.Jacob, 1999); other models have been considered but the box model still can be applied, and is very simple.\n\nThe conversion formula is E = NO2 w^2 / (K f tau), where:\n- E is the NO emissions in the Power Plant (kg/h)\n- NO2 is the average NO2 detected in an S5P image (mol/m2)\n- w is the width of the box; in my case, the box rotates each time a measure is taken, so I am taking an average\n- K is a conversion factor, embedded in the rest of conversions\n- f is the rate between NO2 and NO; let's say S5P detects 100 parts of NO2; actually, the Power Plant emitted 132 of NO, out of which 32 did not convert to NO2 and kept as NO. So we have to multiply by 1.32 the NO2 detected to infer the NO emitted. This is a parameter justified by [[Beirle, 2011](http://projects.knmi.nl/publications/fulltexts/1737.full.pdf)]\n- tau is a combination of times: 1/tau = 1/tau_loss + 1/tau_wind; tau_loss is the chemical lifetime, or rate at which the NO2 converts to other gases; it depends on many things: time of day (it is due mainly to sunlight), season of year, etc. There is an interesting figure showing different tau_loss at different times and seasons [[Siyang, 2018](https://www.sciencedirect.com/science/article/pii/S1001074218327426)]; the interesting thing is that tau_loss varies among 4-8 hours, meaning tht most of the emitted NO2 will be lossed during the day even if there is no wind\n- tau_wind is called sometimes the residence time: wind favors the conversion of NO2 into other things, and in Puerto Rico, which is quite windy, is likely to be much lower than tau_loss => the inverse is bigger and we do not have to care too much about tau_loss. The formula for tau_wind is w/2U; since \"w\" usually is in km and U in m/s, a conversion factor must be applied to obtain hours, typical unit for these times"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Parameters:\n# - NO2 in mol/m2\n# - U speed of wind, in m/s (sqrt(u^2+v^2))\n# - w width of box, in km\ndef emissions_jacob(NO2, U, w): # inputs (mol/m2) (m/s) (km); output kg/h\n NO2 = NO2 * 46.0/1000 # 1 mol_NO2/m2 = 46gr/m2 = 46/1000 kg/m2\n w = w * 1000 # km => m\n U = U/(1.0/3600) # m/s => m/h \n tau = w/(2*U) # m/(m/h) => (h)\n f = 1.32\n E = NO2 * w**2 * f /(tau)\n return E","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from haversine import haversine\nimport math\n\ndef emissions_PP(begin_date, end_date, band, PPindex, side, scale):\n w = 21 # km\n\n box, box_rect, u_pp_mean, v_pp_mean = plot_ee_data_on_map(begin_date, end_date, zoom_country=False, \n PPindex=PPindex, side=side, showMode=False)\n boxList = box.reduceRegion(reducer=ee.Reducer.toList(),\\\n geometry=box_rect,\\\n scale=scale);\n NO2_list = boxList.get(band).getInfo()\n NO2_mean = np.mean(NO2_list)\n U = math.sqrt(u_pp_mean**2 + v_pp_mean**2)\n E_kg_h = emissions_jacob(NO2_mean, U, w)\n E_tonne_day = E_kg_h*24/1000\n # print(E_tonne_day, \"tonne/day\")\n return E_tonne_day","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def emissions_slices(PPindex, side, scale):\n emissions_l = [0, 0, 0, 0, 0, 0]\n # last minute problem is not allowing to begin before september 2018\n begin_date = ee.Date(\"2018-09-01\")\n band = \"tropospheric_NO2_column_number_density\"\n\n for m in range(4):\n print(\"Month change\")\n if m==0:\n begin_month=begin_date\n else:\n begin_month = end_month\n end_month = begin_month.advance(1, \"month\").advance(-1, \"day\")\n # print(begin_month.format(\"YYYY-MM-dd\").getInfo(), end_month.format(\"YYYY-MM-dd\").getInfo())\n for slice in range(n_slices):\n if slice==0:\n begin_slice = begin_month\n else:\n begin_slice = end_slice # last end_slice\n if slice==n_slices-1:\n end_slice = end_month\n # TODO end month emissions\n else:\n end_slice = begin_slice.advance(period, \"day\")\n\n emissions_l[m] += period*emissions_PP(begin_slice, end_slice, band, PPindex, side, scale)\n # print(begin_slice.format(\"YYYY-MM-dd\").getInfo(), end_slice.format(\"YYYY-MM-dd\").getInfo())\n return emissions_l","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"PPindex = 1; # 0-7, the 8 largest Power plants, sorted alphabetically\nscale = 3000\nside = 0.2 # This is the default in plot_ee function\nn_months = 6 # just can compare with data from 2018\nn_slices = 6 \nperiod = 5\n\nemissions_l = emissions_slices(PPindex, side, scale)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

Monthly Emissions Factors for Aguirre (sept-dec 2018)

"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"The emissions inferred from the system are\", emissions_l[:4])\nprint(\"The emissions estimated from EIA were\", gppd_gdf.iloc[PPindex][\"Emissions ton September\"],\n gppd_gdf.iloc[PPindex][\"Emissions ton October\"], gppd_gdf.iloc[PPindex][\"Emissions ton November\"],\n gppd_gdf.iloc[PPindex][\"Emissions ton December\"])\nprint(\"The Emissions Factors are\", round(1000*emissions_l[0]/gppd_gdf.iloc[PPindex][\"Netgen September\"], 3),\n round(1000*emissions_l[1]/gppd_gdf.iloc[PPindex][\"Netgen October\"], 3),\n round(1000*emissions_l[2]/gppd_gdf.iloc[PPindex][\"Netgen November\"], 3),\n round(1000*emissions_l[3]/gppd_gdf.iloc[PPindex][\"Netgen December\"], 3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"

The Emissions Factors for Aguirre are 1.008, 0.722, 0.885, 0.72 kg/MWh

\nFor the rest of Power Plants we would proceed the same (but I am running out of time...)"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}