{"cells":[{"metadata":{},"cell_type":"markdown","source":"**Update information**\n\nUpdate1(Commit2): Fixed mask displacement in map.\n\nUpdate2(Commit4): Add plots and summary statistics of data.\n\nUpdate3(Commit5): Fixed distplot in create_mask(img, df, MASK_SIZE).  values_in_mask -> values_in_mask[:400]\n\nUpdate4(Commit7): Change mask shape　square to circular.\n\nUpdate5(Commit8): Add time series analysis.\n\nUpdate6(Commit9): Add Simplified Emissions Factor calculation."},{"metadata":{},"cell_type":"markdown","source":"# About this kernel\n\nI'm interested in how to focus on specific power plants and want to dvance our analysis.\n\nI think, I can calculate more robust EF if I focus on specific region near each power plants.\n\nIn this kernel, I tried to get data only around specific power plants and visualize.\n\nIf we can do so, we can get more detail infomation."},{"metadata":{},"cell_type":"markdown","source":"# Preperation\n\nFor analisys, load library, snnipets and data."},{"metadata":{},"cell_type":"markdown","source":"## Load library"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"from decimal import Decimal, ROUND_HALF_UP, ROUND_HALF_EVEN\nimport folium\nimport math\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport os\nimport pandas as pd\nimport rasterio as rio\nimport seaborn as sns\nimport tifffile as tiff ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Snnipets"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def overlay_image_on_puerto_rico_with_Marker(df,file_name,band_layer,lat,lon,zoom):\n    \"\"\"\n    Visualize map overlayed data & plant markers.\n    \"\"\"\n    band = rio.open(file_name).read(band_layer)\n    m = folium.Map([lat, lon], zoom_start=zoom)\n    color={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\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    \n    for i in range(0,len(df)):\n        popup = folium.Popup(str(df.primary_fuel[i:i+1]))\n        folium.Marker([df[\"latitude\"].iloc[i],df[\"longitude\"].iloc[i]],\n                     icon=folium.Icon(icon_color='red',icon ='bolt',prefix='fa',color=color[df.primary_fuel.iloc[i]])).add_to(m)\n        \n    return m","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def overlay_image_on_puerto_rico_df(df, img, zoom):\n    lat_map=df.iloc[[0]].loc[:,[\"latitude\"]].iat[0,0]\n    lon_map=df.iloc[[0]].loc[:,[\"longitude\"]].iat[0,0]\n    m = folium.Map([lat_map, lon_map], zoom_start=zoom)\n    color={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\n    folium.raster_layers.ImageOverlay(\n        image=img,\n        bounds = [[18.56,-67.32,],[17.90,-65.194]],\n        colormap=lambda x: (1, 0, 0, x),\n    ).add_to(m)\n    \n    for i in range(0,len(df)):\n        popup = folium.Popup(str(df.primary_fuel[i:i+1]))\n        folium.Marker([df[\"latitude\"].iloc[i],df[\"longitude\"].iloc[i]],\n                     icon=folium.Icon(icon_color='red',icon ='bolt',prefix='fa',color=color[df.primary_fuel.iloc[i]])).add_to(m)\n        \n    return m","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def 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    color={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\n    plot = folium.Map(location=location,zoom_start=zoom)\n    for i in range(0,len(df)):\n        popup = folium.Popup(str(df.primary_fuel[i:i+1]))\n        folium.Marker([df[latitude_column].iloc[i],df[longitude_column].iloc[i]],popup=popup,\n                     icon=folium.Icon(icon_color='red',icon ='bolt',prefix='fa',color=color[df.primary_fuel.iloc[i]])).add_to(plot)\n    return(plot)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def 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":"## Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants = pd.read_csv('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\npower_plants = split_column_into_new_columns(power_plants,'.geo','latitude',50,66)\npower_plants = split_column_into_new_columns(power_plants,'.geo','longitude',31,48)\npower_plants['latitude'] = power_plants['latitude'].astype(float)\na = np.array(power_plants['latitude'].values.tolist()) # 18 instead of 8\npower_plants['latitude'] = np.where(a < 10, a+10, a).tolist() \npower_plants_df = power_plants.sort_values('capacity_mw',ascending=False).reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_df.columns","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Overview of whole data of Puerto Rico\n\nThere are 35 plants in Puerto Rico, and we can use NO2 gas data as starter data."},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_df[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh']].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image ='/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180701T161259_20180707T175356.tif'\nlatitude=18.1429005246921; longitude=-65.4440010699994\noverlay_image_on_puerto_rico_with_Marker(power_plants_df,image,band_layer=1,lat=latitude,lon=longitude,zoom=8)\n\n#https://www.kaggle.com/paultimothymooney/explore-image-metadata-s5p-gfs-gldas\n#band1: NO2_column_number_density\n#band2: tropospheric_NO2_column_number_density\n#band3: stratospheric_NO2_column_number_density\n#band4: NO2_slant_column_number_density","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img = tiff.imread(image)[:,:,0]\nprint(\"mean: \", img.mean())\nprint(\"max: \", img.max())\nprint(\"min: \", img.min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(img.reshape([70300,]), kde=False, rug=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Focus on Costa Sur power plants\n\nNow I'll focus on specific power plants. If we find way to analyze each power plants, we may get more important insight."},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plant_name = \"Costa Sur\"\npower_plants_specific = power_plants_df[power_plants_df[\"name\"] == power_plant_name]\npower_plants_specific","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lat_map=power_plants_specific.iloc[[0]].loc[:,[\"latitude\"]].iat[0,0]\nlon_map=power_plants_specific.iloc[[0]].loc[:,[\"longitude\"]].iat[0,0]\nplot_points_on_map(power_plants_df[power_plants_df[\"name\"] == power_plant_name],0,425,'latitude',lat_map,'longitude',lon_map,13)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create Mask\n\nGiven tiff data are widespread, so I try to extract data near focused power plant. \n\nWe thought that this operation could be realized by replacing the data other than around the focused power plant with zero.\n\nIn order to do this operation, \n\n・1st: Create zero matrix which size is same as the whole data.\n\n・2nd： Create implicit mask ( this express \"region around the focused pwoer plant\") as a two-dimensional range.\n\n・3nd: Replace elements with orginal data if the element in mask.\n\nIn this strategy we have to define MASK_SIZE for create mask. For now I do this by eye.\n\nIn \"MASK_SIZE\", 1 corresponds to 0.0045 in longitude or latitude. I calculated this value for the following reasons. Original data is in bounds [[18.6,-67.3,],[17.9,-65.2]], and the shape is (148, 475). "},{"metadata":{"trusted":true},"cell_type":"code","source":"MASK_SIZE = 10","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = tiff.imread('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180701T161259_20180707T175356.tif')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_mask(img, df, MASK_SIZE):\n    \"\"\"\n    input: \n       img: orginal image \n       df: dataframe of focused power plants\n       MASK_SIZE\n    output:\n       numpyarray same size as img.\n    ----------\n    ・1st: Create zero matrix which size is same as the whole data.\n    ・2nd： Create implicit mask ( this express \"region around the focused plant\") as a two-dimensional range for each power plants.\n    ・3nd: Replace elements with orginal data if the element in mask.\n    \"\"\"\n    \n    mask = np.zeros((img.shape[0], img.shape[1]))\n    values_in_mask = np.array([])\n    \n    for i in range(len(df)):\n        lat = float(df.iloc[[i]].loc[:,[\"latitude\"]].iat[0,0])\n        lon = float(df.iloc[[i]].loc[:,[\"longitude\"]].iat[0,0])\n\n    \n        f_lat = (lat - 17.903121359128956)*img.shape[0]/(18.563112930177304 - 17.903121359128956)\n        f_lon = (lon + 67.32297404549217)*img.shape[1]/(-65.19437297127342 + 67.32297404549217)\n        f_lat_int = int(Decimal(str(f_lat -1)).quantize(Decimal('0'), rounding=ROUND_HALF_UP))\n        f_lon_int = int(Decimal(str(f_lon -1)).quantize(Decimal('0'), rounding=ROUND_HALF_UP))\n        \n        mask_lat_min = 148  - f_lat_int - MASK_SIZE\n        mask_lat_max = 148  - f_lat_int + MASK_SIZE\n        mask_lon_min = f_lon_int - MASK_SIZE \n        mask_lon_max = f_lon_int + MASK_SIZE \n        #print(mask_lat_min, mask_lat_max, mask_lon_min, mask_lon_max, f_lat, f_lon)\n        \n        for i in range(img.shape[0]):\n            for j in range(img.shape[1]):\n                if math.sqrt( (i- (147 - f_lat_int))**2 + (j  - f_lon_int )**2) <= MASK_SIZE:\n                    \n                    if image[:,:,0][i][j] < 0:\n                        continue\n                        \n                    mask[i,j] = image[:,:,0][i][j]\n                    values_in_mask = np.append(values_in_mask, mask[i,j])\n                #if i in range(mask_lat_min,mask_lat_max) and j in range(mask_lon_min,mask_lon_max):\n                #    mask[i,j] = image[:,:,0][i][j]\n                #    values_in_mask = np.append(values_in_mask, mask[i,j])\n\n    #print(len(values_in_mask))\n    print(power_plant_name, \" is:\")\n    sns.distplot(values_in_mask[:400], kde=False, rug=True)\n\n    return mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_specific = power_plants_df[power_plants_df[\"name\"] == power_plant_name]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_specific = power_plants_specific.reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mask = create_mask(image, power_plants_specific, MASK_SIZE)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"overlay_image_on_puerto_rico_df(power_plants_specific, mask,zoom=12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(\"mean: \", mask.mean())\nprint(\"max: \", mask.max())\nprint(\"min: \", mask.min())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Focus on other power plants"},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plant_name = \"San Juan CC\"\npower_plants_specific = power_plants_df[power_plants_df[\"name\"] == power_plant_name]\npower_plants_specific = power_plants_specific.reset_index(drop=True)\nmask = create_mask(image, power_plants_df[power_plants_df[\"name\"] == power_plant_name], MASK_SIZE)\noverlay_image_on_puerto_rico_df(power_plants_df[power_plants_df[\"name\"] == power_plant_name], mask,zoom=12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(\"mean: \", mask.mean())\nprint(\"max: \", mask.max())\nprint(\"min: \", mask.min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plant_name = \"Palo Seco\"\npower_plants_specific = power_plants_df[power_plants_df[\"name\"] == power_plant_name]\npower_plants_specific = power_plants_specific.reset_index(drop=True)\nmask = create_mask(image, power_plants_df[power_plants_df[\"name\"] == power_plant_name], MASK_SIZE)\noverlay_image_on_puerto_rico_df(power_plants_specific, mask, zoom=12)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(\"mean: \", mask.mean())\nprint(\"max: \", mask.max())\nprint(\"min: \", mask.min())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Apply to Time Series\n\nNow we can extract data around specific power plants, I apply this method to time series data.\n\nWe have 387 sample data of Sentinel-5P OFFL NO2.\n"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"import glob\nimport re\nimport os\nfrom datetime import datetime, timedelta\n\ns5p_no2_timeseries = glob.glob('../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/*')\nprint(\"There are \", len(s5p_no2_timeseries), \"s5p_no2 time series data.\")\nprint(\"  --------------------  \")\nprint(s5p_no2_timeseries)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Referring to https://www.kaggle.com/gpoulain/eda-ef-with-n2o-time-series-earth-engine, I try to plot specific power plant data with mask.\n\nI try to calculate NO2_column_number_density of whole Puerto Rico near each power plants. "},{"metadata":{"trusted":true},"cell_type":"code","source":"dates = [datetime.strptime(i[:79], '../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_%Y%m%d') for i in  s5p_no2_timeseries]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#just stopped plotting graph version.\ndef create_mask_without_plot(img, df, MASK_SIZE):\n    \"\"\"\n    input: \n       img: orginal image \n       df: dataframe of focused power plants\n       MASK_SIZE\n    output:\n       numpyarray same size as img.\n    ----------\n    ・1st: Create zero matrix which size is same as the whole data.\n    ・2nd： Create implicit mask ( this express \"region around the focused plant\") as a two-dimensional range for each power plants.\n    ・3nd: Replace elements with orginal data if the element in mask.\n    \"\"\"\n    \n    mask = np.zeros((img.shape[0], img.shape[1]))\n    values_in_mask = np.array([])\n    \n    for i in range(len(df)):\n        lat = float(df.iloc[[i]].loc[:,[\"latitude\"]].iat[0,0])\n        lon = float(df.iloc[[i]].loc[:,[\"longitude\"]].iat[0,0])\n\n    \n        f_lat = (lat - 17.903121359128956)*img.shape[0]/(18.563112930177304 - 17.903121359128956)\n        f_lon = (lon + 67.32297404549217)*img.shape[1]/(-65.19437297127342 + 67.32297404549217)\n        f_lat_int = int(Decimal(str(f_lat -1)).quantize(Decimal('0'), rounding=ROUND_HALF_UP))\n        f_lon_int = int(Decimal(str(f_lon -1)).quantize(Decimal('0'), rounding=ROUND_HALF_UP))\n        \n        mask_lat_min = 148  - f_lat_int - MASK_SIZE\n        mask_lat_max = 148  - f_lat_int + MASK_SIZE\n        mask_lon_min = f_lon_int - MASK_SIZE \n        mask_lon_max = f_lon_int + MASK_SIZE \n        #print(mask_lat_min, mask_lat_max, mask_lon_min, mask_lon_max, f_lat, f_lon)\n        \n        for i in range(img.shape[0]):\n            for j in range(img.shape[1]):\n                if math.sqrt( (i- (147 - f_lat_int))**2 + (j  - f_lon_int )**2) <= MASK_SIZE:\n                    \n                    if image[:,:,0][i][j] < 0:\n                        continue\n\n                    mask[i,j] = image[:,:,0][i][j]\n                    values_in_mask = np.append(values_in_mask, mask[i,j])\n                #if i in range(mask_lat_min,mask_lat_max) and j in range(mask_lon_min,mask_lon_max):\n                #    mask[i,j] = image[:,:,0][i][j]\n                #    values_in_mask = np.append(values_in_mask, mask[i,j])\n\n    #print(len(values_in_mask))\n    #print(power_plant_name, \" is:\")\n    #sns.distplot(values_in_mask[:400], kde=False, rug=True)\n\n    return mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stat_max = np.array([])\nstat_mean = np.array([])\nstat_min = np.array([])\n\nfor name in set(power_plants_df[\"name\"]):\n    power_plant_name = name\n    power_plants_specific = power_plants_df[power_plants_df[\"name\"] == power_plant_name]\n    power_plants_specific = power_plants_specific.reset_index(drop=True)\n\n    stat_max_tmp = []\n    stat_mean_tmp = []\n    stat_min_tmp = []\n\n    for i in s5p_no2_timeseries:  \n        image = tiff.imread(i)\n        mask = create_mask_without_plot(image, power_plants_df[power_plants_df[\"name\"] == power_plant_name], MASK_SIZE)\n        stat_max_tmp.append(mask.max())\n        stat_mean_tmp.append(mask.mean())\n        stat_min_tmp.append(mask.min())\n        \n    stat_max = np.append(stat_max, stat_max_tmp)\n    stat_mean = np.append(stat_mean, stat_mean_tmp)\n    stat_min = np.append(stat_min, stat_min_tmp)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stat_max_reshape = stat_max.reshape(32, int(len(stat_max)/32))\nstat_mean_reshape = stat_mean.reshape(32, int(len(stat_mean)/32))\nstat_min_reshape = stat_min.reshape(32,int(len(stat_min)/32))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stat_max_nanmean = np.nanmean(stat_max_reshape, axis=0)\nstat_mean_nanmean = np.nanmean(stat_mean_reshape, axis=0)\nstat_min_nanmean = np.nanmean(stat_min_reshape, axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#results = pd.DataFrame(index=dates, data=stat[:387], columns=['San Juan CC'])\nresults = pd.DataFrame(columns=['max', 'mean', 'min'])\nresults = pd.DataFrame({'max': stat_max_nanmean,\n                    'mean': stat_mean_nanmean,\n                    'min': stat_min_nanmean},\n                    index=dates)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results[['max','min']].plot()\nplt.title('Max of NO2_column_number_density in Puerto Rico')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results[['mean','min']].plot()\nplt.title('Mean of NO2_column_number_density in Puerto Rico')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Calculate Simplified Emissions Factor"},{"metadata":{},"cell_type":"markdown","source":"Refer from https://www.kaggle.com/paultimothymooney/overview-of-the-eie-analytics-challenge#E-=-A-x-EF-x-(1-ER/100),\nwe can calculate simplified emissions factor by below formula:\n\nSimplified Emissions Factor = (Measure of NO2 emissions) / (Quanity of electricity generated)\n\n\nUnit of NO2_column_number_density is mol/m^2, so I think I should multiply the area by the NO2_column_number_density. But for simplicity, I use NO2_column_number_density as is."},{"metadata":{"trusted":true},"cell_type":"code","source":"#power_plant_name = \"San Juan CC\"\n#power_plants_specific = power_plants_df[power_plants_df[\"name\"] == power_plant_name]\n#Quanity_of_electricity_generated = power_plants_specific.loc[:,[\"estimated_generation_gwh\"]].iat[0,0]\nQuanity_of_electricity_generated = np.sum(power_plants_df.loc[:,[\"estimated_generation_gwh\"]])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"(results[['mean','min']]/Quanity_of_electricity_generated[0]).plot()\nplt.title('Mean Simplified Emissions Factor in Puerto Rico')\n#plt.title('Mean Simplified Emissions Factor in San Juan CC')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Simplified_Emissions_Factor = (results.loc[:,['mean']]/Quanity_of_electricity_generated[0]).mean().at['mean']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f\"Simplified Emissions Factor of {power_plant_name}:\",Simplified_Emissions_Factor,\"mol/gwh・m^2\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Can we extract emissions by mask?"},{"metadata":{},"cell_type":"markdown","source":"Simply thinking, the NO2 concentration around the energy plant seems to be the highest. And by determining the threshold and looking at the concentration in the certain region, you should be able to determine the NO2 emissions.\n\n\nBy calculating the average value at each data point and looking at their distribution, I tried to see if it would be better to process them into masks rather than raw."},{"metadata":{"trusted":true},"cell_type":"code","source":"images_avg =  np.array([])\nfor i in s5p_no2_timeseries:\n    image = tiff.imread(i)\n    \n    images_avg = np.append(np.nanmean(image[:,:,0]), images_avg)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"masks_nonzero_avg = np.array([])\n\nfor name in set(power_plants_df[\"name\"]):\n    power_plant_name = name\n    power_plants_specific = power_plants_df[power_plants_df[\"name\"] == power_plant_name]\n    power_plants_specific = power_plants_specific.reset_index(drop=True)\n\n\n    for i in s5p_no2_timeseries:  \n        image = tiff.imread(i)\n        mask = create_mask_without_plot(image, power_plants_df[power_plants_df[\"name\"] == power_plant_name], MASK_SIZE)\n        mask_nonzero = [i for i in mask.reshape([70300,]) if i > 0]\n        masks_nonzero_avg = np.append(masks_nonzero_avg, np.nanmean(mask_nonzero))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(images_avg, kde=False, rug=True)\nprint(\"s5p_no2 data mean:\", np.nanmean(images_avg) )\nprint(\"s5p_no2 data variable:\", np.nanvar(images_avg) )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(masks_nonzero_avg, kde=False, rug=True)\nprint(\"masked s5p_no2 data mean:\", np.nanmean(masks_nonzero_avg) )\nprint(\"masked s5p_no2 data variable:\", np.nanvar(masks_nonzero_avg) )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With mask, distribution seems simple gaussian.　We might say this distribution is as we expected.\n\nBut if we calculate simple average of whole data, distribution seems have two peeks.\n\nThe peek near 0.000055 mol/m^2 may be same, but what is one near 0.000045 mol/m^2 ?\n\n"},{"metadata":{},"cell_type":"markdown","source":"Thank you for your reading!"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}