{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "# import cartopy.feature as cf\n",
    "import matplotlib.pyplot as plt\n",
    "# import numpy as np\n",
    "\n",
    "# import matplotlib.pyplot as plt\n",
    "# from matplotlib import cm\n",
    "import numpy as np\n",
    "import h5py as h5py\n",
    "# import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "from cartopy.io.shapereader import Reader\n",
    "from cartopy.feature import ShapelyFeature\n",
    "\n",
    "%matplotlib notebook\n",
    "\n",
    "projection = ccrs.PlateCarree()\n",
    "\n",
    "fig = plt.figure(dpi=150)\n",
    "\n",
    "ax = plt.axes(projection = projection)\n",
    "\n",
    "ax.coastlines('10m')\n",
    "ax.gridlines()\n",
    "# ax.add_feature(cf.BORDERS)\n",
    "# ax.add_feature(cf.STATES.with_scale('10m'))\n",
    "\n",
    "# The extent bounds are specified as an array [[x0, y0], [x1, y1]], \n",
    "# where x0 is the left side of the extent, y0 is the top, x1 is the right and y1 is the bottom.\n",
    "# extent (x0, x1, y0, y1)\n",
    "# extent = [-180,180, -90,90] # world\n",
    "# extent = [-90, -30, 20, -60] # south america\n",
    "# extent = [-74, -31, 5.5, -33] # brazil\n",
    "# extent = [-53.5, -45, -11, -20] # brazil\n",
    "# extent = [-48, -47, -16, -17] # brazil\n",
    "extent = [-44.996, -44.009, -2.805, -1.809] # brazil\n",
    "# extent = [-90, -30, 20, -60]\n",
    "# extent = [-100, 30, 0, 80]\n",
    "# extent = [-87.35, -79.5, 24.1, 30.8]\n",
    "\n",
    "\n",
    "ax.set_extent(extent, crs=ccrs.PlateCarree())\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First: read the geotiff image with GDAL.\n",
    "from osgeo import gdal, osr\n",
    "\n",
    "gdal.UseExceptions()\n",
    "\n",
    "fname = '/notebooks/resources/T23MNT_20190525T132241_TCI_60m.jp2'\n",
    "\n",
    "ds = gdal.Open(fname)\n",
    "data = ds.ReadAsArray()\n",
    "gt = ds.GetGeoTransform()\n",
    "proj = ds.GetProjection()\n",
    "\n",
    "inproj = osr.SpatialReference()\n",
    "inproj.ImportFromWkt(proj)\n",
    "\n",
    "print('\\n\\n## ds ##:\\n\\n' + str(ds))\n",
    "print('\\n\\n## data ##:\\n\\n' + str(data))\n",
    "print('\\n\\n## gt ##:\\n\\n' + str(gt))\n",
    "print('\\n\\n## proj ##:\\n\\n' + str(proj))\n",
    "print('\\n\\n## inproj ##:\\n\\n' + str(inproj))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projcs = inproj.GetAuthorityCode('PROJCS')\n",
    "print('\\n\\n## projcs ##:\\n\\n' + str(projcs))\n",
    "\n",
    "image_projection = ccrs.epsg(projcs)\n",
    "print('\\n\\n## image_projection ##:\\n\\n' + str(image_projection))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = plt.imread(fname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import cartopy.crs as ccrs\n",
    "# import cartopy.feature as cf\n",
    "# import matplotlib.pyplot as plt\n",
    "# %matplotlib inline\n",
    "\n",
    "projection = ccrs.PlateCarree()\n",
    "\n",
    "fig = plt.figure(dpi=100)\n",
    "\n",
    "ax = plt.axes(projection = projection)\n",
    "\n",
    "ax.coastlines('10m')\n",
    "ax.gridlines()\n",
    "\n",
    "# The extent bounds are specified as an array [[x0, y0], [x1, y1]], \n",
    "# where x0 is the left side of the extent, y0 is the top, x1 is the right and y1 is the bottom.\n",
    "extent = [-44.996, -44.009, -2.805, -1.809] # são luíz do maranhão\n",
    "\n",
    "ax.set_extent(extent, crs=ccrs.PlateCarree())\n",
    "img = plt.imread(fname)\n",
    "ax.imshow(img, extent=extent, origin='upper', transform=ccrs.PlateCarree())\n",
    "\n",
    "# plt.show()\n",
    "\n",
    "print('\\n\\n## get_extent() ##:\\n\\n' + str(ax.get_extent()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import matplotlib.pyplot as plt\n",
    "# from matplotlib import cm\n",
    "# import numpy as np\n",
    "# import h5py as h5py\n",
    "# import cartopy.crs as ccrs\n",
    "# import cartopy.feature as cfeature\n",
    "# from cartopy.io.shapereader import Reader\n",
    "# from cartopy.feature import ShapelyFeature\n",
    "# %matplotlib inline\n",
    "\n",
    "hdf5 = '/notebooks/resources/3B-MO.MS.MRG.3IMERG.20190101-S000000-E235959.01.V06A.HDF5'\n",
    "dataset = h5py.File(hdf5,'r')\n",
    "\n",
    "precip = dataset['Grid/precipitation'][:]\n",
    "precip = np.transpose(precip[0])\n",
    "\n",
    "theLats = dataset['Grid/lat'][:]\n",
    "theLons = dataset['Grid/lon'][:]\n",
    "\n",
    "# [-44.996, -44.009, -2.805, -1.809]\n",
    "\n",
    "condition = ((theLats < -44) & (theLats > -45))\n",
    "extractedLats = np.extract(condition, theLats)\n",
    "\n",
    "# print(theLats[(theLats < -44) & (theLats > -45)])\n",
    "# print(theLats[condition])\n",
    "# print(extractedLats)\n",
    "\n",
    "condition = ((theLons < -1.809) & (theLons > -2.805))\n",
    "extractedLons = np.extract(condition, theLons)\n",
    "\n",
    "# print(theLons[(theLons < -1.809) & (theLons > -2.805)])\n",
    "# print(theLons[condition])\n",
    "# print(extractedLons)\n",
    "\n",
    "clevs = np.arange(0,1.26,0.125)\n",
    "\n",
    "x, y = np.float32(np.meshgrid(theLons, theLats))\n",
    "# x, y = np.float32(np.meshgrid(extractedLons, extractedLats))\n",
    "\n",
    "masked_array = np.ma.masked_where(precip < 0,precip)\n",
    "\n",
    "cmap = 'nipy_spectral'\n",
    "\n",
    "plt.title('January 2019 Monthly Average Rain Rate')\n",
    "\n",
    "font = {'weight' : 'bold', 'size' : 3}\n",
    "\n",
    "plt.rc('font', **font)\n",
    "\n",
    "plt.xlim(-2.805, -1.809)\n",
    "plt.ylim(-45, -44)\n",
    "\n",
    "# cs = ax.contourf(np.reshape(x, precip), np.reshape(y, precip), precip, clevs, transform=ccrs.PlateCarree(), cmap=cmap)\n",
    "cs = ax.contourf(x, y, precip, clevs, transform=ccrs.PlateCarree(), cmap=cmap)\n",
    "\n",
    "sm = plt.cm.ScalarMappable(cmap=cmap,norm=plt.Normalize(0,1))\n",
    "sm._A = []\n",
    "plt.colorbar(sm, ax=ax, label='mm/h', shrink=0.5)\n",
    "\n",
    "# plt.show()"
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}