01_GB_Tensors_20210317.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create tensors sampled from Sentinel 2\n",
    "\n",
    "This notebook generates a tensor set from set of points on a Sentinel satellite manifest. It is implemented using cuDF, cuPy and numba and requires a GPU to run.\n",
    "\n",
    "The NVIDIA Docker container RAPIDS 21.02 was used to generate a computational environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cudf\n",
    "from numba import cuda\n",
    "import cupy as cp\n",
    "import numpy as np\n",
    "import gdal\n",
    "import time\n",
    "import csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Disable memory pool for device memory (GPU)\n",
    "cp.cuda.set_allocator(None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = 3\n",
    "clip = False\n",
    "suffix = str(factor)\n",
    "if clip:\n",
    "    suffix += 'c'\n",
    "print(\"Factor:\",factor,\"clip:\",clip,\"suffix:\",suffix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'onspd_sentinel_20210209_tiles.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df = cudf.read_csv(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert filenames to categories for query\n",
    "df['TCI_CAT'] = df['TCI_FILENAME'].astype('category').cat.codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class cp_manifest:\n",
    "    \n",
    "    def __init__(self, src_filename):        # Instantiation - pass database\n",
    "        self.tci_filename = src_filename\n",
    "        self.nir_filename = src_filename.replace('TCI','B08')\n",
    "        self.red_filename = src_filename.replace('TCI','B04')\n",
    "        self.grn_filename = src_filename.replace('TCI','B03')\n",
    "        self.blu_filename = src_filename.replace('TCI','B02')\n",
    "        self.scl_filename = src_filename.replace('TCI_10m','SCL_20m')\n",
    "        self.fci_filename = src_filename.replace('TCI','FCI').replace('jp2','TIFF')\n",
    "        self.rbg_filename = src_filename.replace('TCI','RGB').replace('jp2','TIFF')\n",
    "        self.nvi_filename = src_filename.replace('TCI','NVI').replace('jp2','TIFF')\n",
    "        self.nvi_filename = src_filename.replace('TCI_10m','SCL_20m').replace('jp2','TIFF')\n",
    "        self.npz_filename = src_filename.replace('TCI','NPZ').replace('jp2','npz')\n",
    "        src_ds = gdal.Open(src_filename)\n",
    "        self.metadata = src_ds.GetMetadata()\n",
    "        self.geotransform = src_ds.GetGeoTransform()\n",
    "        self.pixel_size = (self.geotransform[1], self.geotransform[5])\n",
    "        self.tile_size = (src_ds.RasterXSize, src_ds.RasterYSize)\n",
    "        self.origin = cp.array((self.geotransform[0], self.geotransform[3] + self.tile_size[1] * self.pixel_size[1]),dtype=cp.float64)\n",
    "        src_ds = None\n",
    "\n",
    "    def gen_opt_bands(self,type=\"uint8\",red_threshold=4500,grn_threshold=4500,blu_threshold=4500,nir_threshold=6500):\n",
    "        \n",
    "        bands = []\n",
    "        \n",
    "        ds = gdal.Open(self.scl_filename)\n",
    "        scl = cp.array(ds.GetRasterBand(1).ReadAsArray())\n",
    "        scl = cp.array(cp.kron(scl,cp.ones((2,2))))\n",
    "        ds = None\n",
    "        mask = scl==4\n",
    "        scl = None\n",
    "        #print(mask,mask.shape,scl.shape)\n",
    "        \n",
    "        for i, file in enumerate((self.nir_filename,self.red_filename,self.grn_filename,self.blu_filename)):\n",
    "            ds = gdal.Open(file)\n",
    "            bands += [cp.asarray(ds.GetRasterBand(1).ReadAsArray())]\n",
    "            \n",
    "            scene_median = cp.median(bands[i][bands[i]>0])\n",
    "            scene_mean = cp.mean(bands[i][bands[i]>0])\n",
    "            scene_std = cp.std(bands[i][bands[i]>0])\n",
    "            \n",
    "            mask_median = cp.median(bands[i][mask])\n",
    "            mask_mean = cp.mean(bands[i][mask])\n",
    "            mask_std = cp.std(bands[i][mask])\n",
    "\n",
    "            pc99_scene = cp.percentile(bands[i][bands[i]>0],99)\n",
    "            pc99_mask = cp.percentile(bands[i][mask],99)\n",
    "            pc991_scene = cp.percentile(bands[i][bands[i]>0],99.1)\n",
    "            pc991_mask = cp.percentile(bands[i][mask],99.1)            \n",
    "            pc999_scene = cp.percentile(bands[i][bands[i]>0],99.9)\n",
    "            pc999_mask = cp.percentile(bands[i][mask],99.9)\n",
    "           \n",
    "            ds = None\n",
    "            \n",
    "            with open('spectral_analysis.csv', 'a') as x:\n",
    "                csv_out = csv.writer(x)\n",
    "                csv_out.writerow([file, i, scene_median, mask_median, scene_mean, mask_mean, scene_std, mask_std, pc99_scene, pc99_mask, pc991_scene, pc991_mask, pc999_scene, pc999_mask])\n",
    "        \n",
    "        self.bands = cp.stack(bands)\n",
    "        \n",
    "        start = time.time()\n",
    "        self.channels =  cp.stack((cp.clip((self.bands[0]*(256/nir_threshold)),0,255).astype(cp.uint8),\n",
    "                         cp.clip((self.bands[1]*(256/red_threshold)),0,255).astype(cp.uint8),\n",
    "                         cp.clip((self.bands[2]*(256/grn_threshold)),0,255).astype(cp.uint8),\n",
    "                         cp.clip((self.bands[3]*(256/blu_threshold)),0,255).astype(cp.uint8)))\n",
    "        print(\"Computation time:\",(time.time()-start)*1000,\"ms\")\n",
    "        \n",
    "        #self.ndvi_bands = (self.bands[0].astype(cp.float32) - self.bands[1].astype(cp.float32)) / (self.bands[0].astype(cp.float32) + self.bands[1].astype(cp.float32))\n",
    "        \n",
    "    def gen_auto_bands(self,factor,clip=False):\n",
    "        \n",
    "        bands = []\n",
    "        band_factor = [1.0,factor,factor,factor]\n",
    "        threshold = cp.zeros(4,dtype=cp.uint16)\n",
    "        \n",
    "        ds = gdal.Open(self.scl_filename)\n",
    "        scl = cp.array(ds.GetRasterBand(1).ReadAsArray())\n",
    "        scl = cp.array(cp.kron(scl,cp.ones((2,2))))\n",
    "        ds = None\n",
    "        mask = scl==4\n",
    "        scl = None\n",
    "        #print(mask,mask.shape,scl.shape)\n",
    "        \n",
    "        for i, file in enumerate((self.nir_filename,self.red_filename,self.grn_filename,self.blu_filename)):\n",
    "            ds = gdal.Open(file)\n",
    "            bands += [cp.asarray(ds.GetRasterBand(1).ReadAsArray())]\n",
    "\n",
    "            threshold[i] = cp.percentile(bands[i][mask],99.5) * band_factor[i]\n",
    "                \n",
    "            print(\"Band:\",i,\"threshold:\",threshold[i])\n",
    "           \n",
    "            ds = None\n",
    "        \n",
    "        self.bands = cp.stack(bands)\n",
    "        \n",
    "        start = time.time()\n",
    "        \n",
    "        if clip:\n",
    "            self.channels =  cp.stack([cp.clip((self.bands[0].astype(cp.float32) / threshold[0]).astype(cp.float32),0,1),\n",
    "                                       cp.clip((self.bands[1].astype(cp.float32) / threshold[1]).astype(cp.float32),0,1),\n",
    "                                       cp.clip((self.bands[2].astype(cp.float32) / threshold[2]).astype(cp.float32),0,1),\n",
    "                                       cp.clip((self.bands[3].astype(cp.float32) / threshold[3]).astype(cp.float32),0,1)])\n",
    "        else:\n",
    "            self.channels =  cp.stack([(self.bands[0].astype(cp.float32) / threshold[0]).astype(cp.float32),\n",
    "                                       (self.bands[1].astype(cp.float32) / threshold[1]).astype(cp.float32),\n",
    "                                       (self.bands[2].astype(cp.float32) / threshold[2]).astype(cp.float32),\n",
    "                                       (self.bands[3].astype(cp.float32) / threshold[3]).astype(cp.float32)])            \n",
    "        \n",
    "        print(\"Computation time:\",(time.time()-start)*1000,\"ms\")\n",
    "        \n",
    "        #self.ndvi_bands = (self.bands[0].astype(cp.float32) - self.bands[1].astype(cp.float32)) / (self.bands[0].astype(cp.float32) + self.bands[1].astype(cp.float32))\n",
    "\n",
    "    def gen_ndvi_image(self,filename=None):\n",
    "                \n",
    "        if filename is None:\n",
    "            filename = self.nvi_filename\n",
    "            \n",
    "        ramp = cp.stack((cp.hstack((cp.ones(128,dtype=cp.uint8)*255,cp.arange(255,0,-257/128).astype(cp.uint8))),\n",
    "                         cp.hstack((cp.arange(0,255,2).astype(cp.uint8),cp.arange(255,127,-1).astype(cp.uint8))),\n",
    "                         cp.zeros(256,dtype=cp.uint8)))\n",
    "\n",
    "        src_ds = gdal.Open(src_filename)\n",
    "        driver = gdal.GetDriverByName(\"GTiff\")\n",
    "        nvi_ds = driver.CreateCopy(filename, src_ds, 0)\n",
    "        #print(nvi_ds)\n",
    "        #print(self.channels.shape,ramp.shape,out_bands.shape)\n",
    "        start = time.time()\n",
    "        out_bands = self.gen_ndvi_bands(self.channels,ramp)\n",
    "        print(\"Computation time:\",(time.time()-start),\"s\")\n",
    "        print(out_bands.shape)\n",
    "        print(out_bands)\n",
    "        \n",
    "        for i in range(3):\n",
    "            nvi_ds.GetRasterBand(i+1).WriteArray(cp.asnumpy(out_bands[i]))\n",
    "\n",
    "        #print(out_bands[0])\n",
    "        #print(self.bands[0])\n",
    "        #print(self.bands[0])\n",
    "        nvi_ds.FlushCache()\n",
    "        nvi_ds = None\n",
    "        src_ds = None\n",
    "        \n",
    "    def gen_ndvi_bands(self,channels,ramp):\n",
    "        return ramp[:,cp.nan_to_num((((channels[0].astype(cp.float32) - channels[1].astype(cp.float32)) / \n",
    "                                      (channels[0].astype(cp.float32) + channels[1].astype(cp.float32))) / \n",
    "                                     2.0 + 0.5) * 255.0).astype(cp.uint8)]\n",
    "    \n",
    "    def save_bands(self):\n",
    "        cp.savez_compressed(self.npz_filename,bands=self.bands)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@cuda.jit('void(f8[:],f8[:],f8[:],f8,u2,f4[:,:,:],f4[:,:,:,:],u4[:])')\n",
    "def gen_image_tensors(east,north,origin,scale,size,channels,tensors,seq):\n",
    "    \"\"\" Generate a tensor set from set of points on a Sentinel satellite manifest .\n",
    "    inputs\n",
    "    east: UTM Easting Coordinates\n",
    "    north: UTM Northing Coordinates\n",
    "    origin: UTM Origin (tuple)\n",
    "    scale: Coordinate units per pixel\n",
    "    channels: Satellite colour channels (NIR, R, G, B)\n",
    "    seq: offset in tensor file\n",
    "    Outputs\n",
    "    tensors: A set of tensors, one per point, size x size x 4\n",
    "    Inputs\n",
    "    ------\n",
    "    {inputs}\n",
    "    Outputs\n",
    "    -------\n",
    "    {outputs}\n",
    "    \"\"\"   \n",
    "    i= cuda.grid(1)\n",
    "    if i < east.size and i < north.size:\n",
    "        # Calculate relative origins\n",
    "        c_x = int(round(((east[i] - origin[0]) / scale) - (size / 2)))\n",
    "        c_y = int(round(((north[i] - origin[1]) / scale) - (size / 2)))\n",
    "        # Counters for all black or all white pixels\n",
    "        height = channels.shape[1]\n",
    "        for x in range(size):\n",
    "            for y in range(size):\n",
    "                for j in range(4):\n",
    "                    if 0 <= height-(c_y+y)-1 < channels[j].shape[0] and 0 <= c_x+x < channels[j].shape[1]:\n",
    "                        tensors[seq[i]][size-y-1][x][j] = channels[j][height-(c_y+y)][c_x+x]\n",
    "                    else:\n",
    "                        tensors[seq[i]][size-y-1][x][j] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "img_size = 16\n",
    "#tensors = cp.zeros((len(df['PCDS']),img_size,img_size,4),dtype=cp.uint8)\n",
    "#tensors = cp.zeros((len(df['PCDS']),img_size,img_size,4),dtype=cp.float32)\n",
    "tensors = cuda.mapped_array((len(df['PCDS']),img_size,img_size,4),dtype=np.float32)\n",
    "#tensors[:] = np.zeros((len(df['PCDS']),img_size,img_size,4),dtype=np.float32)\n",
    "cp_tensors = cp.asarray(tensors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Manifests:\",df['TCI_CAT'].nunique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df['TCI_CAT'].unique().values_host)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#with open('spectral_analysis.csv', 'w') as x:\n",
    "#    csv_out = csv.writer(x)\n",
    "#    csv_out.writerow(['FILE','BAND','SCENE_MEDIAN','MASK_MEDIAN','SCENE_MEAN','MASK_MEAN','SCENE_STD','MASK_STD','SCENE_99','MASK_99','SCENE_991','MASK_991','SCENE_999','MASK_999'])\n",
    "\n",
    "for tci in df['TCI_CAT'].unique().values_host:\n",
    "    query = 'TCI_CAT=='+str(tci)\n",
    "    qdf = df.query(query)\n",
    "    print(query)\n",
    "    #print(qdf)\n",
    "    src_filename = qdf['TCI_FILENAME'].iloc[0]\n",
    "    print(tci,src_filename)\n",
    "    manifest = cp_manifest(src_filename)\n",
    "    #manifest.gen_opt_bands()\n",
    "    manifest.gen_auto_bands(factor=factor,clip=clip)\n",
    "    cuda.synchronize()\n",
    "    #manifest.save_bands()\n",
    "    #print(manifest.channels,manifest.channels.shape,manifest.origin,manifest.pixel_size[0],tensors.shape)\n",
    "    #print(qdf['UTM_X'].shape,qdf['UTM_Y'].shape,manifest.origin,manifest.pixel_size[0],img_size,manifest.channels.shape,tensors.shape,qdf['SEQ'].shape)\n",
    "    # Generate Keras Tensor Set\n",
    "    gen_image_tensors.forall(len(qdf['PCDS']))(qdf['UTM_X'],qdf['UTM_Y'],manifest.origin,manifest.pixel_size[0],img_size,manifest.channels,tensors,qdf['SEQ'])\n",
    "    cuda.synchronize()\n",
    "    manifest = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nvidia-smi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "cp.savez_compressed(\"sentinel_gb_995_\"+suffix+\"_20210317\",tensors=cp_tensors)\n",
    "print(\"File sentinel_gb_995_\"+suffix+\"_20210317.npz written\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}