06_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b9aaeef7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona\n",
    "import os\n",
    "from ops.ops import load_json, create_exps_paths, load_exp\n",
    "from tqdm.notebook import tqdm_notebook\n",
    "from osgeo import gdal\n",
    "from osgeo import gdalconst\n",
    "import numpy as np\n",
    "from tensorflow.keras.utils import to_categorical\n",
    "from skimage.morphology import area_opening\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "fdb26936",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_n = 5\n",
    "exp = load_exp(exp_n)\n",
    "paths = load_json(os.path.join('conf', 'paths.json'))\n",
    "shp = load_json(os.path.join('conf', 'shp.json'))\n",
    "shp_path = paths['shp']\n",
    "conf = load_json(os.path.join('conf', 'conf.json'))\n",
    "test_tiles_path = os.path.join('img', 'patches')\n",
    "exps_path, exp_path, models_path, results_path, predictions_path, visual_path, logs_path = create_exps_paths(exp_n)\n",
    "\n",
    "img_source = conf['img_source']\n",
    "grid_shp = shp[f'shp_download_{img_source}']\n",
    "grid_save = os.path.join('shp', f\"{grid_shp}.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b732bc7d",
   "metadata": {},
   "source": [
    "## Create Geotiff Predictions Probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "6aad4490",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eee72d1d482b49c88dc5e4450da3d4e1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Tiles evaluated:   0%|          | 0/30 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with fiona.open(grid_save) as grid:\n",
    "    for feat in tqdm_notebook(grid, desc = 'Tiles evaluated'):\n",
    "        if feat['properties']['dataset'] !=2:\n",
    "            continue\n",
    "        feat_id = int(feat['properties']['id'])\n",
    "        \n",
    "        pred = np.load(os.path.join(predictions_path, f'{feat_id}_fus.npy'))\n",
    "        fz = np.load(os.path.join(test_tiles_path, f'{feat_id}.npz'))\n",
    "        label = fz['label'].reshape(fz['shape'])\n",
    "        pred_p = pred[:,:,1:2]\n",
    "        pred_p[label==2] = 0\n",
    "        \n",
    "        in_data = gdal.Open( os.path.join('img', 'labels', f'label_{feat_id}.tif'), gdalconst.GA_ReadOnly)\n",
    "        geo_transform = in_data.GetGeoTransform()\n",
    "        x_min = geo_transform[0]\n",
    "        y_max = geo_transform[3]\n",
    "        x_max = x_min + geo_transform[1] * in_data.RasterXSize\n",
    "        y_min = y_max + geo_transform[5] * in_data.RasterYSize\n",
    "        x_res = in_data.RasterXSize\n",
    "        y_res = in_data.RasterYSize\n",
    "        crs = in_data.GetSpatialRef()\n",
    "        proj = in_data.GetProjection()\n",
    "        \n",
    "        output = os.path.join(visual_path, f'{feat_id}_prediction_prob_exp_{exp_n}.tif')\n",
    "        target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Float32)\n",
    "        \n",
    "        target_ds.SetGeoTransform(geo_transform)\n",
    "        target_ds.SetSpatialRef(crs)\n",
    "        target_ds.SetProjection(proj)\n",
    "        \n",
    "        band = target_ds.GetRasterBand(1)\n",
    "        band.WriteArray(pred_p[:,:,0], 0, 0)\n",
    "        \n",
    "        band.FlushCache()\n",
    "        target_ds = None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dbd1e85",
   "metadata": {},
   "source": [
    "## Create Geotiff Cloud Predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "264ff178",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e7db4c2a98b64d0ca7bd8411a59be842",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Tiles evaluated:   0%|          | 0/30 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'exps\\\\exp_3\\\\results\\\\predictions\\\\975_cloud.npy'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32mC:\\Users\\FELFER~1\\AppData\\Local\\Temp/ipykernel_14876/2334147370.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      5\u001b[0m         \u001b[0mfeat_id\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfeat\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'properties'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'id'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 7\u001b[1;33m         \u001b[0mpred\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mload\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mos\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpredictions_path\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34mf'{feat_id}_cloud.npy'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      8\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      9\u001b[0m         \u001b[0min_data\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgdal\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mOpen\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mos\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'img'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'labels'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34mf'label_{feat_id}.tif'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgdalconst\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mGA_ReadOnly\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\.conda\\envs\\ferrari_1\\lib\\site-packages\\numpy\\lib\\npyio.py\u001b[0m in \u001b[0;36mload\u001b[1;34m(file, mmap_mode, allow_pickle, fix_imports, encoding)\u001b[0m\n\u001b[0;32m    415\u001b[0m             \u001b[0mown_fid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mFalse\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    416\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 417\u001b[1;33m             \u001b[0mfid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mstack\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0menter_context\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mopen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mos_fspath\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfile\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"rb\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    418\u001b[0m             \u001b[0mown_fid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    419\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'exps\\\\exp_3\\\\results\\\\predictions\\\\975_cloud.npy'"
     ]
    }
   ],
   "source": [
    "with fiona.open(grid_save) as grid:\n",
    "    for feat in tqdm_notebook(grid, desc = 'Tiles evaluated'):\n",
    "        if feat['properties']['dataset'] !=2:\n",
    "            continue\n",
    "        feat_id = int(feat['properties']['id'])\n",
    "        \n",
    "        pred = np.load(os.path.join(predictions_path, f'{feat_id}_cloud.npy'))\n",
    "        \n",
    "        in_data = gdal.Open( os.path.join('img', 'labels', f'label_{feat_id}.tif'), gdalconst.GA_ReadOnly)\n",
    "        geo_transform = in_data.GetGeoTransform()\n",
    "        x_min = geo_transform[0]\n",
    "        y_max = geo_transform[3]\n",
    "        x_max = x_min + geo_transform[1] * in_data.RasterXSize\n",
    "        y_min = y_max + geo_transform[5] * in_data.RasterYSize\n",
    "        x_res = in_data.RasterXSize\n",
    "        y_res = in_data.RasterYSize\n",
    "        crs = in_data.GetSpatialRef()\n",
    "        proj = in_data.GetProjection()\n",
    "        \n",
    "        output = os.path.join(visual_path, f'{feat_id}_cloud_prediction_exp_{exp_n}.tif')\n",
    "        target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Float32 )\n",
    "        \n",
    "        target_ds.SetGeoTransform(geo_transform)\n",
    "        target_ds.SetSpatialRef(crs)\n",
    "        target_ds.SetProjection(proj)\n",
    "        \n",
    "        band = target_ds.GetRasterBand(1)\n",
    "        band.WriteArray(np.clip(pred[:,:,0], 0, 1), 0, 0)\n",
    "        \n",
    "        band.FlushCache()\n",
    "        target_ds = None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c79c088c",
   "metadata": {},
   "source": [
    "## Create Geotiff Binary Predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0f18bfd3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f1527cf9863445b0a3ce68d85db08167",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Tiles evaluated:   0%|          | 0/30 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with fiona.open(grid_save) as grid:\n",
    "    for feat in tqdm_notebook(grid, desc = 'Tiles evaluated'):\n",
    "        if feat['properties']['dataset'] !=2:\n",
    "            continue\n",
    "        feat_id = int(feat['properties']['id'])\n",
    "        \n",
    "        pred = np.load(os.path.join(predictions_path, f'{feat_id}.npy'))\n",
    "        fz = np.load(os.path.join(test_tiles_path, f'{feat_id}.npz'))\n",
    "        label = fz['label']\n",
    "        pred_b = np.zeros_like(label)\n",
    "        pred_b[pred[:,:,1]>0.5] = 1\n",
    "        pred_b[label == 2] = 0\n",
    "        pred_b = area_opening(pred_b.squeeze(), 625)\n",
    "        \n",
    "        \n",
    "        \n",
    "        in_data = gdal.Open( os.path.join('img', 'labels', f'label_{feat_id}.tif'), gdalconst.GA_ReadOnly)\n",
    "        geo_transform = in_data.GetGeoTransform()\n",
    "        x_min = geo_transform[0]\n",
    "        y_max = geo_transform[3]\n",
    "        x_max = x_min + geo_transform[1] * in_data.RasterXSize\n",
    "        y_min = y_max + geo_transform[5] * in_data.RasterYSize\n",
    "        x_res = in_data.RasterXSize\n",
    "        y_res = in_data.RasterYSize\n",
    "        crs = in_data.GetSpatialRef()\n",
    "        proj = in_data.GetProjection()\n",
    "        \n",
    "        output = os.path.join(visual_path, f'{feat_id}_binary_pred_binary_exp_{exp_n}.tif')\n",
    "        target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)\n",
    "        \n",
    "        target_ds.SetGeoTransform(geo_transform)\n",
    "        target_ds.SetSpatialRef(crs)\n",
    "        target_ds.SetProjection(proj)\n",
    "        \n",
    "        band = target_ds.GetRasterBand(1)\n",
    "        band.WriteArray(pred_b, 0, 0)\n",
    "        \n",
    "        band.FlushCache()\n",
    "        target_ds = None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "16f369ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2896, 2895, 1)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred_b.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d54e6504",
   "metadata": {},
   "source": [
    "## Create Geotiff error binary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c63fefc4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b9f9e70ae73c44c582a6579cf5ccdcfc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Tiles evaluated:   0%|          | 0/30 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with fiona.open(grid_save) as grid:\n",
    "    for feat in tqdm_notebook(grid, desc = 'Tiles evaluated'):\n",
    "        if feat['properties']['dataset'] !=2:\n",
    "            continue\n",
    "        feat_id = int(feat['properties']['id'])\n",
    "        \n",
    "        pred = np.load(os.path.join(predictions_path, f'{feat_id}.npy'))\n",
    "        fz = np.load(os.path.join(test_tiles_path, f'{feat_id}.npz'))\n",
    "        label = fz['label']\n",
    "        pred_b = np.zeros_like(label)\n",
    "        pred_b[pred[:,:,1]>0.5] = 1\n",
    "        pred_b[label == 2] = 0\n",
    "        pred_b = area_opening(pred_b.squeeze(), 625)\n",
    "        pred_b = np.expand_dims(pred_b, axis=-1)\n",
    "        \n",
    "        error_b = np.zeros_like(label)\n",
    "        error_b[label==2] = 1 #background\n",
    "        error_b[np.logical_and(label==1, pred_b==0)] = 2 #FN\n",
    "        error_b[np.logical_and(label==0, pred_b==1)] = 3 #FP\n",
    "        error_b[np.logical_and(label==1, pred_b==1)] = 4 #TP\n",
    "        error_b[np.logical_and(label==0, pred_b==0)] = 5 #TN\n",
    "        \n",
    "        \n",
    "        \n",
    "        in_data = gdal.Open( os.path.join('img', 'labels', f'label_{feat_id}.tif'), gdalconst.GA_ReadOnly)\n",
    "        geo_transform = in_data.GetGeoTransform()\n",
    "        x_min = geo_transform[0]\n",
    "        y_max = geo_transform[3]\n",
    "        x_max = x_min + geo_transform[1] * in_data.RasterXSize\n",
    "        y_min = y_max + geo_transform[5] * in_data.RasterYSize\n",
    "        x_res = in_data.RasterXSize\n",
    "        y_res = in_data.RasterYSize\n",
    "        crs = in_data.GetSpatialRef()\n",
    "        proj = in_data.GetProjection()\n",
    "        \n",
    "        output = os.path.join(visual_path, f'{feat_id}_binary_error_b_exp_{exp_n}.tif')\n",
    "        target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)\n",
    "        \n",
    "        target_ds.SetGeoTransform(geo_transform)\n",
    "        target_ds.SetSpatialRef(crs)\n",
    "        target_ds.SetProjection(proj)\n",
    "        \n",
    "        band = target_ds.GetRasterBand(1)\n",
    "        band.WriteArray(error_b[:,:,0], 0, 0)\n",
    "        \n",
    "        band.FlushCache()\n",
    "        target_ds = None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca37a9a4",
   "metadata": {},
   "source": [
    "## Generate minimum error map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7382f4e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "be3ff3de698e4b45b0ae591a95d8dd10",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Tiles evaluated:   0%|          | 0/30 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "opt_exp = 1\n",
    "sar_exp = 2\n",
    "fus_exp = 3\n",
    "\n",
    "\n",
    "output_path = os.path.join('exps', 'compile')\n",
    "\n",
    "with fiona.open(grid_save) as grid:\n",
    "    for feat in tqdm_notebook(grid, desc = 'Tiles evaluated'):\n",
    "        if feat['properties']['dataset'] !=2:\n",
    "            continue\n",
    "        feat_id = int(feat['properties']['id'])\n",
    "\n",
    "        opt_pred = np.load(os.path.join('exps', f'exp_{opt_exp}', 'results', 'predictions', f'{feat_id}.npy'))[:,:,1:2]\n",
    "        sar_pred = np.load(os.path.join('exps', f'exp_{sar_exp}', 'results', 'predictions', f'{feat_id}.npy'))[:,:,1:2]\n",
    "        fus_pred = np.load(os.path.join('exps', f'exp_{fus_exp}', 'results', 'predictions', f'{feat_id}.npy'))[:,:,1:2]\n",
    "        fz = np.load(os.path.join(test_tiles_path, f'{feat_id}.npz'))\n",
    "        label = fz['label']\n",
    "\n",
    "        error_min = np.zeros_like(label)\n",
    "        gt = label ==1\n",
    "        opt_error = np.abs(opt_pred - gt)\n",
    "        sar_error = np.abs(sar_pred - gt)\n",
    "        fus_error = np.abs(fus_pred - gt)\n",
    "\n",
    "        errors = np.concatenate([opt_error, sar_error, fus_error], axis=2)\n",
    "        errors_arg = np.expand_dims(np.argmin(errors, axis=2), axis=-1)+1\n",
    "        errors_arg[label==2] = 0\n",
    "\n",
    "\n",
    "        in_data = gdal.Open( os.path.join('img', 'labels', f'label_{feat_id}.tif'), gdalconst.GA_ReadOnly)\n",
    "        geo_transform = in_data.GetGeoTransform()\n",
    "        x_min = geo_transform[0]\n",
    "        y_max = geo_transform[3]\n",
    "        x_max = x_min + geo_transform[1] * in_data.RasterXSize\n",
    "        y_min = y_max + geo_transform[5] * in_data.RasterYSize\n",
    "        x_res = in_data.RasterXSize\n",
    "        y_res = in_data.RasterYSize\n",
    "        crs = in_data.GetSpatialRef()\n",
    "        proj = in_data.GetProjection()\n",
    "        \n",
    "        output = os.path.join(output_path, f'{feat_id}_error.tif')\n",
    "        target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)\n",
    "        \n",
    "        target_ds.SetGeoTransform(geo_transform)\n",
    "        target_ds.SetSpatialRef(crs)\n",
    "        target_ds.SetProjection(proj)\n",
    "        \n",
    "        band = target_ds.GetRasterBand(1)\n",
    "        band.WriteArray(np.squeeze(errors_arg), 0, 0)\n",
    "        \n",
    "        band.FlushCache()\n",
    "        target_ds = None\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e3ac214",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f4b8082",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4af1d599",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "057bcf1b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "6984327cd4a671f1c687c03fe0b5d8b7594868635f0d5dbeb7c4ac1a044793af"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}