1-preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare the synthetic test cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from geone.img import readImageGslib, writeImageGslib, sampleFromImage, pointSetToImage\n",
    "from geone.img import writePointSetGslib, readPointSetGslib\n",
    "from geone.imgplot import drawImage2D\n",
    "from geone.deesseinterface import DeesseInput, deesseRun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_DIR = 'data/'\n",
    "SAMPLES_DIR = 'samples/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "COLOR_SCHEME_BINARY = [ \n",
    "        [x/255 for x in [166,206,227]], \n",
    "        [x/255 for x in [31,120,180]],\n",
    "      ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmark: training image selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training images\n",
    "\n",
    "The training images originally apeared in the publication: Pérez, C., Mariethoz, G., Ortiz, J. M., 2014. Verifying the high-order consistency of training images with data for multiple-point geostatistics. Computers & Geosciences 70, 190 – 205. https://doi.org/10.1016/j.cageo.2014.06.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_and_plot_image(filename):\n",
    "    image = readImageGslib(filename)\n",
    "    drawImage2D(image, categ=True, categCol=COLOR_SCHEME_BINARY)\n",
    "    return image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_A = get_and_plot_image(DATA_DIR+'A.gslib')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_B = get_and_plot_image(DATA_DIR+'B.gslib')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_C = get_and_plot_image(DATA_DIR+'C.gslib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Synthetic realities "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_and_plot_DeeSse_simulation(ti):\n",
    "    # build DeeSse input\n",
    "    deesse_input = DeesseInput(\n",
    "        nx=100, ny=100, nz=1,     # dimension of the simulation grid (number of cells)\n",
    "        sx=1.0, sy=1.0, sz=1.0,   # cells units in the simulation grid (here are the default values)\n",
    "        ox=0.0, oy=0.0, oz=0.0,   # origin of the simulation grid (here are the default values)\n",
    "        nv=1, varname='code',     # number of variable(s), name of the variable(s)\n",
    "        nTI=1, TI=ti,             # number of TI(s), TI (class dsi.Img)\n",
    "        distanceType=0,           # distance type: proportion of mismatching nodes (categorical var., default)\n",
    "        nneighboringNode=60,      # max. number of neighbors (for the patterns)\n",
    "        distanceThreshold=0.05,   # acceptation threshold (for distance between patterns)\n",
    "        maxScanFraction=0.25,     # max. scanned fraction of the TI (for simulation of each cell)\n",
    "        npostProcessingPathMax=1, # number of post-processing path(s)\n",
    "        seed=201912,                 # seed (initialization of the random number generator)\n",
    "        nrealization=1)       # number of realization(s)\n",
    "    \n",
    "    # run the simulation\n",
    "    image = deesseRun(deesse_input)['sim'][0]\n",
    "    \n",
    "    #plot the result\n",
    "    drawImage2D(image, categ=True, categCol=COLOR_SCHEME_BINARY)\n",
    "    return image\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reality_A = get_and_plot_DeeSse_simulation(ti=ti_A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reality_B = get_and_plot_DeeSse_simulation(ti=ti_B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reality_C = get_and_plot_DeeSse_simulation(ti=ti_C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample with different sampling rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SAMPLE_SIZES = [25, 50, 75, 100, 125, 150, 175, 200, 400, 1600]\n",
    "SEED = 201912"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_and_write_gslib(image, size, seed, name):\n",
    "        sampled_point_set = sampleFromImage(image, size=size, seed=SEED+i)\n",
    "        writePointSetGslib(sampled_point_set,\n",
    "                                     SAMPLES_DIR+'sample_{name}_{size}.gslib'.format(name=name, size=size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for reality, name in zip([reality_A, reality_B, reality_C], ['A', 'B', 'C']):\n",
    "    for i, size, in enumerate(SAMPLE_SIZES):\n",
    "        sample_and_write_gslib(reality, size, SEED+i, name)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NX, NY, NZ = 100, 100, 1\n",
    "\n",
    "def get_and_plot_point_set(filename):\n",
    "    point_set = readPointSetGslib(filename)\n",
    "    image = pointSetToImage(point_set, NX, NY, NZ)\n",
    "    drawImage2D(image, categ=True, categCol=COLOR_SCHEME_BINARY)\n",
    "    return image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observations_A = get_and_plot_point_set(SAMPLES_DIR+'sample_A_1600.gslib')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observations_B = get_and_plot_point_set(SAMPLES_DIR+'sample_B_1600.gslib')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observations_C = get_and_plot_point_set(SAMPLES_DIR+'sample_C_1600.gslib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Roussillon: synthetic test case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "COLOR_SCHEME_ROUSSILLON = [ \n",
    "        [x/255 for x in [166,206,227]],\n",
    "        [x/255 for x in [178,223,138]],   \n",
    "        [x/255 for x in [31,120,180]],\n",
    "        [x/255 for x in [51,160,44]],\n",
    "      ]\n",
    "LEGEND = ['alluvial fan', 'flood plain', 'splay', 'river bed']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_true = readImageGslib(DATA_DIR+'trueTI.gslib')\n",
    "drawImage2D(ti_true, categ=True, categCol=COLOR_SCHEME_ROUSSILLON, cticklabels=LEGEND)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drawImage2D(ti_true, iv=1, cmap='cividis', vmin=0.0, vmax=1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = readImageGslib(DATA_DIR+'mask.gslib')\n",
    "drawImage2D(mask, categ=True, categCol=COLOR_SCHEME_BINARY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trend = readImageGslib(DATA_DIR+'trend.gslib')\n",
    "drawImage2D(trend, cmap='cividis', vmin=0.0, vmax=1.0, excludedVal=-9999999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Orientation\n",
    "im_angle = readImageGslib(DATA_DIR+'orientation.gslib')\n",
    "drawImage2D(im_angle, cmap='twilight_shifted', vmin=-180, vmax=180)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Synthetic reality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx, ny, nz = mask.nx, mask.ny, mask.nz      # number of cells\n",
    "sx, sy, sz = mask.sx, mask.sy, mask.sz      # cell unit\n",
    "ox, oy, oz = mask.ox, mask.oy, mask.oz      # origin (corner of the \"first\" grid cell)\n",
    "\n",
    "deesse_input_roussilon = DeesseInput(\n",
    "    nx=nx, ny=ny, nz=nz,\n",
    "    sx=sx, sy=sy, sz=sz,\n",
    "    ox=ox, oy=oy, oz=oz,\n",
    "    nv=2, varname=['Facies', 'trend'],\n",
    "    nTI=1, TI=ti_true,\n",
    "    mask=mask.val,\n",
    "    rotationUsage=1,            # use rotation without tolerance\n",
    "    rotationAzimuthLocal=True,  #    rotation according to azimuth: local\n",
    "    rotationAzimuth=im_angle.val[0,:,:,:],      #    rotation azimuth: map of values\n",
    "    dataImage=trend,\n",
    "    outputVarFlag=[True, False],\n",
    "    distanceType=[0,1],\n",
    "    nneighboringNode=[50,1],\n",
    "    distanceThreshold=[0.05, 0.05],\n",
    "    maxScanFraction=0.5,\n",
    "    npostProcessingPathMax=1,\n",
    "    seed=201912,\n",
    "    nrealization=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roussillon_simulation = deesseRun(deesse_input=deesse_input_roussilon, nthreads=1)['sim'][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample synthetic observation set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point_set_roussillon_50 = sampleFromImage(roussillon_simulation, size=50, seed=20191201, mask=mask)\n",
    "point_set_roussillon_150 = sampleFromImage(roussillon_simulation, size=150, seed=20191201, mask=mask)\n",
    "point_set_roussillon_600 = sampleFromImage(roussillon_simulation, size=600, seed=20191201, mask=mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drawImage2D(roussillon_simulation,  categ=True, categColbad='white', categCol=COLOR_SCHEME_ROUSSILLON, cticklabels=LEGEND)\n",
    "plt.scatter(point_set_roussillon_50.x(), point_set_roussillon_50.y(), marker= 'x', s=20, c='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drawImage2D(roussillon_simulation,  categ=True, categColbad='white', categCol=COLOR_SCHEME_ROUSSILLON, cticklabels=LEGEND)\n",
    "plt.scatter(point_set_roussillon_150.x(), point_set_roussillon_150.y(), marker= 'x', s=20, c='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "drawImage2D(roussillon_simulation,  categ=True, categColbad='white', categCol=COLOR_SCHEME_ROUSSILLON, cticklabels=LEGEND)\n",
    "plt.scatter(point_set_roussillon_600.x(), point_set_roussillon_600.y(), marker= 'x', s=10, c='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writePointSetGslib(point_set_roussillon_50, SAMPLES_DIR+'roussillon_observations_50.gslib')\n",
    "writePointSetGslib(point_set_roussillon_150, SAMPLES_DIR+'roussillon_observations_150.gslib')\n",
    "writePointSetGslib(point_set_roussillon_600, SAMPLES_DIR+'roussillon_observations_600.gslib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figures for publication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Benchmark test case\n",
    "FONT_SIZE = 16\n",
    "FIG_DIR = 'figures/'\n",
    "DPI = 600\n",
    "\n",
    "import matplotlib\n",
    "\n",
    "matplotlib.rcParams[\"image.interpolation\"] = None\n",
    "matplotlib.rcParams['pdf.fonttype'] = 42\n",
    "matplotlib.rcParams['ps.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training image selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ti_selection_plot(image, filename, title, colorbar=True):\n",
    "    plt.figure(figsize=(5,6))\n",
    "    drawImage2D(image,   title=title,\n",
    "                removeColorbar=not colorbar,\n",
    "                              categ=True,\n",
    "                              categColbad='white',\n",
    "                              categCol=COLOR_SCHEME_BINARY,\n",
    "                              title_fontsize=FONT_SIZE+2,\n",
    "                              cticklabels_fontsize=FONT_SIZE,\n",
    "                              xlabels_fontsize=FONT_SIZE,\n",
    "                              ylabels_fontsize=FONT_SIZE,\n",
    "                              xticklabels = [0,100],\n",
    "                              yticklabels = [0, 100],\n",
    "                              xticklabels_fontsize=FONT_SIZE,\n",
    "                              yticklabels_fontsize=FONT_SIZE,\n",
    "                              xticks=[0,100],\n",
    "                              yticks=[0,100],\n",
    "                              ylabel_rotation=0,\n",
    "                \n",
    "                             )\n",
    "    plt.savefig(filename, dpi=DPI)\n",
    "    #!convert -trim $filename $filename\n",
    "    !pdfcrop $filename $filename"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Training images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(ti_A, FIG_DIR+'ti_a.pdf', 'a)', colorbar=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(ti_B, FIG_DIR+'ti_b.pdf', 'b)', colorbar=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(ti_C, FIG_DIR+'ti_c.pdf', 'c)', colorbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Synthetic realities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(reality_A, FIG_DIR+'reality_a.pdf', 'a)', colorbar=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(reality_B, FIG_DIR+'reality_b.pdf', 'b)', colorbar=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(reality_C, FIG_DIR+'reality_c.pdf', 'c)', colorbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Observation sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(observations_A, FIG_DIR+'observations_a.pdf', 'a)', colorbar=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(observations_B, FIG_DIR+'observations_b.pdf', 'b)', colorbar=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_selection_plot(observations_C, FIG_DIR+'observations_c.pdf', 'c)', colorbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Roussillon parameter selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EXCLUDED_VAL = -9999999"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Area with the trend and orientation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roussillon_plot(image, filename, title, cmap, vmin, vmax, cticks, cticklabels, colorbar=True, categ=False, iv=0):\n",
    "    fig = plt.figure(figsize=(5,5))\n",
    "    fig.subplots_adjust(left=0.05, right=0.9)\n",
    "    xmin, xmax = [int(x) for x in [image.xmin(), image.xmax()]]\n",
    "    ymin, ymax = [int(y) for y in [image.ymin(), image.ymax()]]\n",
    "    drawImage2D(image, iv=iv, excludedVal=EXCLUDED_VAL, cmap=cmap, vmin=vmin, vmax=vmax, title=title,\n",
    "                removeColorbar=not colorbar,\n",
    "                categ=categ,\n",
    "                categCol=COLOR_SCHEME_ROUSSILLON,\n",
    "                title_fontsize=FONT_SIZE+2,\n",
    "                cticklabels_fontsize=FONT_SIZE,\n",
    "                xlabels_fontsize=FONT_SIZE,\n",
    "                ylabels_fontsize=FONT_SIZE,\n",
    "                xticklabels = [xmin, xmax],\n",
    "                yticklabels = [ymin, ymax],\n",
    "                xticklabels_fontsize=FONT_SIZE,\n",
    "                yticklabels_fontsize=FONT_SIZE,\n",
    "                xticks=[xmin, xmax],\n",
    "                yticks=[ymin, ymax],\n",
    "                cticks=cticks,\n",
    "                cticklabels=cticklabels,\n",
    "                ylabel_rotation=0,\n",
    "               )\n",
    "    plt.savefig(filename, dpi=DPI)\n",
    "    #!convert -trim $filename $filename\n",
    "    !pdfcrop $filename $filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roussillon_plot(trend, filename=FIG_DIR+'trend.pdf', title='a)', cmap='cividis', vmin=0.0, vmax=1.0, cticks=[0.0, 1.0], cticklabels=[0.0, 1.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In order to limit the orientation over the domain, we must tweak the array values\n",
    "# We create copy of the orientation and put values outside mask to exclude_value\n",
    "im_angle_plot = copy.deepcopy(im_angle)\n",
    "\n",
    "# Iterate over the orientation array and the mask simultaneously and modify in-place\n",
    "with np.nditer(im_angle_plot.val[0:1,:], op_flags=['readwrite']) as iterator:\n",
    "    for angle_el, mask_el in zip(iterator, np.nditer(mask.val.flatten())):\n",
    "        if mask_el == 0:\n",
    "            angle_el[...] = EXCLUDED_VAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roussillon_plot(im_angle_plot, filename=FIG_DIR+'orientation.pdf', title='b)', cmap='twilight_shifted', vmin=-180, vmax=180, cticks=[-180, 0, 180], cticklabels=['$-180^\\circ$','$0^\\circ$', '$180^\\circ$'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Training images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roussillon_cat_plot(image, filename, title, colorbar=True):\n",
    "    fig = plt.figure(figsize=(5,5))\n",
    "    fig.subplots_adjust(left=0.05, right=0.9)\n",
    "    xmin, xmax = [int(x) for x in [image.xmin(), image.xmax()]]\n",
    "    ymin, ymax = [int(y) for y in [image.ymin(), image.ymax()]]\n",
    "    drawImage2D(image, excludedVal=EXCLUDED_VAL,\n",
    "                title=title,\n",
    "                removeColorbar=not colorbar,\n",
    "                categ=True,\n",
    "                categCol=COLOR_SCHEME_ROUSSILLON,\n",
    "                cticklabels=LEGEND,\n",
    "                title_fontsize=FONT_SIZE+2,\n",
    "                cticklabels_fontsize=FONT_SIZE,\n",
    "                xlabels_fontsize=FONT_SIZE,\n",
    "                ylabels_fontsize=FONT_SIZE,\n",
    "                xticklabels = [xmin, xmax],\n",
    "                yticklabels = [ymin, ymax],\n",
    "                xticklabels_fontsize=FONT_SIZE,\n",
    "                yticklabels_fontsize=FONT_SIZE,\n",
    "                xticks=[xmin, xmax],\n",
    "                yticks=[ymin, ymax],\n",
    "                ylabel_rotation=0,\n",
    "               )\n",
    "    plt.savefig(filename, dpi=DPI, bbox_inches=\"tight\")\n",
    "    !pdfcrop $filename $filename\n",
    "    #!convert -trim $filename $filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roussillon_cat_plot(ti_true, filename=FIG_DIR+'ti_true.pdf', title='a)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roussillon_plot(ti_true, filename=FIG_DIR+'ti_true_trend.pdf', title='b)', cmap='cividis', vmin=0.0, vmax=1.0, cticks=[0.0, 1.0], cticklabels=[0.0, 1.0], iv=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti_analog = readImageGslib(DATA_DIR+'analogTI.gslib')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roussillon_cat_plot(ti_analog, filename=FIG_DIR+'ti_analog.pdf', title='a)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roussillon_plot(ti_analog, filename=FIG_DIR+'ti_analog_trend.pdf', title='b)', cmap='cividis', vmin=0.0, vmax=1.0, cticks=[0.0, 1.0], cticklabels=[0.0, 1.0], iv=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = roussillon_simulation\n",
    "\n",
    "for nsamples, title in zip(['50', '150', '600'], ['a)', 'b)', 'c)']):\n",
    "    filename = FIG_DIR + \"roussillon_observations_{}.pdf\".format(nsamples)\n",
    "    point_set_roussillon = readPointSetGslib(SAMPLES_DIR+\"roussillon_observations_{}.gslib\".format(nsamples))\n",
    "    #image = pointSetToImage(point_set_roussillon, image.nx, image.ny, image.nz)\n",
    "    fig = plt.figure(figsize=(5,5))\n",
    "    fig.subplots_adjust(left=0.05, right=0.9)\n",
    "    xmin, xmax = [int(x) for x in [image.xmin(), image.xmax()]]\n",
    "    ymin, ymax = [int(y) for y in [image.ymin(), image.ymax()]]\n",
    "    if nsamples == '600':\n",
    "        # set to False if you want colorbar with 3rd inage\n",
    "        removeColorbar = True\n",
    "    else:\n",
    "        removeColorbar = True\n",
    "    drawImage2D(image, excludedVal=EXCLUDED_VAL,\n",
    "                title=title,\n",
    "                removeColorbar=removeColorbar,\n",
    "                categ=True,\n",
    "                categColbad='white',\n",
    "                categCol=COLOR_SCHEME_ROUSSILLON,\n",
    "                cticklabels=LEGEND,\n",
    "                title_fontsize=FONT_SIZE+2,\n",
    "                cticklabels_fontsize=FONT_SIZE,\n",
    "                xlabels_fontsize=FONT_SIZE,\n",
    "                ylabels_fontsize=FONT_SIZE,\n",
    "                xticklabels = [xmin, xmax],\n",
    "                yticklabels = [ymin, ymax],\n",
    "                xticklabels_fontsize=FONT_SIZE,\n",
    "                yticklabels_fontsize=FONT_SIZE,\n",
    "                xticks=[xmin, xmax],\n",
    "                yticks=[ymin, ymax],\n",
    "                ylabel_rotation=0,\n",
    "               )\n",
    "    plt.scatter(point_set_roussillon.x(), point_set_roussillon.y(), marker= 'o', s=5, c='black')\n",
    "    plt.savefig(filename, dpi=DPI, bbox_inches=\"tight\")\n",
    "    #!convert -trim $filename $filename\n",
    "    !pdfcrop $filename $filename"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}