#! /usr/bin/env python3

import os
import pygrib
import matplotlib.pyplot as plt
import cartopy, cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import numpy as np
from cartopy.util import add_cyclic_point

grib = pygrib.open('/lustre_xc50/ioper/models/SMNA-Oper/SMG/datainout/bam/pos/dataout/TQ0299L064/2025100100/GPOSCPT20251001002025100400P.fct.TQ0299L064.grb')

nx = 900
lon0 = 0.0
dx = 0.4
lons = lon0 + dx * np.arange(nx)

lats = np.array([
-89.69415, -89.29794, -88.89940, -88.50032, -88.10105, -87.70167, -87.30225, -86.90279,
-86.50331, -86.10381, -85.70430, -85.30479, -84.90526, -84.50574, -84.10621, -83.70667,
-83.30714, -82.90760, -82.50806, -82.10851, -81.70897, -81.30943, -80.90988, -80.51033,
-80.11079, -79.71124, -79.31169, -78.91214, -78.51259, -78.11304, -77.71349, -77.31394,
-76.91439, -76.51484, -76.11528, -75.71573, -75.31618, -74.91663, -74.51708, -74.11752,
-73.71797, -73.31842, -72.91886, -72.51931, -72.11976, -71.72020, -71.32065, -70.92110,
-70.52154, -70.12199, -69.72244, -69.32288, -68.92333, -68.52377, -68.12422, -67.72466,
-67.32511, -66.92556, -66.52600, -66.12645, -65.72689, -65.32734, -64.92778, -64.52823,
-64.12867, -63.72912, -63.32956, -62.93001, -62.53045, -62.13090, -61.73134, -61.33179,
-60.93223, -60.53268, -60.13312, -59.73357, -59.33401, -58.93446, -58.53490, -58.13535,
-57.73579, -57.33624, -56.93668, -56.53713, -56.13757, -55.73802, -55.33846, -54.93891,
-54.53935, -54.13980, -53.74024, -53.34069, -52.94113, -52.54157, -52.14202, -51.74246,
-51.34291, -50.94335, -50.54380, -50.14424, -49.74469, -49.34513, -48.94558, -48.54602,
-48.14646, -47.74691, -47.34735, -46.94780, -46.54824, -46.14869, -45.74913, -45.34958,
-44.95002, -44.55046, -44.15091, -43.75135, -43.35180, -42.95224, -42.55269, -42.15313,
-41.75358, -41.35402, -40.95446, -40.55491, -40.15535, -39.75580, -39.35624, -38.95669,
-38.55713, -38.15757, -37.75802, -37.35846, -36.95891, -36.55935, -36.15980, -35.76024,
-35.36069, -34.96113, -34.56157, -34.16202, -33.76246, -33.36291, -32.96335, -32.56380,
-32.16424, -31.76468, -31.36513, -30.96557, -30.56602, -30.16646, -29.76691, -29.36735,
-28.96779, -28.56824, -28.16868, -27.76913, -27.36957, -26.97002, -26.57046, -26.17090,
-25.77135, -25.37179, -24.97224, -24.57268, -24.17313, -23.77357, -23.37401, -22.97446,
-22.57490, -22.17535, -21.77579, -21.37623, -20.97668, -20.57712, -20.17757, -19.77801,
-19.37846, -18.97890, -18.57934, -18.17979, -17.78023, -17.38068, -16.98112, -16.58157,
-16.18201, -15.78245, -15.38290, -14.98334, -14.58379, -14.18423, -13.78468, -13.38512,
-12.98556, -12.58601, -12.18645, -11.78690, -11.38734, -10.98778, -10.58823, -10.18867,
-9.78912,  -9.38956, -8.99001,  -8.59045,  -8.19089,  -7.79134, -7.39178,  -6.99223,
-6.59267,  -6.19311,  -5.79356,  -5.39400,  -4.99445,  -4.59489,  -4.19534,  -3.79578,
-3.39622,  -2.99667,  -2.59711,  -2.19756,  -1.79800,  -1.39845,  -0.99889,  -0.59933,
-0.19978,   0.19978,   0.59933,   0.99889,   1.39845,   1.79800,   2.19756,   2.59711,
2.99667,   3.39622,  3.79578,   4.19534,   4.59489,   4.99445,   5.39400,   5.79356,
6.19311,   6.59267,   6.99223,   7.39178,   7.79134,   8.19089,   8.59045,   8.99001,
9.38956,   9.78912, 10.18867,  10.58823,  10.98778,  11.38734,  11.78690, 12.18645,
12.58601,  12.98556,  13.38512,  13.78468,  14.18423,  14.58379,  14.98334,  15.38290,
15.78245, 16.18201,  16.58157, 16.98112,  17.38068,  17.78023,  18.17979,  18.57934,
18.97890,  19.37846,  19.77801,  20.17757,  20.57712,  20.97668,  21.37623,  21.77579,
22.17535,  22.57490,  22.97446,  23.37401,  23.77357,  24.17313,  24.57268,  24.97224,
25.37179,  25.77135,  26.17090,  26.57046,  26.97002,  27.36957,  27.76913,  28.16868,
28.56824,  28.96779,  29.36735,  29.76691,  30.16646,  30.56602,  30.96557,  31.36513,
31.76468,  32.16424,  32.56380,  32.96335,  33.36291,  33.76246,  34.16202,  34.56157,
34.96113,  35.36069,  35.76024,  36.15980,  36.55935,  36.95891,  37.35846,  37.75802,
38.15757,  38.55713,  38.95669,  39.35624,  39.75580,  40.15535,  40.55491,  40.95446,
41.35402,  41.75358,  42.15313,  42.55269,  42.95224,  43.35180,  43.75135,  44.15091,
44.55046,  44.95002,  45.34958,  45.74913,  46.14869,  46.54824,  46.94780,  47.34735,
47.74691,  48.14646,  48.54602,  48.94558,  49.34513,  49.74469,  50.14424,  50.54380,
50.94335,  51.34291,  51.74246,  52.14202,  52.54157,  52.94113,  53.34069,  53.74024,
54.13980,  54.53935,  54.93891,  55.33846,  55.73802,  56.13757,  56.53713,  56.93668,
57.33624,  57.73579,  58.13535,  58.53490,  58.93446,  59.33401,  59.73357,  60.13312,
60.53268,  60.93223,  61.33179,  61.73134,  62.13090,  62.53045,  62.93001,  63.32956,
63.72912,  64.12867,  64.52823,  64.92778,  65.32734,  65.72689,  66.12645,  66.52600,
66.92556,  67.32511,  67.72466,  68.12422,  68.52377,  68.92333,  69.32288,  69.72244,
70.12199,  70.52154,  70.92110,  71.32065,  71.72020,  72.11976,  72.51931,  72.91886,
73.31842,  73.71797,  74.11752,  74.51708,  74.91663,  75.31618,  75.71573,  76.11528,
76.51484,  76.91439,  77.31394,  77.71349,  78.11304,  78.51259,  78.91214,  79.31169,
79.71124,  80.11079,  80.51033,  80.90988,  81.30943,  81.70897,  82.10851,  82.50806,
82.90760,  83.30714,  83.70667,  84.10621,  84.50574,  84.90526,  85.30479,  85.70430,
86.10381,  86.50331,  86.90279,  87.30225,  87.70167,  88.10105,  88.50032,  88.89940,
89.29794,  89.69415])

