﻿"""
Created on Thu Dec  6 19:21:36 2018

@author: allenea
Eric Allen, University of Delaware

Last Updated: 2/6/2020

Written in Python3(.6)
        Python/3.6
        Cartopy/0.17
        Matplotlib/3.0
CLASS: wps_info ("namelist.wps") -- > Calls get_wps_info

METHODS:
    - get_wps_info(filename,plot_info=False) -- Default
        -- Intialized when calling wps_info(filename)... Optional to print info right away

    - print_info()

    - get_proj_lambert_conformal()
        +(RETURNS: projection)

    - get_proj_mercator()
        +(RETURNS: projection)

    - get_proj_stereographic()
        +(RETURNS: projection)

    - get_proj_latlon()
        +(RETURNS: projection)

    - plot_domain_number(domain_number)  -- default -- 1
        + (RETURNS: domain_number)

    - calc_wps_domain_info() **
        - (RETURNS: wpsproj, latlonproj, corner_lat_full, corner_lon_full, length_x, length_y)

    - calc_corner_point_latlon(center_lat, center_lon, e_we, e_ns, dx, dy, wpsproj, latlonproj, loc)
       +(RETURNS: corner_lon, corner_lat)

    - calc_center_point_latlon(corner_lat_parent, corner_lon_parent,\
                    dx_parent, dy_parent, e_we, e_ns, dx, dy, i, j, wpsproj, latlonproj) **
        +(RETURNS: center_lon_child, center_lat_child)

    - reproject_corners(corner_lons, corner_lats, wpsproj, latlonproj) **
        +(RETURNS: return corner_x, corner_y)

    - plot_zoom_out(corner_x, corner_y, ns_zoom,we_zoom)
        +(RETURNS: return corner_x, corner_y)

    - plot_zoom_in(corner_x, corner_y, ns_zoom,we_zoom)
        +(RETURNS: return corner_x, corner_y)

    - plot_all_domains(ax1, fig1, wpsproj, latlonproj, corner_lat_full,\
                           corner_lon_full, length_x, length_y)
        +(RETURNS: ax1, fig1 )


OUTSIDE FUNCTION:
    - _ismissing(val, islat=True)  -- FROM WRF-PYTHON
        +(RETURNS: (bool) True/False)


** Adapted from Xiaodong Chen @lucas-uw on github

Code heavily adapted and corrected from https://github.com/lucas-uw/WRF-tools/blob/master/
            WRF_input_tools/Visualize_WPS_domain.ipynb and combined with my previous
attempts at reading/using/plotting wrf's namelist.wps file-data.


He had a great calc_wps_domain_info function.
    I modified it to work with my wps_info class replacing my reprojection code.
    - Decided not to implement the lambert projection gridlines work-around
        - Wait until cartopy implements it in a future version or
                let the user implement at their own peril

I was having trouble recalculating the centerpoint and then getting the corners
 perfect in all nested domains for all projection types.
Instead of just plotting the bounds of each domain, I have made it  so that you
 can create a cartopy projection and use only that domain,
 then from there you can choose to zoom in and out as you wish.

You do not have to use this class for plotting, but that is one of the features...
 The other main feature is it hold namelist.wps data
    - You can also use it to recreate your grid for other uses such as spatial
        interpolation or plotting observation data.


#!! Problems/Improvements:
     - Gridlines
     - X/Y Ticks and Labels for geographic reference
     - Is it possible to improve zoom or shift center point of zoom?
     - No Testing
     - How plot_all_domains applies gridlines

"""
#IMPORTS
from __future__ import print_function
import os
import os.path
#import sys
import math
from copy import copy
import matplotlib as mpl
from matplotlib.font_manager import FontProperties
import matplotlib.ticker as ticker
import numpy as np
import cartopy
import cartopy.crs as ccrs
from cartopy.feature import OCEAN, LAKES, LAND  #,NaturalEarthFeature
#from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import shapely.geometry as sgeom


