{"cells":[{"metadata":{},"cell_type":"markdown","source":"<h1>Questions for DS4G-EIE AMA</h1>\nFirst, some helper code."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport folium\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)\n\ndef add_ee_layer(self, ee_image_object, vis_params, name):\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    ).add_to(self)\n    \n# modified from https://www.kaggle.com/paultimothymooney/how-to-get-started-with-the-earth-engine-data\n# zoom_country=True when showing Puerto Rico\n# takes default lat and long (unless they are overwritten) and zoom=8\ndef plot_ee_data_on_map(add_ee_layer, begin_date, end_date,\n                        min_value, max_value, opacity=1.0, lat=18.233, long=-66.279, \n                        zoom_country=True, res=0.01):\n    \n    product = \"NO2\"\n    dataset = \"COPERNICUS/S5P/OFFL/L3_NO2\"\n    column = 'tropospheric_NO2_column_number_density'\n    \n    if zoom_country: # zoom at the country level, default lat and long\n        zoom_start = 8\n        lat1 = lat-0.33; long1 = long-1.06\n        lat2 = lat+0.33; long2 = long+1.06\n    else: # zoom at a Power Plant level\n        zoom_start = 13\n        lat1 = lat-res/2; long1 = long-res/2\n        lat2 = lat+res/2; long2 = long+res/2\n    rectangle = ee.Geometry.Rectangle([long1, lat1, long2, lat2]) # (x, y) math style   \n        \n    Map = folium.Map(location=[lat, long], zoom_start=zoom_start) # (y, x) geo style\n    folium.Map.add_ee_layer = add_ee_layer\n\n    sat_image = (ee.ImageCollection(dataset)\n           .select(column)\n           .filterDate(begin_date, end_date)\n           .mean()\n          )\n    \n    vis_params = {\n      'min': min_value,\n      'max': max_value,\n      'opacity': opacity,\n      'palette': ['green', 'blue', 'yellow', 'red']}\n        \n    Map.add_ee_layer(sat_image.clip(rectangle), vis_params, product)\n    Map.add_child(folium.LayerControl())\n    display(Map)\n    return sat_image    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<h2>1.- Resolution</h2>\n<h3>1.1 EE Resolution</h3>\nNO2 product has a resolution of 0.01 arc degrees (in other places it is said to be 7 x 7 km, or 5 x 5 for images after oct 2019 which is not our case). Let's see how a 0.01 arc degrees \"pixel\" looks like:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# A \"pixel\" around San Juan Power Plant\nres = 0.01 # pixel resolution in arc degrees\nlong = -66.1045; lat = 18.427 # San Juan\n\nbegin_date = '2019-05-03'; end_date = '2019-05-04'\nmin_value = 0.00001; max_value = 0.000025\n\nsat_image = plot_ee_data_on_map(add_ee_layer, begin_date, end_date, min_value, max_value, \n                        zoom_country=False, long=long, lat=lat, res=res)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check its size with haversine:"},{"metadata":{"trusted":true},"cell_type":"code","source":"res = 0.01\nlat1 = lat-res/2; long1 = long-res/2\nlat2 = lat+res/2; long2 = long+res/2\n\nfrom haversine import haversine # distance in km (default) between points in UTM coordinates\np1 = (lat1, long1) # (y, x) geo style\np2 = (lat1, long2)\np3 = (lat2, long1)\np4 = (lat2, long2)\nprint(\"Horizontal pixel side\", haversine(p1, p2), \"km\")\nprint(\"Vertical pixel side\", haversine(p1, p3), \"km\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<h3>*Question 1.1: Shouldn't the Horizontal and Vertical pixel sides be larger, like 7 x 7 km?*</h3>\n<h2>1.2 TIFF images resolution</h2>\n(from https://www.kaggle.com/paultimothymooney/explore-image-metadata-s5p-gfs-gldas)"},{"metadata":{"trusted":true},"cell_type":"code","source":"import rasterio as rio\nimport os\n\ns5p_file = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20190501T161114_20190507T174400.tif'\ndef preview_meta_data(file_name):\n    with rio.open(file_name) as img_filename:\n        print('Metadata for: ',file_name)\n        print('Bounding Box:',img_filename.bounds)\n        print('Resolution:',img_filename.res)\n        print('Tags:',img_filename.tags())\n        print('More Tags:',img_filename.tags(ns='IMAGE_STRUCTURE'))\n        print('Number of Channels =',img_filename.count,'\\n')\n\npreview_meta_data(s5p_file)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<h3>*Question 1.2: Metadata in the s5p TIFF image says its resolution is 0.004 x 0.004??, \nwhile in EE it was 0.01*</h3>\n<h2>2.- From EE to numpy</h2>\n<h3>Now let's check the shape of the pixel array (should be 1x1)</h3>\n(This is a workaround, I hope there are better ways)\n\nRather than taking the mean, as before, I first want to know how many images are there in a given time period and have contributed to the mean."},{"metadata":{"trusted":true},"cell_type":"code","source":"rectangle = ee.Geometry.Rectangle([long1, lat1, long2, lat2]) # the pixel\n\ncollection = (ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')\n  .filterDate(begin_date, end_date))\n\ncount = collection.size()\nprint('Count: ', str(count.getInfo())+'\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It looks there are 14 images in just one day period. Let's get the first image."},{"metadata":{"trusted":true},"cell_type":"code","source":"image = collection.first()\nband_arrs = image.sampleRectangle(rectangle);\n# Get individual band arrays.\nband_arr = band_arrs.get('tropospheric_NO2_column_number_density')\nnp_arr = np.array(band_arr.getInfo()) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<h3>Big crash!</h3>\nI took a look at the metadata in the EE code editor but didn't see anything wrong in the 14 Images (features) in ImageCollection.\n\nThen I developed a workaround to check which images were REAL images: looping through the orbits until I find an image which does not throw an exception when converting to numpy."},{"metadata":{"trusted":true},"cell_type":"code","source":"orbitStats = collection.aggregate_stats(\"ORBIT\")\nminOrbit = orbitStats.getInfo()['values']['min']\nmaxOrbit = orbitStats.getInfo()['values']['max']\narrayList = []\nfor orbit in range(minOrbit, maxOrbit+1):\n    index = orbit - minOrbit\n    filtered = collection.filterMetadata('ORBIT', 'equals', orbit);\n    image = filtered.first()\n    try:\n        date = image.date()\n        # only arrives here in case of no error: incomplete image throws exception\n        band_arrs = image.sampleRectangle(rectangle);\n        # Get individual band arrays.\n        band_arr = band_arrs.get('tropospheric_NO2_column_number_density')\n\n        # Transfer the arrays from server to client and cast as np array.\n        try:\n            np_arr = np.array(band_arr.getInfo()) \n            # only arrives here in case of no error\n            arrayList.append(np_arr)\n            print(\"Catch!\", date.format().getInfo())\n        except ee.EEException:\n            print(\"Bad luck\", date.format().getInfo())\n    except ee.EEException:\n            # incomplete image\n            print(\"Very bad luck, incomplete image\", orbit)\n        \nprint(\"Fin\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<h3>*Question 2.1: Is there another way to catch the non-empty images in an ImageCollection than this orbit workaround?*</h3>\nAmong the 14 \"images\", just one turned to have data. Let's check the size is 1x1"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(arrayList[0].shape, arrayList[0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<h3>*Question 2.2: Why the size of the array extracted from the pixel Image has a (2, 2) shape rather than (1, 1) if it's supposed to have the minimum resolution (0.01 arc degree)?*</h3>"}],"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}