07_rep_modis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "import os\n",
    "from glob import glob\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024001.h13v10.061.2024019001517.hdf',\n",
       " '/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024009.h13v10.061.2024027091412.hdf',\n",
       " '/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024017.h13v10.061.2024034011201.hdf',\n",
       " '/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024025.h13v10.061.2024041004244.hdf',\n",
       " '/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024033.h13v10.061.2024049000729.hdf']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fo = '/home/geodata/Raster/modis/MYD16A2/raw/'\n",
    "ff = sorted(glob(fo+'*.hdf'))\n",
    "ff[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024001.h13v10.061.2024019001517.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024009.h13v10.061.2024027091412.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024017.h13v10.061.2024034011201.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024025.h13v10.061.2024041004244.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024033.h13v10.061.2024049000729.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024041.h13v10.061.2024061020456.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024049.h13v10.061.2024066061718.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024057.h13v10.061.2024073143618.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024065.h13v10.061.2024086100101.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024065.h13v10.061.2024088224613.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024073.h13v10.061.2024094180246.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024073.h13v10.061.2024099150458.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024081.h13v10.061.2024097021437.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024081.h13v10.061.2024099165306.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024089.h13v10.061.2024104235232.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024097.h13v10.061.2024116112327.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024105.h13v10.061.2024121001923.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024113.h13v10.061.2024130230347.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "/home/geodata/Raster/modis/MYD16A2/raw/MYD16A2.A2024121.h13v10.061.2024137000530.hdf\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n",
      "Input file size is 3304, 2624\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    }
   ],
   "source": [
    "ofo = '/home/geodata/Clientes/0FARMS/FabioRicardi-Barreiras_BA/modis/'\n",
    "#extent = '-46.200579418 -11.200452249 -44.927499105 -12.210613802' # around the Bahia farm\n",
    "extent = '-45.6861267089844034 -11.4570760726928995 -45.3921318054199006 -11.7891721725463992'\n",
    "\n",
    "Path(ofo).mkdir( parents = True, exist_ok = True)\n",
    "\n",
    "\n",
    "kwargs = {'format': 'GTiff', 'dstSRS': 'EPSG:4326'}\n",
    "for f in ff:\n",
    "    print(f)\n",
    "    fn = f.split('.')[-5][1:]+'.tif'\n",
    "    dataset = gdal.Open(f,gdal.GA_ReadOnly)\n",
    "\n",
    "    # et and pet\n",
    "    ds_et =  gdal.Open(dataset.GetSubDatasets()[0][0], gdal.GA_ReadOnly)\n",
    "    ds_pet =  gdal.Open(dataset.GetSubDatasets()[2][0], gdal.GA_ReadOnly)\n",
    "\n",
    "    # gdalwarp and clip\n",
    "    ds = gdal.Warp(destNameOrDestDS=f'{ofo}Xet_{fn}',outputType=gdal.GDT_UInt16, srcDSOrSrcDSTab=ds_et, **kwargs)\n",
    "    os.system(f'gdal_translate -projwin {extent} -of GTiff {ofo}Xet_{fn} {ofo}et_{fn}')\n",
    "    os.system(f'rm {ofo}Xet_{fn}')\n",
    "    del ds\n",
    "\n",
    "    ds = gdal.Warp(destNameOrDestDS=f'{ofo}Xpet_{fn}',outputType=gdal.GDT_UInt16, srcDSOrSrcDSTab=ds_pet, **kwargs)\n",
    "    os.system(f'gdal_translate -projwin {extent} -of GTiff {ofo}Xpet_{fn} {ofo}pet_{fn}')\n",
    "    os.system(f'rm {ofo}Xpet_{fn}')\n",
    "    del ds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ArrayToTif(Array, Filename, Folder, Metadata, Type = 1):\n",
    "    '''\n",
    "    Converts a array like to .tif file.\n",
    "    Array: Single or Multiband array\n",
    "    Folder: Destination folder of files. Ex. /tmp/\n",
    "    Filename: Name of output tif. Ex. MyTifFile\n",
    "    Metadata: Dictionary with keys: GeoTrasform, projection, rows, cols of original\n",
    "    raster.\n",
    "    Type: 1 = gdal.GDT_Int16\n",
    "          2 = gdal.GDT_Int32\n",
    "          3 = gdal.GDT_Float32\n",
    "          4 = gdal.GDT_Float64\n",
    "    '''\n",
    "    driver = gdal.GetDriverByName(\"GTiff\")\n",
    "    if not os.path.exists(Folder): os.makedirs(Folder)\n",
    "    geot = Metadata['Geo']\n",
    "    proj = Metadata['proj']\n",
    "    rows = Metadata['rows']\n",
    "    cols = Metadata['cols']\n",
    "    if Type == 0:\n",
    "        TypeIF = gdal.GDT_Byte\n",
    "    if Type == 1:\n",
    "        TypeIF = gdal.GDT_Int16\n",
    "    if Type == 2:\n",
    "        TypeIF = gdal.GDT_Int32\n",
    "    if Type == 3:\n",
    "        TypeIF = gdal.GDT_Float32\n",
    "    if Type == 4:\n",
    "        TypeIF = gdal.GDT_Float64\n",
    "    if Type == 5:\n",
    "        TypeIF = gdal.GDT_UInt16\n",
    "    if Array.ndim == 2:\n",
    "        Raster = driver.Create(Folder + '/' + Filename, cols, rows, 1, TypeIF)\n",
    "        Raster.GetRasterBand(1).WriteArray(Array)\n",
    "    else:\n",
    "        Raster = driver.Create(Folder + '/' + Filename, cols, rows, Array.shape[2], TypeIF, [ '-co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=9 -multi' ])\n",
    "        for band in range(Array.shape[2]):\n",
    "            print('Writing band' +str(band+1)+'of '+Filename)\n",
    "            Raster.GetRasterBand(band+1).WriteArray(Array[:,:,band])\n",
    "    Raster.SetGeoTransform(geot)\n",
    "    Raster.SetProjection(proj)\n",
    "    Raster.FlushCache()\n",
    "    print('%s saved.\\n '%(Filename))\n",
    "    del geot, proj, rows, cols, Raster\n",
    "    return\n",
    "\n",
    "\n",
    "def TifToArray(Raster):\n",
    "    '''\n",
    "    This function opens a Raster data as Array to use numpy manipulation\n",
    "    Return: Array, Metadata\n",
    "    Array: A array like variable\n",
    "    Metadata: Metadata of Raster file. Could be needed for a future\\n\n",
    "    ArrayToRaster transformation.\n",
    "    '''\n",
    "\n",
    "    openGdal = gdal.Open(Raster)\n",
    "    geot = openGdal.GetGeoTransform()\n",
    "    proj = openGdal.GetProjection()\n",
    "    rows = openGdal.RasterYSize\n",
    "    cols = openGdal.RasterXSize\n",
    "    Metadata = {'Geo':geot, 'proj':proj, 'rows':rows, 'cols':cols}\n",
    "    Array = openGdal.ReadAsArray()\n",
    "    return Array, Metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geo",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}