class wps_info:
    """
    Developed by Eric Allen, University of Delaware

    wps = wps_info(infile)
    plt_idx = wps.plot_domain_number(plot_domains)
    wpsproj, latlonproj, corner_lat_full, corner_lon_full,
            length_x, length_y = wps.calc_wps_domain_info()

    """
    def __init__(self, file, print_info=False):
        self.filelocname = file
        self._printinfo = print_info
        if os.path.exists(self.filelocname):
            wps_info.get_wps_info(self.filelocname)
            if self._printinfo:
                wps_info.print_info()
        else:
            print("FILE NOT FOUND IN PROVIDED PATH AND/OR FILENAME")
            raise ValueError


    @classmethod
    def get_wps_info(cls, namelist_wps_filename="namelist.wps"):
        """
        Reads the namelist file and stores the important information. Default = namelist.wps

        """

        newdx = []
        newdy = []
        cls.pole_lat = 90.0 # these get changed if there is something to change it with
        cls.pole_lon = 0.0  # these get changed if there is something to change it with
        file = open(namelist_wps_filename, 'r').readlines()
        for row in file:
            row = " ".join(row.split(","))
            row = " ".join(row.split("\n"))
            row = row.lstrip()

            if len(row) != 0 and row[0] == "!":
                pass

            elif 'wrf_core' in row:
                tmprow = row.split("=")[1]
                tmprow = " ".join(tmprow.split("'"))
                cls.wrf_core = str(tmprow.strip())

            elif 'max_dom' in row:
                tmprow = row.split("=")[1]
                cls.max_dom = int(tmprow)

            elif 'start_date' in row:
                tmprow = row.split("=")[1] #All Start Dates for each domain should be the same
                tmprow = " ".join(tmprow.split("'"))
                tmprow2 = tmprow.split(" ")
                start1 = [item for item in tmprow2 if item != '']
                cls.start_date = start1[0]

            elif 'end_date' in row:
                tmprow = row.split("=")[1] #All End Dates for each domain should be the same
                tmprow = " ".join(tmprow.split("'"))
                tmprow2 = tmprow.split(" ")
                end1 = [item for item in tmprow2 if item != '']
                cls.end_date = end1[0]

            elif 'interval_seconds' in row:
                tmprow = row.split("=")[1]
                cls.interval_seconds = int(tmprow)

            elif 'parent_id' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                pid = [int(item) for item in tmprow2 if item != '']
                cls.parent_id = pid

            elif 'parent_grid_ratio' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                pgr = [int(item) for item in tmprow2 if item != '']
                cls.parent_grid_ratio = pgr

            elif 'i_parent_start' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                ipar = []
                for item in tmprow2:
                    if item != '':
                        num = int(item)
                        ipar.append(num-1)


                ipar2 = [int(item) for item in tmprow2 if item != '']
                cls.i_parent_start = ipar
                cls.i_parent_start_original = ipar2

            elif 'j_parent_start' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                jpar = []
                for item in tmprow2:
                    if item != '':
                        num = int(item)
                        jpar.append(num-1)


                jpar2 = [int(item) for item in tmprow2 if item != '']
                cls.j_parent_start = jpar
                cls.j_parent_start_original = jpar2


            elif 'e_we' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                ewe = []
                for item in tmprow2:
                    if item != '':
                        num = int(item)
                        ewe.append(num-1)


                ewe2 = [int(item) for item in tmprow2 if item != '']
                cls.e_we = ewe
                cls.e_we_original = ewe2

            elif 'e_sn' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                esn = []
                for item in tmprow2:
                    if item != '':
                        num = int(item)
                        esn.append(num-1)
                esn2 = [int(item) for item in tmprow2 if item != '']
                cls.e_sn = esn
                cls.e_sn_original = esn2


            elif 'geog_data_res' in row:
                tmprow = row.split("=")[1]
                tmprow = " ".join(tmprow.split("'"))
                tmprow2 = tmprow.split(" ")
                gdr = [item for item in tmprow2 if item != '']
                cls.geog_data_res = gdr

            elif 'dx' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                dxval = [float(item) for item in tmprow2 if item != '']
                dx1 = dxval[0]
                for dom_idx in range(cls.max_dom):
                    if dom_idx == 0:
                        tmpdx = dx1
                        newdx.append(tmpdx)
                    else:
                        tmpdx = tmpdx/cls.parent_grid_ratio[dom_idx]
                        newdx.append(tmpdx)
                cls.dx0 = newdx[0]
                cls.dx = newdx


            elif 'dy' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                dyval = [float(item) for item in tmprow2 if item != '']
                dy1 = dyval[0]
                for dom_idx in range(cls.max_dom):
                    if dom_idx == 0:
                        tmpdy = dy1
                        newdy.append(tmpdy)
                    else:
                        tmpdy = tmpdy/cls.parent_grid_ratio[dom_idx]
                        newdy.append(tmpdy)
                cls.dy0 = newdy[0]
                cls.dy = newdy

            elif 'map_proj' in row:
                tmprow = row.split("=")[1]
                tmprow = "".join(tmprow.split("'"))
                cls.map_proj = str(tmprow.strip())

            elif 'ref_lat' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                rlat = [float(item) for item in tmprow2 if item != '']
                cls.ref_lat = rlat[0]

            elif 'ref_lon' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                rlon = [float(item) for item in tmprow2 if item != '']
                cls.ref_lon = rlon[0]

            elif 'truelat1' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                tl1 = [float(item) for item in tmprow2 if item != '']
                cls.truelat1 = tl1[0]

            elif 'truelat2' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                tl2 = [float(item) for item in tmprow2 if item != '']
                cls.truelat2 = tl2[0]

            elif 'stand_lon' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                stdlon = [float(item) for item in tmprow2 if item != '']
                cls.stand_lon = stdlon[0]

            elif 'ref_x' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                refx = [float(item) for item in tmprow2 if item != '']
                cls.ref_x = refx[0]

            elif 'ref_y' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                refy = [float(item) for item in tmprow2 if item != '']
                cls.ref_y = refy[0]

            elif 'pole_lat' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                polelat = [float(item) for item in tmprow2 if item != '']
                cls.pole_lat = polelat[0]

            elif 'pole_lon' in row:
                tmprow = row.split("=")[1]
                tmprow2 = tmprow.split(" ")
                polelon = [float(item) for item in tmprow2 if item != '']
                cls.pole_lon = polelon[0]


    @classmethod
    def print_info(cls):
        """Print the output from the wps file"""
        print("wrf_core:", cls.wrf_core)
        print("max_dom:", cls.max_dom)
        print("start_date:", cls.start_date)
        print("end_date:", cls.end_date)
        print("interval_seconds:", cls.interval_seconds)
        print("parent_id:", cls.parent_id)
        print("parent_grid_ratio:", cls.parent_grid_ratio)
        print("i_parent_start:", cls.i_parent_start_original)
        print("j_parent_start:", cls.j_parent_start_original)
        print("e_we:", cls.e_we_original)
        print("e_sn:", cls.e_sn_original)
        print("geog_data_res:", cls.geog_data_res)
        print("dx0 =", cls.dx0)
        print("dx =", cls.dx)
        print("dy0 =", cls.dy0)
        print("dy =", cls.dy)
        print("map_proj:", cls.map_proj)
        print("ref_lat:", cls.ref_lat)
        print("ref_lon:", cls.ref_lon)
        print("truelat1:", cls.truelat1)
        print("truelat2:", cls.truelat2)
        print("stand_lon:", cls.stand_lon)
        print("ref_x:", cls.ref_x)
        print("ref_y:", cls.ref_y)
        print("pole_lat:", cls.pole_lat)
        print("pole_lon:", cls.pole_lon)

    #%% Projection Class-Methods
    @classmethod
    def get_proj_lambert_conformal(cls):
        """doc-string"""
        lccproj = ccrs.LambertConformal(central_longitude=cls.ref_lon,\
                            central_latitude=cls.ref_lat,\
                            standard_parallels=(cls.truelat1, cls.truelat2),\
                            globe=None, cutoff=10)
        return lccproj

    @classmethod #UNTESTED
    def get_proj_mercator(cls):
        """doc-string"""
        merproj = ccrs.Mercator(central_longitude=cls.ref_lon, globe=None)
        return merproj


    @classmethod #UNTESTED
    def get_proj_stereographic(cls):
        """doc-string"""
        hemi = -90. if cls.truelat1 < 0 else 90.

        lat_ts = (None
                  if wps_info._ismissing(cls.truelat1)
                  else cls.truelat1)

        polarproj = ccrs.Stereographic(central_latitude=hemi,\
                                          central_longitude=cls.stand_lon,\
                                          true_scale_latitude=lat_ts,\
                                          globe=None)
        return polarproj


    @classmethod #UNTESTED
    def get_proj_latlon(cls):
        """doc-string"""
        latlonproj = ccrs.PlateCarree(central_longitude=cls.stand_lon,\
                                            globe=None)
        return latlonproj


    #%% Calculating and Plotting WPS Domain
    @classmethod
    def plot_domain_number(cls, number=1):
        """ Developed by Eric Allen to plot certain domains from namelist.wps file

        Optional Input - Default to coarsest domain, but configurable to the max_domain
        """
        if cls.max_dom < number:
            print("Domain Number Out Of Range: Setting to coarsest domain")
            cls._plot_domain_number = 0
        else:
            cls._plot_domain_number = number - 1

        cls._maxdom_LIMIT = 10
        if cls.max_dom > cls._maxdom_LIMIT:
            print("USING TOO MANY DOMAINS. Re-evaluate your model setup ")
            cls._plot_domain_number = 0

        return cls._plot_domain_number

    @classmethod
    def calc_wps_domain_info(cls):
        """ MODIFIED FROM Xiaodong Chen @lucas-uw on github"""
        # get WPS projection info
        if cls.map_proj == 'lambert':
            cls.wpsproj = wps_info.get_proj_lambert_conformal()
        elif cls.map_proj == "mercator":
            cls.wpsproj = wps_info.get_proj_mercator()
        elif cls.map_proj == "polar":
            cls.wpsproj = wps_info.get_proj_stereographic()
        elif cls.map_proj == "lat-lon":
            cls.wpsproj = wps_info.get_proj_latlon()
        else:
            print("WARNING: PROJECTION UNKNOWN")
            raise ValueError


        # Geodetic, for lat/lon projection
        cls.latlonproj = ccrs.Geodetic()
        #cls.latlonproj = ccrs.PlateCarree()

        length_x = np.zeros((cls.max_dom, 1))
        length_y = np.zeros((cls.max_dom, 1))

        for i in np.arange(0, cls.max_dom):
            length_x[i] = cls.dx[i]*cls.e_we[i] #Units = Meters
            length_y[i] = cls.dy[i]*cls.e_sn[i] #Units = Meters


        center_lat_full = np.zeros((cls.max_dom, 1))
        center_lon_full = np.zeros((cls.max_dom, 1))
        corner_lat_full = np.zeros((cls.max_dom, 4)) # ll, lr, uf, ur
        corner_lon_full = np.zeros((cls.max_dom, 4)) # ll, lr, uf, ur

        center_lat_full[0] = cls.ref_lat
        center_lon_full[0] = cls.ref_lon

        corner_lon_full[0, 0], corner_lat_full[0, 0] =\
            wps_info.calc_corner_point_latlon(float(center_lat_full[0]), float(center_lon_full[0]),\
                                                               cls.e_we[0], cls.e_sn[0],\
                                                               float(cls.dx[0]), float(cls.dy[0]),\
                                                               cls.wpsproj, cls.latlonproj, 'll')
        corner_lon_full[0, 1], corner_lat_full[0, 1] =\
            wps_info.calc_corner_point_latlon(center_lat_full[0], center_lon_full[0],\
                                                   cls.e_we[0], cls.e_sn[0],\
                                                   cls.dx[0], cls.dy[0],\
                                                   cls.wpsproj, cls.latlonproj, 'lr')
        corner_lon_full[0, 2], corner_lat_full[0, 2] =\
            wps_info.calc_corner_point_latlon(center_lat_full[0], center_lon_full[0],\
                                                   cls.e_we[0], cls.e_sn[0],\
                                                   cls.dx[0], cls.dy[0],\
                                                   cls.wpsproj, cls.latlonproj, 'ul')
        corner_lon_full[0, 3], corner_lat_full[0, 3] =\
            wps_info.calc_corner_point_latlon(center_lat_full[0], center_lon_full[0],\
                                                   cls.e_we[0], cls.e_sn[0],\
                                                   cls.dx[0], cls.dy[0],\
                                                   cls.wpsproj, cls.latlonproj, 'ur')
        # Getting information for subsequent domains
        if cls.max_dom > 1:
            for i in np.arange(1, cls.max_dom):
                center_lon_full[i], center_lat_full[i] =\
                    wps_info.calc_center_point_latlon(corner_lat_full[i-1, 0],\
                                                      corner_lon_full[i-1, 0],\
                                                   cls.dx[i-1], cls.dy[i-1],\
                                                   cls.e_we[i], cls.e_sn[i],\
                                                   cls.dx[i], cls.dy[i],\
                                                   cls.i_parent_start[i], cls.j_parent_start[i],\
                                                   cls.wpsproj, cls.latlonproj)


                corner_lon_full[i, 0], corner_lat_full[i, 0] =\
                    wps_info.calc_corner_point_latlon(center_lat_full[i], center_lon_full[i],\
                                                           cls.e_we[i], cls.e_sn[i],\
                                                           cls.dx[i], cls.dy[i],\
                                                           cls.wpsproj, cls.latlonproj, 'll')
                corner_lon_full[i, 1], corner_lat_full[i, 1] =\
                    wps_info.calc_corner_point_latlon(center_lat_full[i], center_lon_full[i],\
                                                           cls.e_we[i], cls.e_sn[i],\
                                                           cls.dx[i], cls.dy[i],\
                                                           cls.wpsproj, cls.latlonproj, 'lr')
                corner_lon_full[i, 2], corner_lat_full[i, 2] =\
                    wps_info.calc_corner_point_latlon(center_lat_full[i], center_lon_full[i],\
                                                           cls.e_we[i], cls.e_sn[i],\
                                                           cls.dx[i], cls.dy[i],\
                                                           cls.wpsproj, cls.latlonproj, 'ul')
                corner_lon_full[i, 3], corner_lat_full[i, 3] =\
                    wps_info.calc_corner_point_latlon(center_lat_full[i], center_lon_full[i],\
                                                           cls.e_we[i], cls.e_sn[i],\
                                                           cls.dx[i], cls.dy[i],\
                                                           cls.wpsproj, cls.latlonproj, 'ur')

        #DO I REALLY WANT THESE VARIABLES PUBLIC, Is there a use that this class doesn't cover for?
        #Making private since you get these variables in the return
        cls._corner_lat_full = corner_lat_full
        cls._corner_lon_full = corner_lon_full
        cls._length_x = length_x
        cls._length_y = length_y

        return cls.wpsproj, cls.latlonproj, cls._corner_lat_full,\
                cls._corner_lon_full, cls._length_x, cls._length_y




    #%%  FROM Xiaodong Chen for finding accurate map information on inner domains to plot
    @staticmethod
    def calc_corner_point_latlon(center_lat, center_lon, e_we, e_ns, dx_h,\
                                     dy_v, wpsproj, latlonproj, loc):
        """ FROM Xiaodong Chen @lucas-uw on github"""

        center_x, center_y = wpsproj.transform_point(center_lon, center_lat, latlonproj)
        if loc == 'll':
            xpt = center_x - dx_h*e_we/2.0
            ypt = center_y - dy_v*e_ns/2.0
        elif loc == 'lr':
            xpt = center_x + dx_h*e_we/2.0
            ypt = center_y - dy_v*e_ns/2.0
        elif loc == 'ul':
            xpt = center_x - dx_h*e_we/2.0
            ypt = center_y + dy_v*e_ns/2.0
        elif loc == 'ur':
            xpt = center_x + dx_h*e_we/2.0
            ypt = center_y + dy_v*e_ns/2.0
        else:
            print("Invalid corner location (VALID: ll, lr, ul, ur)")
            raise ValueError

        corner_lon, corner_lat = latlonproj.transform_point(xpt, ypt, wpsproj)
        # TRANSFORMS TO DD Lat, DD LON

        return corner_lon, corner_lat

    @staticmethod
    def calc_center_point_latlon(corner_lat_parent, corner_lon_parent, dx_parent, dy_parent,\
                                 e_we, e_ns, dx1, dy1, i, j, wpsproj, latlonproj):
        """ FROM Xiaodong Chen @lucas-uw on github"""
        corner_x_parent, corner_y_parent = wpsproj.transform_point(corner_lon_parent,\
                                                        corner_lat_parent, latlonproj)
        center_x_child = corner_x_parent + dx_parent*i + dx1*e_we/2.0
        center_y_child = corner_y_parent + dy_parent*j + dy1*e_ns/2.0
        center_lon_child, center_lat_child = latlonproj.transform_point(center_x_child,\
                                        center_y_child, wpsproj) # FROM #### to Lat/Lon
        return center_lon_child, center_lat_child

    @staticmethod
    def reproject_corners(corner_lons, corner_lats, wpsproj, latlonproj):
        """ FROM Xiaodong Chen @lucas-uw on github"""
        corner_x = np.zeros((4, 1))
        corner_y = np.zeros((4, 1))
        corner_x[0], corner_y[0] = wpsproj.transform_point(corner_lons[0],\
                                                    corner_lats[0], latlonproj)
        corner_x[1], corner_y[1] = wpsproj.transform_point(corner_lons[1],\
                                                    corner_lats[1], latlonproj)
        corner_x[2], corner_y[2] = wpsproj.transform_point(corner_lons[2],\
                                                    corner_lats[2], latlonproj)
        corner_x[3], corner_y[3] = wpsproj.transform_point(corner_lons[3],\
                                                    corner_lats[3], latlonproj)
        return corner_x, corner_y


