{ "cells": [ { "cell_type": "markdown", "metadata": { "_cell_guid": "2695f074-c7f5-b5a4-8f6a-69e3b1082677" }, "source": [ "## Panchromatic sharpening ##" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "18bd5886-c61c-0c56-dbbb-18b4b2da2381" }, "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" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "ca94efd9-0edc-65aa-9af6-db20387c2751" }, "outputs": [], "source": [ "# libs\n", "import numpy as np\n", "import tifffile as tiff\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "from skimage.transform import rescale\n", "import skimage.color as color" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "201fc77d-30e3-f619-5e64-646ce90131fd" }, "outputs": [], "source": [ "# read m and p data\n", "IM_ID = '6120_2_2'\n", "m = tiff.imread('../input/sixteen_band/{}_M.tif'.format(IM_ID))\n", "pan = tiff.imread('../input/sixteen_band/{}_P.tif'.format(IM_ID))\n", "\n", "print ('m_shape: ', m.shape)\n", "print ('pan_shape: ', pan.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "1cba3813-66a6-f783-e4dd-9cf12b20325c" }, "outputs": [], "source": [ "# functions \n", "def 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[tb] = b\n", " out[:,:,i] =t\n", " return out.astype(np.uint8)\n", "\n", "def 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" }, "outputs": [], "source": [ "rgbn, sharpened_browley_simple, I = pansharpen(m, pan, method='simple_browley', all_data=True)\n", "sharpened_mean = pansharpen(m, pan, method='sample_mean')\n", "sharpened_esri = pansharpen(m, pan, method='esri')\n", "\n", "sharpened_browley2 = pansharpen(m, pan, W=0.1)\n", "sharpened_browley4 = pansharpen(m, pan, W=0.2)\n", "sharpened_browley6 = pansharpen(m, pan, W=0.3)\n", "sharpened_browley8 = pansharpen(m, pan, W=0.4)\n", "\n", "sharpened_ihs2 = pansharpen(m, pan, method='hsv', W=0.1)\n", "sharpened_ihs4 = pansharpen(m, pan, method='hsv', W=0.2)\n", "sharpened_ihs6 = pansharpen(m, pan, method='hsv', W=0.3)\n", "sharpened_ihs8 = pansharpen(m, pan, method='hsv', W=0.4)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "d0829b56-84ab-1cbe-3cc8-aa69dbb5231a" }, "outputs": [], "source": [ "fig, axes = plt.subplots(ncols=3, nrows=4, figsize=(20, 20))\n", "ax = axes.ravel()\n", "\n", "# no IR sharpening\n", "ax[0].imshow(stretch(rgbn[500:1000,500:1000,:3]))\n", "ax[0].set_title('RGB')\n", "ax[0].axis('off')\n", "ax[1].imshow(stretch(sharpened_mean[500:1000,500:1000,:]))\n", "ax[1].set_title('Sharpened Mean') \n", "ax[1].axis('off')\n", "ax[2].imshow(stretch(sharpened_browley_simple[500:1000,500:1000,:]))\n", "ax[2].set_title('Sharpened Browley Simple')\n", "ax[2].axis('off')\n", "ax[3].imshow(stretch(sharpened_esri[500:1000,500:1000,:3]))\n", "ax[3].set_title('Sharpened ESRI')\n", "ax[3].axis('off')\n", "\n", "# IR Browley\n", "ax[4].imshow(stretch(sharpened_browley2[500:1000,500:1000,:3]))\n", "ax[4].set_title('Sharpened Browley(0.1)')\n", "ax[4].axis('off')\n", "ax[5].imshow(stretch(sharpened_browley4[500:1000,500:1000,:3]))\n", "ax[5].set_title('Sharpened Browley(0.2)')\n", "ax[5].axis('off')\n", "ax[6].imshow(stretch(sharpened_browley6[500:1000,500:1000,:3]))\n", "ax[6].set_title('Sharpened Browley(0.3)')\n", "ax[6].axis('off')\n", "ax[7].imshow(stretch(sharpened_browley8[500:1000,500:1000,:3]))\n", "ax[7].set_title('Sharpened Browley(0.4)')\n", "ax[7].axis('off')\n", "\n", "# IR HSV\n", "ax[8].imshow(stretch(sharpened_ihs2[500:1000,500:1000,:3]))\n", "ax[8].set_title('Sharpened IHS(0.1)')\n", "ax[8].axis('off')\n", "ax[9].imshow(stretch(sharpened_ihs4[500:1000,500:1000,:3]))\n", "ax[9].set_title('Sharpened IHS(0.2)')\n", "ax[9].axis('off')\n", "ax[10].imshow(stretch(sharpened_ihs6[500:1000,500:1000,:3]))\n", "ax[10].set_title('Sharpened IHS(0.3)')\n", "ax[10].axis('off')\n", "ax[11].imshow(stretch(sharpened_ihs8[500:1000,500:1000,:3]))\n", "ax[11].set_title('Sharpened IHS(0.4)')\n", "ax[11].axis('off')\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "8aedf050-7792-d53e-891e-9d8ead419d21" }, "outputs": [], "source": "" } ], "metadata": { "_change_revision": 0, "_is_fork": false, "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 }