{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.0"}},"nbformat":4,"nbformat_minor":0,"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"2695f074-c7f5-b5a4-8f6a-69e3b1082677","_active":false},"source":"## Panchromatic sharpening  ##","outputs":[]},{"cell_type":"markdown","metadata":{"_cell_guid":"18bd5886-c61c-0c56-dbbb-18b4b2da2381","_active":false},"source":"Info: [Fundamentals of panchromatic sharpening][1]\n\n\n  [1]: http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/fundamentals-of-panchromatic-sharpening.htm","outputs":[]},{"cell_type":"code","execution_count":1,"metadata":{"_cell_guid":"ca94efd9-0edc-65aa-9af6-db20387c2751","_active":false},"outputs":[],"source":"# libs\nimport numpy as np\nimport tifffile as tiff\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom skimage.transform import rescale\nimport skimage.color as color","execution_state":"idle"},{"cell_type":"code","execution_count":2,"metadata":{"_cell_guid":"201fc77d-30e3-f619-5e64-646ce90131fd","_active":false,"collapsed":false},"outputs":[],"source":"# read m and p data\nIM_ID = '6010_1_2'\nm = tiff.imread('../input/sixteen_band/{}_M.tif'.format(IM_ID))\npan = tiff.imread('../input/sixteen_band/{}_P.tif'.format(IM_ID))\n\nprint ('m_shape: ', m.shape)\nprint ('pan_shape: ', pan.shape)","execution_state":"idle"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"1cba3813-66a6-f783-e4dd-9cf12b20325c","_active":false},"outputs":[],"source":"# functions \ndef stretch(bands, lower_percent=2, higher_percent=98):\n    out = np.zeros_like(bands)\n    for i in range(3):\n        a = 0 \n        b = 255 \n        c = np.percentile(bands[:,:,i], lower_percent)\n        d = np.percentile(bands[:,:,i], higher_percent)        \n        t = a + (bands[:,:,i] - c) * (b - a) / (d - c)    \n        t[t<a] = a\n        t[t>b] = b\n        out[:,:,i] =t\n    return out.astype(np.uint8)\n\ndef pansharpen(m, pan, method='browley', W=0.1, all_data=False):\n    \n    # get m_bands\n    rgbn = np.empty((m.shape[1], m.shape[2], 4)) \n    rgbn[:,:,0] = m[4,:,:] # red\n    rgbn[:,:,1] = m[2,:,:] # green\n    rgbn[:,:,2] = m[1,:,:] # blue\n    rgbn[:,:,3] = m[6,:,:] # NIR-1\n    \n    # scaled them\n    rgbn_scaled = np.empty((m.shape[1]*4, m.shape[2]*4, 4))\n\n    for i in range(4):\n        img = rgbn[:,:,i]\n        scaled = rescale(img, (4,4))\n        rgbn_scaled[:,:,i] = scaled\n    \n    # check size and crop for pan band\n    if pan.shape[0] < rgbn_scaled.shape[0]:\n        rgbn_scaled = rgbn_scaled[:pan.shape[0],:, :]\n    else:\n        pan = pan[:rgbn_scaled.shape[0], :]\n        \n    if pan.shape[1] < rgbn_scaled.shape[1]:\n        rgbn_scaled = rgbn_scaled[:,:pan.shape[1], :]\n    else:\n        pan = pan[:,:rgbn_scaled.shape[1]]\n\n    R = rgbn_scaled[:,:,0]\n    G = rgbn_scaled[:,:,1]\n    B = rgbn_scaled[:,:,2]\n    I = rgbn_scaled[:,:,3]\n\n    image = None\n\n    if method == 'simple_browley':\n        all_in = R+G+B\n        prod = np.multiply(all_in, pan)\n\n        r = np.multiply(R, pan/all_in)[:, :, np.newaxis]\n        g = np.multiply(G, pan/all_in)[:, :, np.newaxis]\n        b = np.multiply(B, pan/all_in)[:, :, np.newaxis]\n\n        image = np.concatenate([r,g,b], axis=2)\n    \n    if method == 'sample_mean':\n        r = 0.5 * (R + pan)[:, :, np.newaxis]\n        g = 0.5 * (G + pan)[:, :, np.newaxis]\n        b = 0.5 * (B + pan)[:, :, np.newaxis]\n        \n        image = np.concatenate([r,g,b], axis=2)\n        \n    if method == 'esri':\n        ADJ = pan-rgbn_scaled.mean(axis=2) \n        r = (R + ADJ)[:, :, np.newaxis]\n        g = (G + ADJ)[:, :, np.newaxis]\n        b = (B + ADJ)[:, :, np.newaxis]\n        i = (I + ADJ)[:, :, np.newaxis]\n        \n        image = np.concatenate([r,g,b,i], axis=2)\n    \n    if method == 'browley':\n        DNF = (pan - W*I)/(W*R+W*G+W*B)\n\n        r = (R * DNF)[:, :, np.newaxis]\n        g = (G * DNF)[:, :, np.newaxis]\n        b = (B * DNF)[:, :, np.newaxis]\n        i = (I * DNF)[:, :, np.newaxis]\n\n        image = np.concatenate([r,g,b,i], axis=2)\n    \n    if method == 'hsv':\n        hsv = color.rgb2hsv(rgbn_scaled[:,:,:3])\n        hsv[:,:,2] = pan - I*W\n        image = color.hsv2rgb(hsv)\n        \n    if all_data:\n        return rgbn_scaled, image, I\n    else:\n        return image"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"11c3973d-0fcc-7055-c61b-c792ff186b13","_active":false},"outputs":[],"source":"rgbn, sharpened_browley_simple, I = pansharpen(m, pan, method='simple_browley', all_data=True)\nsharpened_mean = pansharpen(m, pan, method='sample_mean')\nsharpened_esri = pansharpen(m, pan, method='esri')\n\nsharpened_browley2 = pansharpen(m, pan, W=0.1)\nsharpened_browley4 = pansharpen(m, pan, W=0.2)\nsharpened_browley6 = pansharpen(m, pan, W=0.3)\nsharpened_browley8 = pansharpen(m, pan, W=0.4)\n\nsharpened_ihs2 = pansharpen(m, pan, method='hsv', W=0.1)\nsharpened_ihs4 = pansharpen(m, pan, method='hsv', W=0.2)\nsharpened_ihs6 = pansharpen(m, pan, method='hsv', W=0.3)\nsharpened_ihs8 = pansharpen(m, pan, method='hsv', W=0.4)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"d0829b56-84ab-1cbe-3cc8-aa69dbb5231a","_active":false},"outputs":[],"source":"fig, axes = plt.subplots(ncols=3, nrows=4, figsize=(20, 20))\nax = axes.ravel()\n\n# no IR sharpening\nax[0].imshow(stretch(rgbn[500:1000,500:1000,:3]))\nax[0].set_title('RGB')\nax[0].axis('off')\nax[1].imshow(stretch(sharpened_mean[500:1000,500:1000,:]))\nax[1].set_title('Sharpened Mean')    \nax[1].axis('off')\nax[2].imshow(stretch(sharpened_browley_simple[500:1000,500:1000,:]))\nax[2].set_title('Sharpened Browley Simple')\nax[2].axis('off')\nax[3].imshow(stretch(sharpened_esri[500:1000,500:1000,:3]))\nax[3].set_title('Sharpened ESRI')\nax[3].axis('off')\n\n# IR Browley\nax[4].imshow(stretch(sharpened_browley2[500:1000,500:1000,:3]))\nax[4].set_title('Sharpened Browley(0.1)')\nax[4].axis('off')\nax[5].imshow(stretch(sharpened_browley4[500:1000,500:1000,:3]))\nax[5].set_title('Sharpened Browley(0.2)')\nax[5].axis('off')\nax[6].imshow(stretch(sharpened_browley6[500:1000,500:1000,:3]))\nax[6].set_title('Sharpened Browley(0.3)')\nax[6].axis('off')\nax[7].imshow(stretch(sharpened_browley8[500:1000,500:1000,:3]))\nax[7].set_title('Sharpened Browley(0.4)')\nax[7].axis('off')\n\n# IR HSV\nax[8].imshow(stretch(sharpened_ihs2[500:1000,500:1000,:3]))\nax[8].set_title('Sharpened IHS(0.1)')\nax[8].axis('off')\nax[9].imshow(stretch(sharpened_ihs4[500:1000,500:1000,:3]))\nax[9].set_title('Sharpened IHS(0.2)')\nax[9].axis('off')\nax[10].imshow(stretch(sharpened_ihs6[500:1000,500:1000,:3]))\nax[10].set_title('Sharpened IHS(0.3)')\nax[10].axis('off')\nax[11].imshow(stretch(sharpened_ihs8[500:1000,500:1000,:3]))\nax[11].set_title('Sharpened IHS(0.4)')\nax[11].axis('off')\n\nplt.tight_layout()\nplt.show()"}]}