#grb = grib.select(name='Time mean surface relative humidity',level=3)[0]
grb = grib.select(name='Time mean surface relative humidity',level=3)[0]

init  = str(grb.analDate)      # Init date / time
run   = str(grb.hour).zfill(2) # Run
valid = str(grb.validDate)     # Valid date / time

tmtmp = grb.values 

tmtmp_plot = np.flipud(tmtmp)

# remove linha branca vertical no meio da figura
data, lons = add_cyclic_point(tmtmp_plot, coord=lons)

plt.figure(figsize=(30,10))
ax = plt.axes(projection=ccrs.PlateCarree())

lon2d, lat2d = np.meshgrid(lons, lats)

#data_min = 600
#data_max = 1000
#interval = 50
#levels = np.arange(data_min, data_max+interval, interval)

vmin = np.nanmin(data)
vmax = np.nanmax(data)
interval = 50
levels = np.linspace(vmin, vmax, interval)

img1 = ax.contourf(lon2d, lat2d, data, cmap='jet', levels=levels, extend='both')

#img2 = ax.contour(lon2d, lat2d, data, colors='white', linewidths=0.3, levels=levels)
#ax.clabel(img2, inline=1, inline_spacing=15, fontsize=10, fmt='%1.0f', colors='white')

shapefile = list(shpreader.Reader('/lustre_xc50/carlos_bastarz/SMNAMonitoringApp/cron_scripts/anls_img/br_unidades_da_federacao/BR_UF_2019.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='black', facecolor='none', linewidth=0.3)

ax.coastlines(resolution='10m', color='black', linewidth=0.8)
ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(), color='gray', alpha=1.0, linestyle='--', linewidth=0.25,
xlocs=np.arange(-180, 180, 5), ylocs=np.arange(-90, 90, 5), draw_labels=True)
gl.top_labels = False
gl.right_labels = False

plt.colorbar(img1, label='Time mean surface relative humidity', orientation='vertical', pad=0.01, fraction=0.05)

plt.title('SMNA: Time mean surface relative humidity @ 3 hPa, FCT 72 h', fontweight='bold', fontsize=12, loc='left')
plt.title('Valid: ' + valid, fontsize=12, loc='right')

fig_dir = os.path.join('/lustre_xc50/carlos_bastarz/SMNAMonitoringApp/cron_scripts/anls_img/imgs', 'SMNA', str(2025100100), str(grb.name).replace(' ','_').replace('(','').replace(')',''), str(3))

os.makedirs(fig_dir, exist_ok=True)

#fname_fig = 'SMNA_' + str(grb.name).replace(' ','_').replace('(','').replace(')','') + '_level_' + str(3) + 'hPa_valid_for_' + str(valid).replace(' ','_').replace('-','_').replace(':00:00','Z') + '_fct_' + str(72) + 'h'
#fname_fig = 'SMNA/' + str(grb.name).replace(' ','_').replace('(','').replace(')','') + '/' + str(3) + '/' + str(72)
#fname_fig = os.path.join('/lustre_xc50/carlos_bastarz/SMNAMonitoringApp/cron_scripts/anls_img/imgs', 'SMNA', str(grb.name).replace(' ','_').replace('(','').replace(')',''), str(3), str(72))
fname_fig = os.path.join(fig_dir, str(72))

#plt.savefig(fname_fig + '.png', dpi=150, bbox_inches='tight')
plt.savefig(fname_fig + '.png', bbox_inches='tight')

grib.close()