#%%   PLOT ZOOM(s) and PLOT_ALL_DOMAINS
    @staticmethod
    def plot_zoom_out(corner_x, corner_y, ns_zoom, we_zoom):
        """NOT FULLY TESTED.... FOR ALL PROJECTION TYPES... trial and error"""
        if 0 <= ns_zoom < 40 and 0 <= we_zoom < 40:
            corner_x1 = corner_x + corner_x * we_zoom
            corner_y1 = corner_y + corner_y * ns_zoom
            return corner_x1, corner_y1
        print("Invalid Zoom Level, No Zoom Applied")
        return corner_x, corner_y

    @staticmethod
    def plot_zoom_in(corner_x, corner_y, ns_zoom, we_zoom):
        """HAVE HAD NO ISSUES BUT NOT FULLY TESTED FOR ALL PROJECTION TYPES... trial and error"""
        if 0 <= ns_zoom < 100 and 0 <= we_zoom < 100:
            corner_x1 = corner_x - corner_x * we_zoom/100.
            corner_y1 = corner_y - corner_y * ns_zoom/100.
            return corner_x1, corner_y1
        print("Invalid Zoom Level, No Zoom Applied")
        return corner_x, corner_y

    @staticmethod
    def plot_all_domains(ax1, fig1, max_dom, wpsproj, latlonproj, corner_lat_full,\
                         corner_lon_full, length_x, length_y):
        """
        Plot all the Domains

        Location of dNum string to be dynamic relative to domain size
        and figure size would be much better.

        """
        # d01
        font0 = FontProperties()
        font0.set_weight('bold')
        #print("Number of Domains: ", max_dom)
        colors = ['blue', "white", "red", "cyan", "magenta", "gold",\
                  "black", "green", "yellow", 'pink']
        dnum = ['D01', 'D02', 'D03', 'D04', 'D05', 'D06', 'D07', 'D08', 'D09', 'D10']
        ## ORIGINAL VALUES
        #xbuff = [0.05, 0.05, 0.1]
        #ybuff = [0.9, 1.1, 0.8]
        for i in range(max_dom):
            if i == 0:
                corner_x, corner_y = wps_info.reproject_corners(corner_lon_full[i, :],\
                                            corner_lat_full[i, :], wpsproj, latlonproj)
                ## In Geodetic Coordinates # Want Lat/Lon so I can add gridlines and labels on axis
                ax1.set_xlim([corner_x[0]-length_x[0]/15, corner_x[3]+length_x[0]/15])
                ax1.set_ylim([corner_y[0]-length_y[0]/15, corner_y[3]+length_y[0]/15])
            elif i < 10:
                corner_x, corner_y = wps_info.reproject_corners(corner_lon_full[i, :],\
                                            corner_lat_full[i, :], wpsproj, latlonproj)
            else:
                print("Maximum Domains To Plot is 9")
                raise IndexError

            ax1.add_patch(mpl.patches.Rectangle((corner_x[0], corner_y[0]), length_x[i],\
                                length_y[i], fill=None, lw=3, edgecolor=colors[i], zorder=10))

            ax1.text(corner_x[0]+length_x[i]*0.05, corner_y[0]+length_y[i]*0.85, dnum[i],
                     fontproperties=font0, size=15, color=colors[i], zorder=10) ## Can be improved

        fig1.canvas.draw()

        ax1.add_feature(LAND, edgecolor='k', facecolor='limegreen')
        ax1.add_feature(OCEAN, edgecolor='k', facecolor='deepskyblue')#, facecolor='deepskyblue')
        ax1.add_feature(LAKES, edgecolor='k', facecolor='deepskyblue')
        ax1.coastlines('10m', 'black')
        ax1.add_feature(cartopy.feature.STATES.with_scale('10m'))

        # I DO NOT LIKE THIS LAMBERT PROJ WORK AROUND.
        #   WAIT UNTIL CARTOPY SUPPORTS IT OR DON'T USE OR IMPLEMENT
        #if "Lambert" in str(wpsproj):
            # Add tick marks for lambert projection (only)
            #wps_info.set_lambert_ticks(ax1)
        #    wps_info.set_lambert_ticks(ax1,xskip=5.,yskip=5.,x_thickness=14,y_thickness=14)
        if "Mercator" in str(wpsproj) or "PlateCarree" in str(wpsproj):
            ax1.gridlines(color='lightgrey', linestyle='-', draw_labels=True)
        elif "Polar" in str(wpsproj) or "Stereo" in str(wpsproj) or "lat-lon" in str(wpsproj):
            pass
        else:
            pass

        # Don't think this is needed
        return ax1, fig1

    @staticmethod
    def _ismissing(val, islat=True):
        """
        From WRF-PYTHON....
        Return True if a value is None or out of bounds.
        This function is used to check for invalid latitude/longitude values.
        Args
            val (numeric): A numeric value.
            islat (:obj:`bool`): Set to False if checking for longitude values.

        Returns
            :obj:`bool`: True if the value is None, or an out of bounds value.

        """
        if islat:
            if val is None:
                return True

            if math.fabs(val) > 90.:
                return True
        else:
            if val is None:
                return True

            if math.fabs(val) > 360.:
                return True

        return False
