{"cells":[{"metadata":{"id":"ZqQ0WNGv608i"},"cell_type":"markdown","source":"# **Explore Earth Engine Data, Hi Hi**","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"This is my first try on kaggle. At first, I wanted to use R language but later I realized that it does not support Google Earth Engine API. So, I had to start learning both Google Earth Engine API and Python. It is really fun to have a try and here is my findings.\n\nElectricity shutdown in Puerto Rico:\n1. In 2017 September (20/9/2017), [Hurricane Maria](https://www.mercycorps.org/blog/quick-facts-hurricane-maria-puerto-rico) tore through the island and electricity was cut off to 100 percent of the island, and access to clean water and food became limited for most. It wasn’t until August 2018 — nearly a year after the storm — that Puerto Rico Electric Power Authority (PREPA) announced that 100 percent of customers have power restored. \n2. According to [EIA](https://www.eia.gov/state/analysis.php?sid=RQ), Hurricane Maria damaged the one coal-fired electricity generating plant, at Guayama and transmission grid, but the plant resumed generating electricity in February 2018.\n3. According to [cnn news](https://edition.cnn.com/2020/01/09/us/puerto-rico-earthquake-power-outages-satellite-images-trnd/index.html) dated 09/01/2020, there was serious earthquark very closed to Puerto Rico causing power outages on the island. \n\nAssumptions:\n1. It is expected that declining of NO2 emission level could be found during abovementioned events in September 2017 and January 2020. If this is the case, the emission level could be used as reference for background emission rate on the island.\n2. Hydro, solar, and wind power generations are emission-free and will be ignored.\n\nModels:\n1. [Box model](http://acmg.seas.harvard.edu/people/faculty/djj/book/bookchap3.html#pgfId=112721) will be adopted for calculating emission. This will require to determine dimensions of a box around a power plant (pp) for sampling purpose. To calculate emission, I need to know lifetime of NO2 and background NO2 emission on the island. NO2 lifetime can be calculated by wind speed and dimensions of bounding box around powerplants. On the other hand, finding background NO2 emission is tricky.\n\nDefintion:\n1. [Power capacity](https://www.eia.gov/tools/faqs/faq.php?id=101&t=3) - the maximum electric output an electricity generator can produce under specific conditions. \n2. Electricity generation - a generator with 1 megawatt (MW) power generation capacity means the power generator can operate at that capacity consistently for one hour. Alternatively, it means the power generator can produce 1 megawatthour (MWh) of electricity in every one hour.\n3. Estimated yearly generation by plant (in GDPP) - this is estimation done by two methods: scaling information on aggregate geneation by plant size and machine learning appraoch.\n4. [Capacity Factor](https://www.eia.gov/tools/faqs/faq.php?id=101&t=3) is a measure (expressed as a percentage) of how often an electricity generator operates during a specific period of time using a ratio of the actual output to the maximum possible output during that time period. \n5. Emission is mass of emission over time.\n6. Emission Factor is the total emission with specified period (e.g. hour) divided by total electricity generation with the same period of time.\n7. [Watt](https://en.wikipedia.org/wiki/Watt) is unit of power, or energy rate, Joules per second (J/s)\n8. kWh, kilo watt hour is power generated or used in one hour.\n9. Mariginal Emission Factor is the rate of change between change of power generation and change of emission.\n10. Calm wind conditions: wind speed <= 2 meters per second below 500m.\n11. Windy conditions: Wind speed > 2 meters per second.\n12. Wind direction: u is the ZONAL VELOCITY, i.e. the component of the horizontal wind TOWARDS EAST. v is the MERIDIONAL VELOCITY, i.e. the component of the horizontal wind TOWARDS NORTH.\n\nDuration of Data To Be Considered:\nFrom July 2018 to July 2019, however, imagery data may not be available in some time slots.\n\nReading:\n\n[de Foy et al. (2015) Estimates of power plant NOx emissions and lifetimes from OMI NO2 satellite retrievals, Atmoshperic Environment 116(2015)1-11.]]((http://dx.doi.org/10.1016/j.atmosenv.2015.05.056))\n[Liu et al. (2020) A methodology to constrain carbon dixoide emissions from coal-fired power plants using satellite obesrevations of co-emitted nitrogen dioxide, Atmospheric Chemistry and Physics, 20, 99-116.](http://doi.org/10.5194/acp-20-99-2020)\n","execution_count":null},{"metadata":{"id":"6DZldYs2608k"},"cell_type":"markdown","source":"![](http://)**Step 1: Set up mapping functions and start to explore powerplant data**","execution_count":null},{"metadata":{"_kg_hide-input":true,"id":"weDvKvFh608l","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport os\nimport rasterio as rio\nimport folium\n        \n# These are standard functions from folium\n# This use folium to plot a map showing location of powerplants classified by fuel type.\ndef plot_points_on_map(dataframe,begin_index,end_index,latitude_column,latitude_value,longitude_column,longitude_value,zoom):\n    df = dataframe[begin_index:end_index]\n    location = [latitude_value,longitude_value]\n    plot = folium.Map(location=location,zoom_start=zoom)\n    for i in range(0,len(df)):\n        pf = str(df.primary_fuel[i])\n        popup = folium.Popup(pf+'\\n<b>'+str(df.name[i])+'</b>')\n        if (pf == 'Gas'):\n            folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup,icon=folium.Icon(color='green')).add_to(plot)\n        elif (pf == 'Coal'):\n            folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup,icon=folium.Icon(color='black')).add_to(plot)\n        elif (pf == 'Hydro'):\n            folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup,icon=folium.Icon(color='blue')).add_to(plot)\n        elif (pf == 'Solar'):\n            folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup,icon=folium.Icon(color='red')).add_to(plot)\n        elif (pf == 'Wind'):\n            folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup,icon=folium.Icon(color='beige')).add_to(plot)\n        else:\n            folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup,icon=folium.Icon(color='orange')).add_to(plot)\n    return(plot)\n\ndef overlay_image_on_puerto_rico(file_name,band_layer):\n    band = rio.open(file_name).read(band_layer)\n    m = folium.Map([lat, lon], zoom_start=8)\n    folium.raster_layers.ImageOverlay(\n        image=band,\n        bounds = [[18.6,-67.3,],[17.9,-65.2]],\n        colormap=lambda x: (1, 0, 0, x),\n    ).add_to(m)\n    return m\n\ndef plot_scaled(file_name):\n    vmin, vmax = np.nanpercentile(file_name, (5,95))  # 5-95% stretch\n    img_plt = plt.imshow(file_name, cmap='gray', vmin=vmin, vmax=vmax)\n    plt.show()\n\ndef split_column_into_new_columns(dataframe,column_to_split,new_column_one,begin_column_one,end_column_one):\n    for i in range(0, len(dataframe)):\n        dataframe.loc[i, new_column_one] = dataframe.loc[i, column_to_split][begin_column_one:end_column_one]\n    return dataframe","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Step 1.1: Import GPPD dataset using geopandas**","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"I use geopandas and shapely.geometry because it can handle json encoded locational information in a natural way. After creating geometry for each powerplant, I can extract latitude and longitude and store them in separated series.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import geopandas as gpd\nfrom shapely.geometry import mapping, shape\nimport json\nfrom shapely.geometry import shape\npp = gpd.read_file('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\npp['geometry']=pp['.geo'].apply(lambda x: shape(json.loads(x))) #convert str to shape.point and store it in a geometry field.\nprint('converted data type:', pp['geometry'][0]) #Print out to see whether above conversion is okay\npp['latitude']=pp['geometry'].apply(lambda x:(x.y)) #Generate latitude using the point feature in geometry field.\npp['longitude']=pp['geometry'].apply(lambda x:(x.x)) #Generate longitude using the point feature in geometry field","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot the powerplants on map\nlat=18.200178; lon=-66.664513 # set center of map\nend=len(pp.index)\nplot_points_on_map(pp,0,end,'latitude',lat,'longitude',lon,9) # befin index: 0, end index: set to larger than no. of rows is okay","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Examine data type for each variable\npp.dtypes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pp.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#capacity_mw means electrical cap\npower_plants_df = pp.copy()\npower_plants_df['capacity_mw'] = power_plants_df['capacity_mw'].astype('float')\npower_plants_df['commissioning_year'] = power_plants_df['commissioning_year'].apply(lambda x: x.replace('.0', '')).astype('int')\npower_plants_df['estimated_generation_gwh'] = power_plants_df['estimated_generation_gwh'].astype('float')\npower_plants_df[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh','commissioning_year']]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"id":"0eDVlNlC608r","outputId":"72f1995e-1183-4ad9-8659-02ed8806a31b","trusted":true},"cell_type":"code","source":"#sorting by capacity_mw by descending order\npower_plants_df = power_plants_df.sort_values(by=['capacity_mw'],ascending=False).reset_index()\npower_plants_df[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh','commissioning_year']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#power_plants_df['capacity_mw']=power_plants_df['capacity_mw'].astype(float)\n#power_plants_df['estimated_generation_gwh']=power_plants_df['estimated_generation_gwh'].astype(float)\npower_plants_df['capacity_factor'] = power_plants_df['estimated_generation_gwh']/(power_plants_df['capacity_mw']*24*365/1000)\npower_plants_df[['name', 'capacity_mw', 'primary_fuel', 'estimated_generation_gwh', 'capacity_factor']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Need to fix the capacity factor for AES, as it is over 100%, which is impossible. I adopted the ratio suggested by https://www.kaggle.com/c/ds4g-environmental-insights-explorer/discussion/130537","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#Adopt from https://www.kaggle.com/c/ds4g-environmental-insights-explorer/discussion/130537\nsource_capacity_factors = {\"Coal\": 0.55, \"Hydro\": 0.40, \"Gas\": 0.80, \"Oil\": 0.64, \"Solar\": 0.25, \"Wind\": 0.30, \"Nuclear\": 0.85}\n\n# \"force_fix\": - if False, the source_capacity_factors dictionary values are applied only to the \n#                \"estimated_generation_gwh\" values whose Capacity Factor is > 1\n#              - if True, all the \"estimated_generation_gwh\" values are fixed with the source_capacity_factors dictionary values\ndef fix_estimated_generation(gpp_df, source_capacity_factors, force_fix=False):\n    gpp_df[\"capacity_factor\"] = np.where(gpp_df[\"capacity_mw\"] > 0, gpp_df[\"estimated_generation_gwh\"] / (gpp_df[\"capacity_mw\"]*24*365/1000), 0)\n    for idx in range(gpp_df.shape[0]):\n        if (gpp_df.loc[idx, 'capacity_factor'] > 1) or force_fix: \n            gpp_df.loc[idx, 'capacity_factor'] = source_capacity_factors[gpp_df.loc[idx, \"primary_fuel\"]]\n            gpp_df.loc[idx, 'estimated_generation_gwh'] = gpp_df.loc[idx, \"capacity_factor\"] * gpp_df.loc[idx, \"capacity_mw\"] * 24*365/1000\n    return gpp_df\n\nfix_estimated_generation(power_plants_df,source_capacity_factors)\npower_plants_df[['name', 'capacity_mw', 'primary_fuel', 'estimated_generation_gwh', 'capacity_factor']]","execution_count":null,"outputs":[]},{"metadata":{"id":"-13q1iOg608z"},"cell_type":"markdown","source":"## **Step 2: Connect to the Google Earth Engine API**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#ee.Authenticate()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"id":"TFXRp_Hy6082","trusted":true},"cell_type":"code","source":"from kaggle_secrets import UserSecretsClient\nfrom google.oauth2.credentials import Credentials\nimport ee\nimport folium\n\n# Define a method for displaying Earth Engine image tiles on a folium map.\ndef add_ee_layer(self, ee_object, vis_params, name):\n    \n    try:    \n        # display ee.Image()\n        if isinstance(ee_object, ee.image.Image):    \n            map_id_dict = ee.Image(ee_object).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n            ).add_to(self)\n        # display ee.ImageCollection()\n        elif isinstance(ee_object, ee.imagecollection.ImageCollection):    \n            ee_object_new = ee_object.mosaic()\n            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n            ).add_to(self)\n        # display ee.Geometry()\n        elif isinstance(ee_object, ee.geometry.Geometry):    \n            folium.GeoJson(\n            data = ee_object.getInfo(),\n            name = name,\n            overlay = True,\n            control = True\n        ).add_to(self)\n        # display ee.FeatureCollection()\n        elif isinstance(ee_object, ee.featurecollection.FeatureCollection):  \n            ee_object_new = ee.Image().paint(ee_object, 0, 2)\n            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n        ).add_to(self)\n    \n    except:\n        print(\"Could not display {}\".format(name))\n\ndef plot_ee_data_on_map(dataset,column,begin_date,end_date,minimum_value,maximum_value,latitude,longitude,zoom):\n    # https://github.com/google/earthengine-api/blob/master/python/examples/ipynb/ee-api-colab-setup.ipynb\n    folium.Map.add_ee_layer = add_ee_layer\n    vis_params = {\n      'min': minimum_value,\n      'max': maximum_value,\n      'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5']}\n    my_map = folium.Map(location=[latitude,longitude], zoom_start=zoom, height=500)\n    s5p = ee.ImageCollection(dataset).filterDate(\n        begin_date, end_date)\n    my_map.add_ee_layer(s5p.first().select(column), vis_params, 'Color')\n    my_map.add_child(folium.LayerControl())\n    display(my_map)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#!cat ~/.config/earthengine/credentials","execution_count":null,"outputs":[]},{"metadata":{"id":"aqCw05AR6086","trusted":true},"cell_type":"code","source":"from kaggle_secrets import UserSecretsClient\nuser_secret = \"earth_engine\" # 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)\nee.Initialize(credentials=credentials)","execution_count":null,"outputs":[]},{"metadata":{"id":"5XUi-MLC6089"},"cell_type":"markdown","source":"# **Step 3: Preparing data for analysis using the Google Earth Engine API**\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Here I would like to create the following functions for further analysis\n1. A function to accept powerplant location and generate bounding box. This bounding box would be used for sampling different satellite measurements such as NO2 VCD, temperature, wind direction and so on. For NO2 VCD, I would use troposperhic NO2 VCD direct as it can reflect the situation within sensible elevation level.\n2. A function to generate reduce images by whatever summarizing functions such as mean and median by whatever fixed duration such as month, week or day. This will help for comparison to identify difference in emission between days with and without power plant operation. Of course, this is under an assumption that the powerplants were shutdown for some reasons and such scenarios were recorded by the satellites.\n3. Some comparison functions, and I would compare data from two different levels. One is Island-wise scale and another is power-plant scale. By observing island-wise change of NO2 VCD over time, I would like to identify background emission level. Consequently, I could apply the background emission to each power-plant, do comparison over time and find out relevant emission. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#Filter points from power_plants_df\n#Args:\n#        gdf (geoPandas.GeoDataFrame) : the input geodataframe\n#        fuel type : 'Hydro', 'Oil', 'Solar', 'Gas', 'Coal', 'Wind'\n#    Returns:\n#        fc (ee.FeatureCollection) : feature collection (server side)  \ndef FilterPointbyFuelType(df,fuelcolumn,fueltype):\n    #'primary_fuel'\n    fy = df[df[fuelcolumn].isin(fueltype)]\n    features=[]\n    for index, row in fy.iterrows():\n        g=ee.Geometry.Point([row['geometry'].x,row['geometry'].y])\n        feature = ee.Feature(g,{'name':ee.String(row['name']), 'primary_fuel':ee.String(row['primary_fuel']),\\\n                                'capacity_mw':ee.Number(row['capacity_mw']), 'estimated_generation_gwh':ee.Number(row['estimated_generation_gwh']),\\\n                                'capacity_factor':ee.Number(row['capacity_factor'])})\n        features.append(feature)\n    return ee.FeatureCollection(features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Select power plants based on primary fuel type 'Coal' only\nftype = ['Coal']\n# Convert geopandas dataframe to ee.featurecollection\nsel = FilterPointbyFuelType(power_plants_df,'primary_fuel',ftype)\n#sel.getInfo()\nprint('No.of power plants selected:', sel.size().getInfo(), ' for ', ftype)\ntype(sel)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#print(sel.size().getInfo())\n#print(type(sel))\n#Create a bounding box for each powerplant for sampling purpose\n#create a buffer zone for given point type geodataframe, then find its bounding box\n#use the bounding box to define feature\n\n#Define buffer distance, i.e. 1000 m\nsize=1000\n#Define a bounding boxes with dimensions, 2x size * 2x size, i.e. 2000 m x 2000 m\nboundingBoxes = sel.map(lambda f: f.buffer(size).bounds())\n#print(boundingBoxes.)\n#apply the bounding box to extract image\n#Calculate mean value by applying chart by region        \n    \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here I would perform data prepation, i.e. obtain images, reduction, and generate images.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import folium\nfrom folium import plugins\n\n# Define a bounding box for the island and it should be a bit larger to cover enough sea area at South.\nPR_geometry = ee.Geometry.Rectangle([-67.32, 17.70, -65.19, 18.56])\n\ntestimg = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2').filter(ee.Filter.calendarRange(2020, 2020, 'year'))\\\n        .filter(ee.Filter.calendarRange(1, 1, 'month')).filterBounds(PR_geometry)\n\ntestimg.select('tropospheric_NO2_column_number_density')\nprint('no.of images in test', testimg.size().getInfo())\n\nreducedtest = testimg.reduce(ee.Reducer.mean())\n\n# Add EE drawing method to folium.\nfolium.Map.add_ee_layer = add_ee_layer\n\n# Set visualization parameters.\nband_viz = {\n  'min': 0,\n  'max': 0.00010,\n  'palette': ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']\n}\n\n# Create a folium map object.\nmy_map = folium.Map(location=[17.95,-66.15], zoom_start=9, height=500)\n\n# Add the exampleImg model to the map object.\nmy_map.add_ee_layer(reducedtest.select('NO2_column_number_density_mean').clip(PR_geometry), band_viz, 'reducedtest')\n#Add the boundingBoxes\nmy_map.add_ee_layer(boundingBoxes, {'palette': ['FF0000']},'boundingBoxes') #FF0000: red, FFFF00: yellow\n\n# Add a layer control panel to the map.\nmy_map.add_child(folium.LayerControl())\nplugins.Fullscreen().add_to(my_map)\n\n# Display the map.\ndisplay(my_map)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Here I will select images by a date range, and bounding box, then reduce it by month. However, not all the months have images,\n# so I need to remove entry without image before further processing.\n\nimport folium\nfrom folium import plugins\n\n# Define a bounding box for the island and it should be a bit larger to cover enough sea area at South.\nPR_geometry = ee.Geometry.Rectangle([-67.32, 17.70, -65.19, 18.56])\n\n# Define time range\nstartyear = 2017\nendyear = 2020\n \n# Set date in ee date format\nstartdate = ee.Date.fromYMD(startyear,7,1) # 2017/7/1\nenddate = ee.Date.fromYMD(endyear,2,28) # 2020/2/28\n\n# Select sentinel 5P NO2 images filtered by PR_geometry\nimgcollection = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2').filterBounds(PR_geometry)\n\nprint('total no. of selected images:', imgcollection.size().getInfo())\nprint(imgcollection.first().bandNames().getInfo())\n# There are 12 bands.\n#['NO2_column_number_density', 'tropospheric_NO2_column_number_density', 'stratospheric_NO2_column_number_density', \n# 'NO2_slant_column_number_density', 'tropopause_pressure', 'absorbing_aerosol_index', 'cloud_fraction', 'sensor_altitude', \n# 'sensor_azimuth_angle', 'sensor_zenith_angle', 'solar_azimuth_angle', 'solar_zenith_angle']\n\n# Add EE drawing method to folium.\nfolium.Map.add_ee_layer = add_ee_layer\n\n# Set visualization parameters.\nband_viz = {\n  'min': 0,\n  'max': 0.00010,\n  'palette': ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']\n}\n# Create a folium map object.\nmy_map = folium.Map(location=[17.95,-66.15], zoom_start=9, height=500)\n\n# Create list of dates for time series\n# first, find out number of months between end date and start date\nn_months = enddate.difference(startdate, 'month').round()\n# Make a day of month sequence from 0 to n_months-1 with a 1-month step.\nmonths = ee.List.sequence(0, n_months.subtract(1), 1)\n\n# Use start date as anchor point and generate list of date, each entry will be by n-month advancement\ndef make_datelist(n):\n    return startdate.advance(n, 'month')\n\n# Here apply the make_datelist function and generate dates\ndates =  months.map(make_datelist)\n\n# getm will reduce images by averaging values for images within the same month acccording to given date.\ndef getm(d1):\n    d1 = ee.Date(d1)\n    m = d1.get('month')\n    y = d1.get('year')\n    # Create image by mean reducer\n    s = imgcollection.filter(ee.Filter.calendarRange(y, y, 'year'))\\\n    .filter(ee.Filter.calendarRange(m, m, 'month')).mean()\n    return s\\\n           .set('month', m)\\\n           .set('year', y)\\\n           .set('count', s.bandNames().length())\n\n# This will create a list of images\nlistofimages = dates.map(getm)\n# Create an ImageCollection\nmonthlyCol = ee.ImageCollection.fromImages(listofimages)    \n\n# Remove null band images by selecting count greater than or equal to 1\nmonthlyCol = monthlyCol.filter(ee.Filter.gte('count', 1))\n\n# print out no. of images\nprint('no. of remaindar images with bands:', monthlyCol.size().getInfo())\n\n# print('bandnames: ', monthlyCol.first().bandNames().getInfo())\n\n# Call the first image for visual check\nimg = ee.Image(monthlyCol.first())\n\n# Add img to folium\nmy_map.add_ee_layer(img.select('tropospheric_NO2_column_number_density').clip(PR_geometry), band_viz, 'tropospheric_NO2_column_number_density')\n# Add bounding Boxes\nmy_map.add_ee_layer(boundingBoxes, {'palette': ['FF0000']},'boundingBoxes') #FF0000: red, FFFF00: yellow\n# Add study area\nmy_map.add_ee_layer(PR_geometry, {'palette': ['FFFF00']},'PR') #FF0000: red, FFFF00: yellow\n\n# Add a layer control panel to the map.\nmy_map.add_child(folium.LayerControl())\nplugins.Fullscreen().add_to(my_map)\n\n# Display the map.\ndisplay(my_map)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And now, I have boundingbox for each selected power plant (the red boxes), and monthly average NO2 measurements for the entire island from 7/2017 to 2/2020, i.e. 32 images. Following, I need to find average tropospheric NO2 measurement for each box or you can think of, for each powerplant. Then I need to compile it as a panda dataframe for time series analysis. So, in the dataframe, it has the following fields: Power plant name, NO2 value (mean), Month, and Year.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# ==========================================================================\n# Function to Convert Feature Classes to Pandas Dataframe\n# Adapted from: https://events.hpc.grnet.gr/event/47/material/1/12.py\ndef fc2df(fc):\n    # Convert a FeatureCollection into a pandas DataFrame\n    # Features is a list of dict with the output\n    features = fc.getInfo()['features']\n\n    dictarr = []\n\n    for f in features:\n        # Store all attributes in a dict\n        attr = f['properties']\n        # and treat geometry separately\n        attr['geometry'] = f['geometry']  # GeoJSON Feature!\n        # attr['geometrytype'] = f['geometry']['type']\n        dictarr.append(attr)\n\n    df = gpd.GeoDataFrame(dictarr)\n    # Convert GeoJSON features to shape\n    df['geometry'] = map(lambda s: shape(s), df.geometry)\n    return df\n\n# ==========================================================================\n# Function to iterate over image collection, returning a pandas dataframe\ndef extract_point_values(img_id, polys):\n    image = ee.Image(img_id)\n\n    fc_image_red = image.reduceRegions(collection=polys,\n                                  reducer=ee.Reducer.mean(),\n                                  scale=30)\n\n    # Convert to Pandas Dataframe\n    df_image_red = fc2df(fc_image_red)\n\n    # Add Date as Variable\n    df_image_red['band'] = image.bandNames()\n\n    return df_image_red\n\n# ==========================================================================\n# Function to iterate over image collection, returning a pandas dataframe\ndef _extract_point_values_map(img):\n    image = ee.Image(img)\n    # fc_image_red is feature collection containing mean value after reduceRegions\n    fc_image_mean = image.reduceRegions(collection=boundingBoxes,\n                                  reducer=ee.Reducer.mean(),\n                                  scale=30)\n    # Convert to Pandas Dataframe\n    df_image_mean = fc2df(fc_image_mean)\n    \n    # Select required field only\n    df_image_mean = df_image_mean[['NO2_column_number_density',\n     'NO2_slant_column_number_density',\n     'absorbing_aerosol_index',\n     'capacity_factor',\n     'capacity_mw',\n     'cloud_fraction',\n     'estimated_generation_gwh',\n     'name',\n     'primary_fuel',\n     'stratospheric_NO2_column_number_density',\n     'tropopause_pressure',\n     'tropospheric_NO2_column_number_density',\n     'geometry']]\n    # Add Date as Variable    \n    df_image_mean['month'] = image.getInfo()['properties']['month']\n    df_image_mean['year'] = image.getInfo()['properties']['year']\n    \n    return df_image_mean\n# ==========================================================================\n\n# ==========================================================================\n# Function to iterate over image collection, returning a pandas dataframe\n# mean\ndef extract_point_values_map2(img):\n    image = ee.Image(img)\n    # fc_image_red is feature collection containing mean value after reduceRegions\n    fc_image_mean = image.reduceRegions(collection=boundingBoxes,\n                                  reducer=ee.Reducer.mean(),\n                                  scale=30)\n    # Convert to Pandas Dataframe\n    df_image_mean = fc2df(fc_image_mean) \n    # Add Date as Variable    \n    df_image_mean['month'] = image.getInfo()['properties']['month']\n    df_image_mean['year'] = image.getInfo()['properties']['year']\n    \n    return df_image_mean\n# ==========================================================================\n\n# ==========================================================================\n# Function to iterate over image collection, returning a pandas dataframe\n# min\ndef extract_point_values_map3(img):\n    image = ee.Image(img)\n    # fc_image_red is feature collection containing mean value after reduceRegions\n    fc_image_min = image.reduceRegions(collection=boundingBoxes,\n                                  reducer=ee.Reducer.min(),\n                                  scale=30)\n    # Convert to Pandas Dataframe\n    df_image_min = fc2df(fc_image_min) \n    # Add Date as Variable    \n    df_image_min['month'] = image.getInfo()['properties']['month']\n    df_image_min['year'] = image.getInfo()['properties']['year']\n    \n    return df_image_min\n# ==========================================================================\n\n\n# ==========================================================================\n# Function to iterate over image collection, returning a pandas dataframe\ndef reduceImg(img):\n    # fc_image_red is a mean value after reduceRegions\n    image_mean = image.reduceRegions(collection=boundingBoxes,\n                                  reducer=ee.Reducer.mean(),\n                                  scale=30)\n    image_mean = image_mean.set()\n    return image_mean\n# ==========================================================================","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"varImage = monthlyCol\nlistOfImages = varImage.toList(varImage.size())\n#### Create Initial Pandas Dataframe\n\ndf_all_mean = extract_point_values_map2(listOfImages.get(0))\ndf_all_min = extract_point_values_map3(listOfImages.get(0))\nnumOfImage = varImage.size().getInfo()\nli = range(1,numOfImage)\nfor i in li:\n    df_all_mean = df_all.append(extract_point_values_map2(listOfImages.get(i)))\n    df_all_min = df_all.append(extract_point_values_map3(listOfImages.get(i)))\n\n#Convert ['year','month'] to date\ndf_all_mean['DATE'] = pd.to_datetime(df_all_mean[['year', 'month']].assign(DAY=1)).dt.to_period('M')\ndf_all_min['DATE'] = pd.to_datetime(df_all_min[['year', 'month']].assign(DAY=1)).dt.to_period('M')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#df_all_mean.drop(df_all_mean.index, inplace=True)\n#df_all_min.drop(df_all_mean.index, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#### Display Results\ndf_all_min","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_all_mean","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_all_min[['NO2_column_number_density',\n 'NO2_slant_column_number_density',\n 'absorbing_aerosol_index',\n 'capacity_factor',\n 'capacity_mw',\n 'cloud_fraction',\n 'estimated_generation_gwh',\n 'name',\n 'primary_fuel',\n 'stratospheric_NO2_column_number_density',\n 'tropopause_pressure',\n 'tropospheric_NO2_column_number_density',\n 'geometry',\n 'month',\n 'year']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_all_mean[['NO2_column_number_density',\n 'NO2_slant_column_number_density',\n 'absorbing_aerosol_index',\n 'capacity_factor',\n 'capacity_mw',\n 'cloud_fraction',\n 'estimated_generation_gwh',\n 'name',\n 'primary_fuel',\n 'stratospheric_NO2_column_number_density',\n 'tropopause_pressure',\n 'tropospheric_NO2_column_number_density',\n 'geometry',\n 'month',\n 'year']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Finding Background NO2**\n\nTo determine NO2 emission generating from power plants, it is crucial to understand the background NO2. To do this, one quick way is to observe NO2 measurements when the powerplants were out of services or at very low electricity generation. According to the news, [Hurricane Maria](https://www.mercycorps.org/blog/quick-facts-hurricane-maria-puerto-rico) hitted the island and electricity was cut off to 100 percent of the island in end of September 2017. As such, we would take NO2 samples from October to December 2017 and comparing it with that recorded from July to September 2017. It is expected that there is decreasing NO2 trend between the two periods.\n\nFurther, I would take samples before and after the events to see any difference. \n\nThings to do:\n1. Identify any significant difference of NO2 emission over time.\n    - Preprocess images from S5P and GFS\n2. Calculate true NO2 emission generated by power plants\n3. Calculate Emission Factor\n\n\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"list(df_all.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plotting\n#\nimport matplotlib.pyplot as plt\nimport seaborn as sns\ndf_all_date=df_all.set_index('DATE')\nsns.set(rc={'figure.figsize':(11, 4)})\ncols_plot = ['tropospheric_NO2_column_number_density', 'NO2_column_number_density','NO2_slant_column_number_density']\ndf_all_date[cols_plot].plot(linewidth=0.5);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#gldas gives \n\ngldas = ee.ImageCollection('NASA/GLDAS/V021/NOAH/G025/T3H').filterDate(startdate, enddate).filterBounds(boundingBoxes)\n #.filterBounds(ee.Geometry.Point(-66.15, 17.95)) \\\n\n# Get the number of images.\ncount = gldas.size()\nprint('Counting gldas: ', str(count.getInfo())+'\\n')\n\n# Get information about the bands as a list.\nbandName = gldas.first().bandNames()\nprint('Band name: ', bandName.getInfo()); # ee.List of band names","execution_count":null,"outputs":[]},{"metadata":{"id":"YeLxgsiI609A"},"cell_type":"markdown","source":"** Step 6: Explore the weather data data using the Google Earth Engine API**\n1. Try to get the low wind dates according to the bounding boxes.\n2. Check timeseries for 'u_component_of_wind_10m_above_ground' and 'v_component_of_wind_10m_above_ground'.\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\nimage = ee.Image('NOAA/GFS0P25').filterDate(startdate, enddate)\n\n# Create a time series chart.\ntempTimeSeries = ui.Chart.image.seriesByRegion(\n    image, boundingBoxes, ee.Reducer.mean(), 'temperature_2m_above_ground', 200, 'system:time_start', 'label')\n        .setChartType('ScatterChart')\n        .setOptions({\n          title: 'Temperature over time in regions of the bounding boxes',\n          vAxis: {title: 'Temperature (°C)'},\n          lineWidth: 1,\n          pointSize: 4,\n          series: {\n            0: {color: 'FF0000'}, // urban\n            1: {color: '00FF00'}, // forest\n            2: {color: '0000FF'}  // desert\n}});\n\n// Display.\nprint(tempTimeSeries);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gfs = ee.ImageCollection('NOAA/GFS0P25').filterDate(startdate, enddate).filterBounds(boundingBoxes)\n #.filterBounds(ee.Geometry.Point(-66.15, 17.95)) \\\n\n# Get the number of images.\ncount = gfs.size()\nprint('Counting gfs: ', str(count.getInfo())+'\\n')\n\n# Get information about the bands as a list.\nbandName = gfs.first().bandNames()\nprint('Band name(gfs): ', bandName.getInfo()); # ee.List of band names\n\n#gfsv2 = gfs.filter(ee.Filter.lte('v_component_of_wind_10m_above_ground', 2))\n\n#Get statistics for a property of the images in the collection.\nwindStats = gfs.aggregate_stats('v_component_of_wind_10m_above_ground');\nprint('v_component_of_wind_10m_above_ground statistics: ', windStats.getInfo());\n\nrange = gfs.reduceColumns(ee.Reducer.minMax(), ['v_component_of_wind_10m_above_ground'])\nprint('Wind speed range: ', ee.Number(range.get('min')).getInfo(), ee.Number(range.get('max')).getInfo())\ntype(range)\n\n\n# Get the number of images.\n#count = gfsv2.size()\n#print('Counting gfs_v1: ', str(count.getInfo())+'\\n')\n\n##########\n## Want to know the which date having calm wind for getting error-free NO2 observation.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"","execution_count":null}],"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}