01c_Label_Generation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c9452d15",
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "from osgeo import gdalconst\n",
    "import os\n",
    "from osgeo import ogr\n",
    "from osgeo import osr\n",
    "import fiona\n",
    "from ops.ops import load_json\n",
    "from tqdm.notebook import tqdm_notebook\n",
    "from osgeo import gdal_array\n",
    "from skimage.morphology import disk, dilation, erosion\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "82f4aaa7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#image_folder = 'img/GEE-imgs'\n",
    "\n",
    "paths = load_json(os.path.join('conf', 'paths.json'))\n",
    "shp = load_json(os.path.join('conf', 'shp.json'))\n",
    "conf = load_json(os.path.join('conf', 'conf.json'))\n",
    "shp_path = paths['shp']\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\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "55725037",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rasterize(img_id):\n",
    "    image_folder = os.path.join('img', 'GEE_imgs_opt')\n",
    "    labels_folder = os.path.join('img', 'labels')\n",
    "    \n",
    "    image_p = os.path.join(image_folder, f'{img_id}_opt_2019.tif')\n",
    "\n",
    "    v_def_until_2007 = ogr.Open(os.path.join(shp_path, 'def_2007_pa.shp'))\n",
    "    l_def_until_2007 = v_def_until_2007.GetLayer()\n",
    "    \n",
    "    v_hidro = ogr.Open(os.path.join(shp_path, 'hidro_pa.shp'))\n",
    "    l_hidro = v_hidro.GetLayer()\n",
    "    \n",
    "    v_no_forest = ogr.Open(os.path.join(shp_path, 'no_forest_pa.shp'))\n",
    "    l_no_forest = v_no_forest.GetLayer()\n",
    "    \n",
    "    v_def_2008_2018 = ogr.Open(os.path.join(shp_path, 'def_2008_2018_pa.shp'))\n",
    "    l_def_2008_2018 = v_def_2008_2018.GetLayer()\n",
    "    \n",
    "    v_def_2019 = ogr.Open(os.path.join(shp_path, 'def_2019_pa.shp'))\n",
    "    l_def_2019 = v_def_2019.GetLayer()\n",
    "    \n",
    "    v_def_2020 = ogr.Open(os.path.join(shp_path, 'def_2020_pa.shp'))\n",
    "    l_def_2020 = v_def_2020.GetLayer()\n",
    "    \n",
    "    v_def_2020_l = ogr.Open(os.path.join(shp_path, 'def_2020_l_pa.shp'))\n",
    "    l_def_2020_l = v_def_2020_l.GetLayer()\n",
    "    \n",
    "    data = gdal.Open(image_p, gdalconst.GA_ReadOnly)\n",
    "    \n",
    "    geo_transform = data.GetGeoTransform()\n",
    "    x_min = geo_transform[0]\n",
    "    y_max = geo_transform[3]\n",
    "    x_max = x_min + geo_transform[1] * data.RasterXSize\n",
    "    y_min = y_max + geo_transform[5] * data.RasterYSize\n",
    "    x_res = data.RasterXSize\n",
    "    y_res = data.RasterYSize\n",
    "    \n",
    "    crs = data.GetSpatialRef()\n",
    "    proj = data.GetProjection()\n",
    "\n",
    "    pixel_width = geo_transform[1]\n",
    "    output = os.path.join(labels_folder, f'label_{img_id}.tif')\n",
    "    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)\n",
    "    #target_ds.SetGeoTransform((x_min, geo_transform[1], 0, y_min, 0, geo_transform[5]))\n",
    "    target_ds.SetGeoTransform(geo_transform)\n",
    "    target_ds.SetSpatialRef(crs)\n",
    "    target_ds.SetProjection(proj)\n",
    "    band = target_ds.GetRasterBand(1)\n",
    "    #NoData_value = 0\n",
    "    #band.SetNoDataValue(NoData_value)\n",
    "    band.FlushCache()\n",
    "    \n",
    "    gdal.RasterizeLayer(target_ds, [1], l_def_until_2007, burn_values=[2])\n",
    "    gdal.RasterizeLayer(target_ds, [1], l_hidro, burn_values=[2])\n",
    "    gdal.RasterizeLayer(target_ds, [1], l_no_forest, burn_values=[2])\n",
    "    gdal.RasterizeLayer(target_ds, [1], l_def_2008_2018, burn_values=[2])\n",
    "    gdal.RasterizeLayer(target_ds, [1], l_def_2019, burn_values=[2])\n",
    "    \n",
    "    gdal.RasterizeLayer(target_ds, [1], l_def_2020, burn_values=[1])\n",
    "    rasterized_data = target_ds.ReadAsArray() \n",
    "    defor_data = rasterized_data == 1\n",
    "    defor_data = defor_data.astype(np.uint8)\n",
    "\n",
    "    border_data = dilation(defor_data, disk(1)) - erosion(defor_data, disk(1))\n",
    "    rasterized_data[border_data==1] = 2\n",
    "\n",
    "    target_ds.GetRasterBand(1).WriteArray(rasterized_data)\n",
    "\n",
    "    #gdal.RasterizeLayer(target_ds, [1], l_def_2020_l, burn_values=[2], options = ['ALL_TOUCHED=TRUE'])\n",
    "    \n",
    "    #print(target_ds.GetGeoTransform())\n",
    "    target_ds = None\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "723f1a61",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0fdc3d6bac96452685c51f732041910f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  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):\n",
    "        rasterize(int(feat['properties']['id']))\n",
    "        #print(feat['properties']['id'])"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "6984327cd4a671f1c687c03fe0b5d8b7594868635f0d5dbeb7c4ac1a044793af"
  },
  "kernelspec": {
   "display_name": "Python 3.9.12 64-bit ('ferrari': conda)",
   "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
}