import cartopy
import cartopy.feature as cfeature
import datetime
from goes2go.data import goes_nearesttime
import matplotlib.pyplot as plt
import os
import pandas as pd
import sys

"""
On NCAR's capser
conda activate goes2go
"""


# Assign first argument to itime. 
# Date/time format must be readable by the pandas module to_datetime() function.
# For example 20191025 or 20201001T1215 
itime = pd.to_datetime(sys.argv[1])

# Download nearest time
G = goes_nearesttime(itime, satellite='G16')

fig = plt.figure(figsize=(10,8)) # define figure size (width,height) in inches

# Initialize axes with GOES projection
ax = plt.subplot(projection=G.rgb.crs)

# state borders
ax.add_feature(cfeature.STATES)


# set TrueColor to True or False
TrueColor = True

if TrueColor:
    # gamma > 1 will lighten an image.
    # gamma < 1 will darken an image.
    gamma = 2.2 # 2.2 is default
    rgb = G.rgb.TrueColor(gamma=gamma, night_IR=True)
else:
    # if not TrueColor, choose ABI Band (1-16)
    abi_band = 14 
    rgb = G[f"CMI_C{abi_band:02d}"] # zero-padded ABI Band Number
"""
=============== ================== ============================================== ======================================
ABI Band Number Central Wavelength  Name                                          Type
=============== ================== ============================================== ======================================
      1             0.47 μm        "Blue" Band                                    Visible
      2             0.64 μm        "Red" Band                                     Visible
      3             0.86 μm        "Veggie" Band                                  Near-IR
      4             1.37 μm        "Cirrus" Band                                  Near-IR
      5             1.6  μm        "Snow/Ice" Band                                Near-IR
      6             2.2  μm        "Cloud Particle Size" Band                     Near-IR
      7             3.9  μm        "Shortwave Window" Band                        IR (with reflected daytime component)
      8             6.2  μm        "Upper-Level Tropospheric Water Vapor" Band    IR
      9             6.9  μm        "Mid-Level Tropospheric Water Vapor" Band      IR
      10            7.3  μm        "Lower-level Water Vapor" Band                 IR
      11            8.4  μm        "Cloud-Top Phase" Band                         IR
      12            9.6  μm        "Ozone Band"                                   IR
      13            10.3 μm        "Clean" IR Longwave Window Band                IR
      14            11.2 μm        IR Longwave Window Band                        IR
      15            12.3 μm        "Dirty" Longwave Window Band                   IR
      16            13.3 μm        "CO2" Longwave infrared                        IR
=============== ================== ============================================== ======================================
"""


ax.imshow(rgb, **G.rgb.imshow_kwargs)

ax.set_title(f'{G.attrs["orbital_slot"]} {rgb.long_name} ({rgb.name})')
text = G.t.dt.strftime('%H:%M UTC %d %B %Y').item()
# Write text in upper left corner of image.
ax.text(.01,0.99,text,color="white", verticalalignment="top", transform=ax.transAxes)

# Zoom in
lon0 = -98
lon1 = -80
lat0 = 30
lat1 = 48
ax.set_extent((lon0, lon1, lat0, lat1))


# labels
cities = {
        "DVN":(-90.58, 41.61),
        "TOP":(-95.62,39.07),
        }

for city in cities:
    print(city)
    ax.plot(*cities[city], 'o', markersize=1, label=city, transform=cartopy.crs.PlateCarree())
plt.legend()

# Save to png file. 
# Put type of image and time in filename. 
ofile = f'{rgb.name}.{itime.strftime("%Y%m%d_%H%Mz.png")}'
ofile = os.path.realpath(ofile) # absolute path, not relative path
plt.savefig(ofile)
print(f"created {ofile}")

# pop-up interactive window
plt.show()