{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Tutorial 1\n### Transforming Multi-Spectral and Panchromatic images in a series of patches\n\nThis tutorial has the goal of generating a series of patches starting from the raw data provided by a satellite.\nDuring the excercise, two main tools for the creation of patches will be presented: Gdal, a library for the management of GeoTIFF data, and scipy.io, a sub-library that allow the usage of MATLAB files in python.","metadata":{"pycharm":{"name":"#%% md\n"}}},{"cell_type":"markdown","source":"### Importing Libraries","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\nfrom osgeo import gdal\nfrom scipy import io\nfrom matplotlib import pyplot as plt\nimport time\nimport torch\nfrom torch import nn\nimport math\nfrom skimage import transform\nimport imageio\n\n%matplotlib inline","metadata":{"jupyter":{"outputs_hidden":false},"pycharm":{"name":"#%%\n"},"collapsed":false,"execution":{"iopub.status.busy":"2022-10-05T19:03:19.980705Z","iopub.execute_input":"2022-10-05T19:03:19.981072Z","iopub.status.idle":"2022-10-05T19:03:20.771334Z","shell.execute_reply.started":"2022-10-05T19:03:19.981042Z","shell.execute_reply":"2022-10-05T19:03:20.770405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\nsys.path.insert(1, '../input/iadf-school-ds-1')\n","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.773589Z","iopub.execute_input":"2022-10-05T19:03:20.773888Z","iopub.status.idle":"2022-10-05T19:03:20.777959Z","shell.execute_reply.started":"2022-10-05T19:03:20.773858Z","shell.execute_reply":"2022-10-05T19:03:20.777016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import spectral_utils as ut","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.779303Z","iopub.execute_input":"2022-10-05T19:03:20.779611Z","iopub.status.idle":"2022-10-05T19:03:20.793375Z","shell.execute_reply.started":"2022-10-05T19:03:20.779580Z","shell.execute_reply":"2022-10-05T19:03:20.792273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Exploring GDAL Dataset object","metadata":{}},{"cell_type":"code","source":"gdal.AllRegister() # It is a recommended to use this function before opening any data. \n                   # It load in memory all known drivers. ","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.795041Z","iopub.execute_input":"2022-10-05T19:03:20.795489Z","iopub.status.idle":"2022-10-05T19:03:20.801905Z","shell.execute_reply.started":"2022-10-05T19:03:20.795442Z","shell.execute_reply":"2022-10-05T19:03:20.801146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path_to_pan = '../input/iadf-school-ds-1/PAN.tif'","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.806086Z","iopub.execute_input":"2022-10-05T19:03:20.806385Z","iopub.status.idle":"2022-10-05T19:03:20.812690Z","shell.execute_reply.started":"2022-10-05T19:03:20.806354Z","shell.execute_reply":"2022-10-05T19:03:20.811875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pan = gdal.Open(path_to_pan)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.814891Z","iopub.execute_input":"2022-10-05T19:03:20.815206Z","iopub.status.idle":"2022-10-05T19:03:20.824086Z","shell.execute_reply.started":"2022-10-05T19:03:20.815176Z","shell.execute_reply":"2022-10-05T19:03:20.823135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"info = pan.GetMetadata_Dict()\nprint(info)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.825493Z","iopub.execute_input":"2022-10-05T19:03:20.825778Z","iopub.status.idle":"2022-10-05T19:03:20.834550Z","shell.execute_reply.started":"2022-10-05T19:03:20.825750Z","shell.execute_reply":"2022-10-05T19:03:20.833657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"proj = pan.GetSpatialRef()\nprint(proj)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.835622Z","iopub.execute_input":"2022-10-05T19:03:20.836149Z","iopub.status.idle":"2022-10-05T19:03:20.866193Z","shell.execute_reply.started":"2022-10-05T19:03:20.836113Z","shell.execute_reply":"2022-10-05T19:03:20.865046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gt = pan.GetGeoTransform()\n\nprint('Top Left pixel x-coordinate: {}, pixel width: {}, pixel rotation along x {}, Top Left pixel y-coordinate: {}, pixel rotation along y: {}, pixel height: {}'\n.format(gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]))\n\n# The pixel height is negative because the count starts from the top left pixel.","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.867588Z","iopub.execute_input":"2022-10-05T19:03:20.867899Z","iopub.status.idle":"2022-10-05T19:03:20.873713Z","shell.execute_reply.started":"2022-10-05T19:03:20.867868Z","shell.execute_reply":"2022-10-05T19:03:20.872881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cols = pan.RasterXSize\nrows = pan.RasterYSize\nbands = pan.RasterCount\n\nprint(rows, cols, bands)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.875198Z","iopub.execute_input":"2022-10-05T19:03:20.875491Z","iopub.status.idle":"2022-10-05T19:03:20.887860Z","shell.execute_reply.started":"2022-10-05T19:03:20.875460Z","shell.execute_reply":"2022-10-05T19:03:20.887016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pan_numpy = pan.ReadAsArray()\nprint(pan_numpy.shape)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.889376Z","iopub.execute_input":"2022-10-05T19:03:20.889664Z","iopub.status.idle":"2022-10-05T19:03:20.917341Z","shell.execute_reply.started":"2022-10-05T19:03:20.889636Z","shell.execute_reply":"2022-10-05T19:03:20.916369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(20, 10))\nplt.imshow(pan.ReadAsArray(), cmap='gray', clim=[0, 2048])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:20.918699Z","iopub.execute_input":"2022-10-05T19:03:20.919020Z","iopub.status.idle":"2022-10-05T19:03:21.477164Z","shell.execute_reply.started":"2022-10-05T19:03:20.918987Z","shell.execute_reply":"2022-10-05T19:03:21.476067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Starting to cut","metadata":{}},{"cell_type":"code","source":"path_to_pan = '../input/iadf-school-ds-1/PAN.tif'\npath_to_ms = '../input/iadf-school-ds-1/MS_LR.tif'","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:21.478614Z","iopub.execute_input":"2022-10-05T19:03:21.478907Z","iopub.status.idle":"2022-10-05T19:03:21.483134Z","shell.execute_reply.started":"2022-10-05T19:03:21.478877Z","shell.execute_reply":"2022-10-05T19:03:21.482116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pan = gdal.Open(path_to_pan)\nms = gdal.Open(path_to_ms)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:21.484434Z","iopub.execute_input":"2022-10-05T19:03:21.484743Z","iopub.status.idle":"2022-10-05T19:03:21.495471Z","shell.execute_reply.started":"2022-10-05T19:03:21.484714Z","shell.execute_reply":"2022-10-05T19:03:21.494696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p_cols = pan.RasterXSize\np_rows = pan.RasterYSize\np_bands = pan.RasterCount\n\nm_cols = ms.RasterXSize\nm_rows = ms.RasterYSize\nm_bands = ms.RasterCount\n\n\nprint(\"PAN shape:\", p_rows, p_cols, p_bands)\nprint(\"MS shape:\", m_rows, m_cols, m_bands)\n\nratio = p_rows / m_rows\nprint(\"Ratio:\", ratio)\n\nratio = int(ratio)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:21.496706Z","iopub.execute_input":"2022-10-05T19:03:21.497245Z","iopub.status.idle":"2022-10-05T19:03:21.510813Z","shell.execute_reply.started":"2022-10-05T19:03:21.497203Z","shell.execute_reply":"2022-10-05T19:03:21.510033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"RGB = (4, 2, 1) # Selecting the Red, Green and Blue bands for visualization.\n\nm = ms.ReadAsArray()\np = pan.ReadAsArray()\n\nm = np.moveaxis(m, 0, -1) # In general (and for visualization purpose in particular), the images are seen as Height x Width x Depth. \n                                          #GDAL returns an array of dimensions Depth x Height x Width. So a switch of the axes is needed.\nplt.figure(figsize=(20, 10))\nax1 = plt.subplot(1,2,1)\nplt.imshow(p, cmap='gray', clim=[0, 2048])\nax1.set_title('PAN')\nax2 = plt.subplot(1,2,2)\nplt.imshow(m[:, :, RGB] / 2048.0)\nax2.set_title('MS')\nplt.show()\n\n","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:21.512407Z","iopub.execute_input":"2022-10-05T19:03:21.512936Z","iopub.status.idle":"2022-10-05T19:03:22.441221Z","shell.execute_reply.started":"2022-10-05T19:03:21.512859Z","shell.execute_reply":"2022-10-05T19:03:22.440189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cutted_pan_v1 = pan.ReadAsArray()\ncutted_pan_v1 = cutted_pan_v1[256:1024 + 256, 512:1024 + 512]","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:22.442597Z","iopub.execute_input":"2022-10-05T19:03:22.442899Z","iopub.status.idle":"2022-10-05T19:03:22.451006Z","shell.execute_reply.started":"2022-10-05T19:03:22.442869Z","shell.execute_reply":"2022-10-05T19:03:22.450047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cutted_pan = pan.ReadAsArray(xoff=512, yoff=256, xsize=1024, ysize=1024) ## AXES ARE INVERTED\n#cutted_ms_wrong = ms.ReadAsArray(xoff=512, yoff=256, xsize=1024, ysize=1024) # WRONG!!\ncutted_ms = ms.ReadAsArray(xoff=512 // ratio, yoff=256// ratio, xsize=1024 // ratio, ysize=1024 // ratio, interleave='band')","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:22.452506Z","iopub.execute_input":"2022-10-05T19:03:22.452802Z","iopub.status.idle":"2022-10-05T19:03:22.459621Z","shell.execute_reply.started":"2022-10-05T19:03:22.452773Z","shell.execute_reply":"2022-10-05T19:03:22.458864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(20, 10))\nplt.subplot(1,2,1)\nplt.imshow(cutted_pan_v1, cmap='gray')\nplt.subplot(1,2,2)\nplt.imshow(cutted_pan, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:22.460723Z","iopub.execute_input":"2022-10-05T19:03:22.461028Z","iopub.status.idle":"2022-10-05T19:03:23.017129Z","shell.execute_reply.started":"2022-10-05T19:03:22.460998Z","shell.execute_reply":"2022-10-05T19:03:23.016045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cutted_ms = np.moveaxis(cutted_ms, 0, -1) \n\nplt.figure(figsize=(20, 10))\nax1 = plt.subplot(1,2,1)\nplt.imshow(cutted_pan, cmap='gray', clim=[0, 2048])\nax1.set_title('PAN')\nax2 = plt.subplot(1,2,2)\nplt.imshow(cutted_ms[:, :, RGB] / 2048.0)\nax2.set_title('MS')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:23.018572Z","iopub.execute_input":"2022-10-05T19:03:23.018882Z","iopub.status.idle":"2022-10-05T19:03:23.694635Z","shell.execute_reply.started":"2022-10-05T19:03:23.018850Z","shell.execute_reply":"2022-10-05T19:03:23.693760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"NIRYC = (6, 3, 1) # Selecting the NIR1, Yellow and Coastal bands for visualization.\nQ_MS = np.quantile(cutted_ms, (0.02, 0.98), (0, 1), keepdims=True)\nQ_PAN = np.quantile(cutted_pan, (0.02, 0.98), (0, 1), keepdims=True)\nplt.figure(figsize=(20, 10))\nax1 = plt.subplot(1,3,1)\n\nplt.imshow(((cutted_pan - Q_PAN[0, :, :]) / (Q_PAN[1, :, :] - Q_PAN[0, :, :]))[400:600, 100:400], cmap='gray')\nax1.set_title('PAN')\nax2 = plt.subplot(1,3,2)\nplt.imshow(np.clip((cutted_ms - Q_MS[0, :, :]) / (Q_MS[1, :, :] - Q_MS[0, :, :]), 0, 1)[400 // ratio: 600 // ratio, 100 // ratio:400 // ratio,  RGB])\nax2.set_title('MS (RGB)')\nax3 = plt.subplot(1,3,3)\nplt.imshow(np.clip((cutted_ms - Q_MS[0, :, :]) / (Q_MS[1, :, :] - Q_MS[0, :, :]), 0, 1)[400 // ratio: 600 // ratio, 100 // ratio:400 // ratio,  NIRYC])\nax3.set_title('MS (NIR-Y-C)')\nplt.show()\n\n# We can see a misalignment behaviour in the edge transition between the tennis fields.","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:23.695991Z","iopub.execute_input":"2022-10-05T19:03:23.696311Z","iopub.status.idle":"2022-10-05T19:03:24.169849Z","shell.execute_reply.started":"2022-10-05T19:03:23.696279Z","shell.execute_reply":"2022-10-05T19:03:24.168887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The problem of cutting with .ReadAsArray() function\n\nWe will lose all the informations seen above. ReadAsArray function returns only a matrix containing the values of the pixels and no other information about. We should try something different...\nLet's define a function to solve this issue","metadata":{}},{"cell_type":"code","source":"def patchify(img, outpath, starting_point, patch_dim):\n    \n    # Defining auxiliary variables\n\n    start_x = starting_point[0]\n    start_y = starting_point[1]\n\n    end_x = start_x + patch_dim[0]\n    end_y = start_y + patch_dim[1]\n\n    # Extracting information about the image\n    geotransf = img.GetGeoTransform()\n    proj = img.GetProjection()\n    rows = img.RasterYSize\n    cols = img.RasterXSize\n    bands = img.RasterCount\n    img_type = img.GetRasterBand(1).DataType\n\n    img_np = img.ReadAsArray()\n    if len(img_np.shape) == 2:\n        img_np = img_np[None, :, :]\n    \n\n\n    # Creating the patch\n\n    driver = gdal.GetDriverByName(\"GTiff\")\n    driver.Register()\n\n\n\n    cutted_img = img_np[:, start_x:end_x, start_y:end_y]\n\n    cutted_img_driver = driver.Create(outpath, \n                                      xsize=cutted_img.shape[2], \n                                      ysize=cutted_img.shape[1],\n                                      bands=bands,\n                                      eType=img_type\n                                      )\n    for b in range(bands):\n        cutted_img_driver.GetRasterBand(b + 1).WriteArray(cutted_img[b, :, :])\n        cutted_img_driver.GetRasterBand(b + 1).SetDescription(img.GetRasterBand(b+1).GetDescription())\n        cutted_img_driver.GetRasterBand(b + 1).ComputeStatistics(0)\n        cutted_img_driver.GetRasterBand(b + 1).FlushCache()\n\n    geo = np.copy(geotransf)\n    geo[0] += start_y * geo[1] + 0.5 * geo[1]\n    geo[3] += start_x * geo[5] + 0.5 * geo[5]\n    cutted_img_driver.SetGeoTransform(geo)\n    cutted_img_driver.SetProjection(proj)\n    \n    # Writing on disk\n    cutted_img_driver = None\n    \n    return","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:24.171357Z","iopub.execute_input":"2022-10-05T19:03:24.171672Z","iopub.status.idle":"2022-10-05T19:03:24.186068Z","shell.execute_reply.started":"2022-10-05T19:03:24.171641Z","shell.execute_reply":"2022-10-05T19:03:24.185019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Test it\n\noutroot_pan = './output/PairMax/PAN'\noutroot_ms = './output/PairMax/MS'\n\n# Creating the output folders\nif not os.path.exists(outroot_pan):\n    os.makedirs(outroot_pan)\n\nif not os.path.exists(outroot_ms):\n    os.makedirs(outroot_ms)\n\n# Defining the starting \nstarting_point = (256, 512) # BE CAREFUL: The axes are reversed!\npatch_size = (1024, 1024)\nsavepath_pan = os.path.join(outroot_pan, 'PairMax_Test.tif')\nsavepath_ms = os.path.join(outroot_ms, 'PairMax_Test.tif')\n\npatchify(pan, savepath_pan, starting_point, patch_size)\npatchify(ms, savepath_ms, (starting_point[0] // ratio, starting_point[1] // ratio), (patch_size[0] // ratio, patch_size[1] // ratio))\n","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:24.187465Z","iopub.execute_input":"2022-10-05T19:03:24.187778Z","iopub.status.idle":"2022-10-05T19:03:24.227330Z","shell.execute_reply.started":"2022-10-05T19:03:24.187749Z","shell.execute_reply":"2022-10-05T19:03:24.226270Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_pan = gdal.Open('./output/PairMax/PAN/PairMax_Test.tif').ReadAsArray()\ntest_ms = gdal.Open('./output/PairMax/MS/PairMax_Test.tif').ReadAsArray()\n\ntest_ms = np.moveaxis(test_ms, 0, -1)\n\nplt.figure(figsize=(20, 10))\nax1 = plt.subplot(1,2,1)\nplt.imshow(test_pan, cmap='gray', clim=[0,2048])\nax2 = plt.subplot(1,2,2)\nplt.imshow((test_ms / 2048.0)[:,:,RGB])","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:24.228852Z","iopub.execute_input":"2022-10-05T19:03:24.229318Z","iopub.status.idle":"2022-10-05T19:03:24.903363Z","shell.execute_reply.started":"2022-10-05T19:03:24.229269Z","shell.execute_reply":"2022-10-05T19:03:24.902368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Let's Patchify!","metadata":{}},{"cell_type":"code","source":"starting_point_x = 256 // ratio\nstarting_point_y = 256 // ratio\n\npatch_size = 512 // ratio\n\nn_horizontal_patches = 3\nn_vertical_patches = 2\n\nn = 1\nfor i in range(starting_point_x, starting_point_x + (patch_size * n_horizontal_patches), patch_size):\n    for j in range(starting_point_y, starting_point_y + (patch_size * n_vertical_patches), patch_size):\n        savepath_pan = os.path.join(outroot_pan, 'PairMax_{:03}.tif'.format(n))\n        savepath_ms = os.path.join(outroot_ms, 'PairMax_{:03}.tif'.format(n))\n\n        patchify(ms, savepath_ms, (j, i), (patch_size, patch_size))\n        patchify(pan, savepath_pan, (j * ratio, i * ratio), (patch_size * ratio, patch_size * ratio))\n        n += 1","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:24.904494Z","iopub.execute_input":"2022-10-05T19:03:24.904815Z","iopub.status.idle":"2022-10-05T19:03:24.977105Z","shell.execute_reply.started":"2022-10-05T19:03:24.904783Z","shell.execute_reply":"2022-10-05T19:03:24.976185Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Converting GeoTiff into Mat files\nBecause the most of algoritms about Remote Sensing are (still) written in MATLAB. MATLAB can manage correctly GeoTiff data, but a mat file may avoid errors during the analysis and the construction of the your custom pipeline. Furthermore, mat file may contain auxiliary own information and they can make the management of the dataset more orderly","metadata":{}},{"cell_type":"code","source":"outroot_mat = './output/PairMax/Mat'\n\n# Creating the output folders\nif not os.path.exists(outroot_mat):\n    os.makedirs(outroot_mat)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:24.978511Z","iopub.execute_input":"2022-10-05T19:03:24.978966Z","iopub.status.idle":"2022-10-05T19:03:24.984290Z","shell.execute_reply.started":"2022-10-05T19:03:24.978906Z","shell.execute_reply":"2022-10-05T19:03:24.983377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(1, 6 + 1):\n    load_pan_tiff = os.path.join(outroot_pan, 'PairMax_{:03}.tif'.format(i))\n    load_ms_tiff = os.path.join(outroot_ms, 'PairMax_{:03}.tif'.format(i))\n\n    pan = np.asarray(imageio.imread(load_pan_tiff))\n    ms = np.asarray(imageio.imread(load_ms_tiff))\n\n    if i == '6':\n        use_as = 'Validation'\n    elif i == '10':\n        use_as = 'Test'\n    else:\n        use_as = 'Training'\n\n    io.savemat(os.path.join(outroot_mat, 'PairMax_{:03}.mat'.format(i)), {'I_MS': ms, 'I_PAN': pan, 'sensor': 'WV3', 'usage': use_as, 'nbits': 11})\n\n    ","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:24.986059Z","iopub.execute_input":"2022-10-05T19:03:24.986464Z","iopub.status.idle":"2022-10-05T19:03:25.039736Z","shell.execute_reply.started":"2022-10-05T19:03:24.986422Z","shell.execute_reply":"2022-10-05T19:03:25.039001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Checking if everything works","metadata":{}},{"cell_type":"code","source":"load_mat = os.path.join(outroot_mat, 'PairMax_{:03}.mat'.format(5))\n\ndict = io.loadmat(load_mat)\n\nprint(dict.keys())","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:25.040900Z","iopub.execute_input":"2022-10-05T19:03:25.041222Z","iopub.status.idle":"2022-10-05T19:03:25.047708Z","shell.execute_reply.started":"2022-10-05T19:03:25.041192Z","shell.execute_reply":"2022-10-05T19:03:25.046983Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ms_from_mat = dict['I_MS']\npan_from_mat = dict['I_PAN']","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:25.048737Z","iopub.execute_input":"2022-10-05T19:03:25.049026Z","iopub.status.idle":"2022-10-05T19:03:25.056875Z","shell.execute_reply.started":"2022-10-05T19:03:25.048998Z","shell.execute_reply":"2022-10-05T19:03:25.055773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Q_MS = np.quantile(ms_from_mat, (0.02, 0.98), (0, 1), keepdims=True)\nplt.figure(figsize=(20, 10))\nax1 = plt.subplot(1,3,1)\nplt.imshow(pan_from_mat, cmap='gray', clim=[0, 2048])\nax1.set_title('PAN')\nax2 = plt.subplot(1,3,2)\nplt.imshow(np.clip((ms_from_mat - Q_MS[0, :, :]) / (Q_MS[1, :, :] - Q_MS[0, :, :]), 0, 1)[:, :, RGB])\nax2.set_title('MS (RGB)')\nax3 = plt.subplot(1,3,3)\nplt.imshow(np.clip((ms_from_mat - Q_MS[0, :, :]) / (Q_MS[1, :, :] - Q_MS[0, :, :]), 0, 1)[:, :, NIRYC])\nax3.set_title('MS (NIR-Y-C)')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:25.058436Z","iopub.execute_input":"2022-10-05T19:03:25.058756Z","iopub.status.idle":"2022-10-05T19:03:25.723833Z","shell.execute_reply.started":"2022-10-05T19:03:25.058714Z","shell.execute_reply":"2022-10-05T19:03:25.722798Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## It is not enough\n\nFor a supervised deep learning technique, we need a Ground-Truth. We resort to a downgrade of the available data. The original data will be used as a reference, while the downscaled data will be used to feed the network during the training.","metadata":{}},{"cell_type":"code","source":"def resize_images(img_ms, img_pan, ratio, sensor=None, mtf=None, apply_mtf_to_pan=False):\n\n    GNyq = []\n    GNyqPan = []\n    if (sensor is None) & (mtf is None):\n        MS_scale = (math.floor(img_ms.shape[0] / ratio), math.floor(img_ms.shape[1] / ratio), img_ms.shape[2])\n        PAN_scale = (math.floor(img_pan.shape[0] / ratio), math.floor(img_pan.shape[1] / ratio))\n        I_MS_LR = transform.resize(img_ms, MS_scale, order=3)\n        I_PAN_LR = transform.resize(img_pan, PAN_scale, order=3)\n\n        return I_MS_LR, I_PAN_LR\n\n    elif (sensor == 'QB') & (mtf is None):\n        GNyq = np.asarray([0.34, 0.32, 0.30, 0.22])  # Bands Order: B,G,R,NIR\n        GNyqPan = np.asarray([0.15])\n    elif ((sensor == 'Ikonos') or (sensor == 'IKONOS')) & (mtf is None):\n        GNyq = np.asarray([0.26, 0.28, 0.29, 0.28])  # Bands Order: B,G,R,NIR\n        GNyqPan = np.asarray([0.17])\n    elif (sensor == 'GeoEye1' or sensor == 'GE1') & (mtf is None):\n        GNyq = np.asarray([0.23, 0.23, 0.23, 0.23])  # Bands Order: B, G, R, NIR\n        GNyqPan = np.asarray([0.16])\n    elif (sensor == 'WV2') & (mtf is None):\n        GNyq = 0.35 * np.ones((1, 7))\n        GNyq = np.append(GNyq, 0.27)\n        GNyqPan = np.asarray([0.11])\n    elif (sensor == 'WV3') & (mtf is None):\n        GNyq = [0.325, 0.355, 0.360, 0.350, 0.365, 0.360, 0.335, 0.315]\n        GNyqPan = np.asarray([0.5])\n    elif mtf is not None:\n        GNyq = mtf['GNyq']\n        GNyqPan = np.asarray([mtf['GNyqPan']])\n\n    N = 41\n\n    r, c, b = img_ms.shape\n\n    img_ms = np.moveaxis(img_ms, -1, 0)\n    img_ms = np.expand_dims(img_ms, axis=0)\n\n    h = ut.nyquist_filter_generator(GNyq, ratio, N)\n    h = np.moveaxis(h, -1, 0)\n    h = np.expand_dims(h, axis=1)\n    h = h.astype('float32')\n\n    h = torch.from_numpy(h).type(torch.float32)\n\n    conv = nn.Conv2d(in_channels=b, out_channels=b, padding=math.ceil(N / 2),\n                     kernel_size=h.shape, groups=b, bias=False, padding_mode='replicate')\n\n    conv.weight.data = h\n    conv.weight.requires_grad = False\n\n    I_MS_LP = conv(torch.from_numpy(img_ms)).numpy()\n    I_MS_LP = np.squeeze(I_MS_LP)\n    I_MS_LP = np.moveaxis(I_MS_LP, 0, -1)\n    MS_scale = (math.floor(I_MS_LP.shape[0] / ratio), math.floor(I_MS_LP.shape[1] / ratio), I_MS_LP.shape[2])\n    PAN_scale = (math.floor(img_pan.shape[0] / ratio), math.floor(img_pan.shape[1] / ratio))\n\n    I_MS_LR = transform.resize(I_MS_LP, MS_scale, order=0)\n\n    if apply_mtf_to_pan:\n        img_pan = np.expand_dims(img_pan, [0, 1])\n\n        h = ut.nyquist_filter_generator(GNyqPan, ratio, N)\n        h = np.moveaxis(h, -1, 0)\n        h = np.expand_dims(h, axis=1)\n        h = h.astype('float32')\n\n        h = torch.from_numpy(h).type(torch.float32)\n\n        conv = nn.Conv2d(in_channels=1, out_channels=1, padding=math.ceil(N / 2),\n                         kernel_size=h.shape, groups=1, bias=False, padding_mode='replicate')\n\n        conv.weight.data = h\n        conv.weight.requires_grad = False\n\n        I_PAN_LP = conv(torch.from_numpy(img_pan)).numpy()\n        I_PAN_LP = np.squeeze(I_PAN_LP)\n        I_PAN_LR = transform.resize(I_PAN_LP, PAN_scale, order=0)\n\n    else:\n        I_PAN_LR = transform.resize(img_pan, PAN_scale, order=3)\n\n    return I_MS_LR, I_PAN_LR\n","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:25.725302Z","iopub.execute_input":"2022-10-05T19:03:25.725615Z","iopub.status.idle":"2022-10-05T19:03:25.754895Z","shell.execute_reply.started":"2022-10-05T19:03:25.725583Z","shell.execute_reply":"2022-10-05T19:03:25.753821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"outroot_downgraded = './output/PairMax/downgraded'\n\n# Creating the output folders\nif not os.path.exists(outroot_downgraded):\n    os.makedirs(outroot_downgraded)","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:25.759617Z","iopub.execute_input":"2022-10-05T19:03:25.760173Z","iopub.status.idle":"2022-10-05T19:03:25.770612Z","shell.execute_reply.started":"2022-10-05T19:03:25.760136Z","shell.execute_reply":"2022-10-05T19:03:25.769687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from ctypes import resize\n\n\nfor i in range(1, 6 + 1):\n    load_pan_tiff = os.path.join(outroot_pan, 'PairMax_{:03}.tif'.format(i))\n    load_ms_tiff = os.path.join(outroot_ms, 'PairMax_{:03}.tif'.format(i))\n\n    pan = np.asarray(imageio.imread(load_pan_tiff)).astype(np.float32)\n    ms = np.asarray(imageio.imread(load_ms_tiff)).astype(np.float32)\n\n    ms_dwngrd, pan_dwngrd = resize_images(ms, pan, ratio, 'WV3')\n\n    if i == '6':\n        use_as = 'Validation'\n    elif i == '10':\n        use_as = 'Test'\n    else:\n        use_as = 'Training'\n\n    io.savemat(os.path.join(outroot_downgraded, 'PairMax_{:03}.mat'.format(i)), {'I_MS': ms, 'I_PAN': pan, 'I_MS_LR': ms_dwngrd, 'I_PAN_LR': pan_dwngrd, 'sensor': 'WV3', 'usage': use_as, 'nbits': 11})","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:25.772442Z","iopub.execute_input":"2022-10-05T19:03:25.772865Z","iopub.status.idle":"2022-10-05T19:03:26.007129Z","shell.execute_reply.started":"2022-10-05T19:03:25.772824Z","shell.execute_reply":"2022-10-05T19:03:26.006345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Checking if everything works","metadata":{}},{"cell_type":"code","source":"load_mat = os.path.join(outroot_downgraded, 'PairMax_{:03}.mat'.format(5))\n\ndict = io.loadmat(load_mat)\n\nms_from_downgraded = dict['I_MS']\npan_from_downgraded = dict['I_PAN']\ndw_ms_from_downgraded = dict['I_MS_LR']\ndw_pan_from_downgraded = dict['I_PAN_LR']","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:26.008193Z","iopub.execute_input":"2022-10-05T19:03:26.008567Z","iopub.status.idle":"2022-10-05T19:03:26.015076Z","shell.execute_reply.started":"2022-10-05T19:03:26.008537Z","shell.execute_reply":"2022-10-05T19:03:26.014076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Q_MS = np.quantile(ms_from_downgraded, (0.02, 0.98), (0, 1), keepdims=True)\nQ_PAN = np.quantile(pan_from_downgraded, (0.02, 0.98), (0, 1), keepdims=True)\nplt.figure(figsize=(20, 10))\nax1 = plt.subplot(1,4,1)\nplt.imshow((np.clip((pan_from_downgraded - Q_PAN[0, :, :]) / (Q_PAN[1, :, :] - Q_PAN[0, :, :]), 0, 1)), cmap='gray')\nax1.set_title('PAN')\nax2 = plt.subplot(1,4,2)\nplt.imshow((np.clip((dw_pan_from_downgraded - Q_PAN[0, :, :]) / (Q_PAN[1, :, :] - Q_PAN[0, :, :]), 0, 1)), cmap='gray')\nax2.set_title('PAN Downgraded')\nax3 = plt.subplot(1,4,3)\nplt.imshow(np.clip((ms_from_downgraded - Q_MS[0, :, :]) / (Q_MS[1, :, :] - Q_MS[0, :, :]), 0, 1)[:, :, RGB])\nax3.set_title('MS (RGB)')\nax4 = plt.subplot(1,4,4)\nplt.imshow(np.clip((dw_ms_from_downgraded - Q_MS[0, :, :]) / (Q_MS[1, :, :] - Q_MS[0, :, :]), 0, 1)[:, :, RGB])\nax4.set_title('MS (RGB) Downgraded')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-05T19:03:26.016413Z","iopub.execute_input":"2022-10-05T19:03:26.016799Z","iopub.status.idle":"2022-10-05T19:03:26.782157Z","shell.execute_reply.started":"2022-10-05T19:03:26.016768Z","shell.execute_reply":"2022-10-05T19:03:26.781010Z"},"trusted":true},"execution_count":null,"outputs":[]}]}