#%%   Xiaodong Chen - Projection Features For Lambert Confromal Plots

    # all these functions are necessary only when LCC projection is used.
    @staticmethod
    def find_side(lstring, side):
        """
        FROM Xiaodong Chen @lucas-uw on github
        Given a shapely LineString which is assumed to be rectangular, return the
        line corresponding to a given side of the rectangle.
        """
        minx, miny, maxx, maxy = lstring.bounds
        points = {'left': [(minx, miny), (minx, maxy)],\
                   'right': [(maxx, miny), (maxx, maxy)],\
                   'bottom': [(minx, miny), (maxx, miny)],\
                   'top': [(minx, maxy), (maxx, maxy)],}
        return sgeom.LineString(points[side])


    @staticmethod
    def lambert_xticks(axs, ticks, size):
        """FROM Xiaodong Chen @lucas-uw on github
             - Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
        te_lambda = lambda xy: xy[0]
        lc_lambda = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
        xticks, xticklabels = wps_info._lambert_ticks(axs, ticks, 'bottom', lc_lambda, te_lambda)
        axs.xaxis.tick_bottom()
        axs.set_xticks(xticks)
        axs.set_xticklabels([axs.xaxis.get_major_formatter()(xtick)\
                            for xtick in xticklabels], size=size)


    @staticmethod
    def lambert_yticks_left(axs, ticks, size):
        """FROM Xiaodong Chen @lucas-uw on github
            - Draw tricks on the left y-axis of a Lamber Conformal projection."""
        te_lambda = lambda xy: xy[1]
        lc_lambda = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
        yticks, yticklabels = wps_info._lambert_ticks(axs, ticks, 'left', lc_lambda, te_lambda)
        axs.yaxis.tick_left()
        axs.set_yticks(yticks)
        axs.set_yticklabels([axs.yaxis.get_major_formatter()(ytick)\
                            for ytick in yticklabels], size=size)

    @staticmethod
    def lambert_yticks_right(axs, ticks, size):
        """FROM Xiaodong Chen @lucas-uw on github
            - Draw ricks on the left y-axis of a Lamber Conformal projection."""
        te_lambda = lambda xy: xy[1]
        lc_lambda = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
        yticks, yticklabels = wps_info._lambert_ticks(axs, ticks, 'right', lc_lambda, te_lambda)
        axs.yaxis.tick_right()
        axs.set_yticks(yticks)
        axs.set_yticklabels([axs.yaxis.get_major_formatter()(ytick)\
                            for ytick in yticklabels], size=size)

    @staticmethod
    def _lambert_ticks(axs, ticks, tick_location, line_constructor, tick_extractor):
        """FROM Xiaodong Chen @lucas-uw on github
            - Get the tick locations and labels for an axis of a Lambert Conformal projection."""
        outline_patch = sgeom.LineString(axs.outline_patch.get_path().vertices.tolist())
        axis = wps_info.find_side(outline_patch, tick_location)
        n_steps = 30
        extent = axs.get_extent(ccrs.PlateCarree())
        _ticks = []
        for tix in ticks:
            xy_line = line_constructor(tix, n_steps, extent)
            proj_xyz = axs.projection.transform_points(ccrs.Geodetic(), xy_line[:, 0],\
                                                       xy_line[:, 1])
            xyt = proj_xyz[..., :2]
            lstr = sgeom.LineString(xyt.tolist())
            locs = axis.intersection(lstr)
            if not locs:
                tick = [None]
            else:
                tick = tick_extractor(locs.xy_line)
            _ticks.append(tick[0])
        # Remove ticks that aren't visible:
        ticklabels = copy(ticks)
        #ticklabels.tolist()
        while True:
            ##Eric: HAD TO MAKE A CHANGE HERE BECAUSE OFTEN
                # ticklabels is a numpy array and not a list.
            ##np.arrays cannot pop a value like that
            try:
                index = _ticks.index(None)
            except ValueError:
                break
            #print (type(_ticks),type(ticklabels))
            _ticks.pop(index)
            try:
                list(ticklabels).pop(index)
            except:
                ticklabels.pop(index)
        return _ticks, ticklabels

    #!! PROBLEM
    @staticmethod
    def set_lambert_ticks(axs, xskip=10., yskip=10., x_thickness=14, y_thickness=14):
        """
        PROBLEMS HERE.... extent_rounded doesn't match the real extent but if the
        extent is less than a skip increment then it becomes all messed up.
        labels are inconsistent. Not sure why it plots sometimes and not othertimes.
        I want to be able to have gridlines on set increments like at 5N,10N,15N,20N
        or 10N,20N,30N.
        Gridlines sometimes don't cover the entire area.
        """
        extent = list(axs.get_extent(ccrs.PlateCarree()))

        # Incase the extent is less than the x/y skip value.
            # But then that messes with the only working gridlines
        extent_rounded = [round(val/10.0) * 10 for val in extent]
        #print("Rounded",extent_rounded)
        if extent_rounded[0] == extent_rounded[1]:
            extent_rounded[0] -= 10
            extent_rounded[0] += 10
        if extent_rounded[0] == extent_rounded[1]:
            extent_rounded[0] -= 10
            extent_rounded[0] += 10

        #print("Extent",extent)
        #print("Rounded",extent_rounded)

        xticks = np.arange(extent[0], extent[1]+1, xskip) # added one to get the full extent
        yticks = np.arange(extent[2], extent[3]+1, yskip) # added one to get the full extent
        #xticks = np.arange(extent_rounded[0],extent_rounded[1]+1,xskip)
        #yticks = np.arange(extent_rounded[2],extent_rounded[3]+1,yskip)


        #print("x_ticks ",xticks,"\t\t","y_ticks ", yticks)
        axs.gridlines(xlocs=xticks, ylocs=yticks)
        #axs.gridlines(color='lightgrey', linestyle='-', draw_labels=True)

        # OPTION 1: Use Cartopy's LAT/LON FORMATTER
        #axs.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
        #axs.yaxis.set_major_formatter(LATITUDE_FORMATTER)

        # OPTION 2: SET PRECISION
        axs.xaxis.set_major_formatter(ticker.FormatStrFormatter("%3.2f"))
        axs.yaxis.set_major_formatter(ticker.FormatStrFormatter("%3.2f"))

        wps_info.lambert_xticks(axs, xticks, x_thickness)
        wps_info.lambert_yticks_left(axs, yticks, y_thickness)
        return axs
