""" This script is for building figures for papers. We wan't nice figures and this will allow for merging figures in the way that I want. """ import os import numpy as np import scipy.stats as sp import pandas as pd import matplotlib.pyplot as plt import matplotlib import matplotlib.animation as animation from matplotlib.animation import FuncAnimation from mpl_toolkits.axes_grid1 import make_axes_locatable import datetime as dt from pyhdf.HDF import * from pyhdf.VS import * import netCDF4 as cdf import cartopy.crs as ccrs import cartopy.feature as feature from warnings import filterwarnings import jja_global_avg import requests from time import sleep import merra_first import bam_analysis import bam_analysis_seasonal import bam_merra_results import go_stats import mast_plot import oned_corr import indexes import cloud_controls import coarsify import import_eke_data import so_story_2 import so_story import get_ceres import area_calc import modeled_lwcld def dt_centered(reg_slp1, reg_slp2, season1, season2, pval1, pval2): cmap_send = matplotlib.cm.get_cmap('RdBu') shifted_cmap1 = mast_plot.shiftedColorMap(cmap_send, 1, 0.5, 0, 'shifted') shifted_cmap2 = mast_plot.shiftedColorMap(cmap_send, 1, 0.5, 0, 'shifted') fig, ax = plt.subplots(2,1, figsize=(13,15)) proj = ccrs.PlateCarree(central_longitude=-180) raw = ccrs.PlateCarree() ax1 = plt.subplot(2,1,1, projection=proj) g1 = ax1.gridlines(crs=raw, draw_labels = False, linestyle = 'dotted') ax1.coastlines() p = ax1.imshow(reg_slp1, transform=raw, origin='lower', aspect='equal', cmap=shifted_cmap1, vmin=-0.1, vmax=0.1) cb = plt.colorbar(p, ax = ax1, shrink=0.66, pad=0.03) axc = cb.ax axc.set_yticks([-0.1, -0.075, -0.05, -0.025, 0, 0.025, 0.05, 0.075, 0.1]) axc.set_yticklabels([-0.1, -0.075, -0.05, -0.025, 0, 0.025, 0.05, 0.075, 0.1], size=13) ax1.set_xticks([60., 120., 180., -120., -60.], crs=ccrs.PlateCarree()) ax1.set_xticklabels([60., 120., 180., -120., -60.], size = 13) ax1.set_yticks([80.,60.,40.,20.,0.,-20.,-40.,-60.,-80.], crs=ccrs.PlateCarree()) ax1.set_yticklabels([80.,60.,40.,20.,0.,-20.,-40.,-60.,-80.], size=13) # ax1.set_title('Seasonally averaged ENSO index (lag -1),\nregressed with CESM2 model generated low cloud occurrence, slope '+ str(season1)+ ' (1996-2012)', size=15) ax1.set_xlabel('Longitude', size=13) ax1.set_ylabel('Latitude', size=13) ax1.text(-200,107, 'a)', size=17) ax1.text(225,0, 'occurrence frequency/\n1STD ENSO') tstep = np.arange(0, len(pval1)) for t in tstep: lat = pval1[t].lat lon = pval1[t].lon lon_s = go_stats.reformat_for_cp_shift(lon) ax1.scatter(lon_s+0.75, lat+1.25, color='black', alpha=1, s=1) ax2 = plt.subplot(2,1,2, projection=proj) g2 = ax2.gridlines(crs=raw, draw_labels=False, linestyle = 'dotted') ax2.coastlines() p = ax2.imshow(reg_slp2, transform=raw, origin='lower', aspect='equal', cmap=shifted_cmap2, vmin=-0.1, vmax=0.1) cb = plt.colorbar(p, ax = ax2, shrink=0.66, pad=0.03) axc = cb.ax axc.set_yticks([-0.1, -0.075, -0.05, -0.025, 0, 0.025, 0.05, 0.075, 0.1]) axc.set_yticklabels([-0.1, -0.075, -0.05, -0.025, 0, 0.025, 0.05, 0.075, 0.1], size=13) ax2.set_xticks([60., 120., 180., -120., -60.], crs=ccrs.PlateCarree()) ax2.set_xticklabels([60., 120., 180., -120., -60.], size = 13) ax2.set_yticks([80.,60.,40.,20.,0.,-20.,-40.,-60.,-80.], crs=ccrs.PlateCarree()) ax2.set_yticklabels([80.,60.,40.,20.,0.,-20.,-40.,-60.,-80.], size=13) # ax2.set_title('Seasonally averaged ENSO index (lag -1),\nregressed with CESM2 model generated low cloud occurrence, slope '+ str(season2)+ ' (1996-2012)', size=15) ax2.set_xlabel('Longitude', size=13) ax2.set_ylabel('Latitude', size=13) ax2.text(-200, 107, 'b)', size=17) ax2.text(225,0, 'occurrence frequency/\n1STD ENSO') tstep = np.arange(0, len(pval2)) for t in tstep: lat = pval2[t].lat lon = pval2[t].lon lon_s = go_stats.reformat_for_cp_shift(lon) ax2.scatter(lon_s+0.75, lat+1.25, color='black', alpha=1, s=1) plt.tight_layout(pad=1) plt.savefig('/uufs/chpc.utah.edu/common/home/u1113223/grad_research/plots/pub_plots/enso_cld_occ/gisst_cldocc_nt_jjadjf_stats_lag-1_flip.png') if __name__ == '__main__': years_arr1 = np.arange(1996,2013) # years_arr1 = np.delete(years_arr1, np.argwhere(years_arr1 == 2011)) # years_arr = np.delete(years_arr, np.argwhere(years_arr == 2012)) years1 = go_stats.make_np_list(years_arr1) months_arr1 = np.arange(6,9) # months_arr = np.array([12,1,2]) months1 = go_stats.make_np_list(months_arr1) months1 = go_stats.fix_days(months1) months_inx_arr1 = np.arange(5,8) # months_inx_arr = np.array([11,12,1]) months_inx1 = go_stats.make_np_list(months_inx_arr1) months_inx1 = go_stats.fix_days(months_inx1) season1 = 'jja' years_arr2 = np.arange(1996,2013) # years_arr2 = np.delete(years_arr2, np.argwhere(years_arr2 == 2011)) # years_arr = np.delete(years_arr, np.argwhere(years_arr == 2012)) years2 = go_stats.make_np_list(years_arr2) # months_arr2 = np.arange(9,12) months_arr2 = np.array([12,1,2]) months2 = go_stats.make_np_list(months_arr2) months2 = go_stats.fix_days(months2) # months_inx_arr2 = np.arange(8,11) months_inx_arr2 = np.array([11,12,1]) months_inx2 = go_stats.make_np_list(months_inx_arr2) months_inx2 = go_stats.fix_days(months_inx2) season2 = 'djf' fname1 = season1+'_mast.nc4' fname2 = season2+'_mast_w2011_corr.nc4' # fname2 = season2+'_mast.nc4' enso1, etime1 = indexes.enso_full(months_inx1, years_arr1, 0, 0) enso2, etime2 = indexes.enso_full(months_inx2, years_arr2, 1, 0) print(etime2) #***EIS Stuff*** # eis1, time1, lat, lon = cloud_controls.import_eis_wstuff(season1) # print(time1) # print(etime1) # eis2, time2, lat, lon = cloud_controls.import_eis_wstuff(season2) #***Radiation Stuff*** # rad1, lat, lon, time1 = get_ceres.call_main_toa_rad(years_arr1, months_arr1, 0) # rad2, lat, lon, time2 = get_ceres.call_main_toa_rad(years_arr2, months_arr2, 0) #***SST Stuff*** # lat, lon, time1, sst1 = so_story_2.get_sst_data(fname1, 1) # lat, lon, time2, sst2 = so_story_2.get_sst_data(fname2, 1) # sst_m1 = bam_merra_results.make_seas_mean(sst1, 3) # sst_m2 = bam_merra_results.make_seas_mean(sst2, 3) # sst_c1, lat_c, lon_c = so_story_2.make_ts_coarse(sst_m1, lat, lon, 1, 1) # sst_c2, lat_c, lon_c = so_story_2.make_ts_coarse(sst_m2, lat, lon, 1, 1) # print(sst_c1.shape) #***Cloud Occurrence*** # cld_occ1, lat, lon = oned_corr.eof_call_main_2(years1, months1, 0, 0, 0) # print(cld_occ1.shape) # cld_occ2, lat, lon = oned_corr.eof_call_main_2(years2, months2, 0, 1, 0) # cld_occ1_fl, lon_fl = oned_corr.prep_co_data(cld_occ1, lon) # cld_occ2_fl, lon_fl = oned_corr.prep_co_data(cld_occ2, lon) #***MERRA2 Low cloud occurrence*** # lat, lon, time1, lowcld1, olr1, var_lowcld1 = so_story.get_so_rad_data(fname1, 1) # lat, lon, time2, lowcld2, olr2, var_lowcld2 = so_story.get_so_rad_data(fname2, 1) # # # print(time2[:16]) # lowcld_m1 = bam_merra_results.make_seas_mean(lowcld1, 3) # lowcld_m2 = bam_merra_results.make_seas_mean(lowcld2[48:], 3) # lowcld_m1_c, lat_c, lon_c = coarsify.make_me_coarse_3d(lowcld_m1, lat, lon) # lowcld_m2_c, lat_c, lon_c = coarsify.make_me_coarse_3d(lowcld_m2, lat, lon) # print(lowcld_m2_c.shape) #***CESM2 Low cloud occurrence*** lw_cld_s1, lat_c, lon_c, time = modeled_lwcld.run_main_modeled(years_arr1, months_arr1) lw_cld_s2, lat_c, lon_c, time = modeled_lwcld.run_main_modeled(years_arr2, months_arr2) # area_grid = area_calc.get_area_grid(lat) # #***THis is for the pacific total radiation restriction*** # pacific_bounds = [42.5,-42.5,-35,120] # bounds_arr = jja_global_avg.restrict_domain(lat, lon, pacific_bounds[0], pacific_bounds[1], # pacific_bounds[2], pacific_bounds[3]) # area_grid, lat_r, lon_r = jja_global_avg.region_arr(area_grid, lat, lon, bounds_arr) # print(area_grid[3,5:10]) # tstep = np.arange(0,rad1.shape[0]) # for t in tstep: # data_adj_hold = rad1[t]*area_grid # if t == 0: # data_adj1 = np.array([data_adj_hold]) # continue # data_adj1 = np.append(data_adj1, np.array([data_adj_hold]), axis=0) # data_adj1 = data_adj1/(5.1*10**14) # tstep = np.arange(0,rad2.shape[0]) # for t in tstep: # data_adj_hold = rad2[t]*area_grid # if t == 0: # data_adj2 = np.array([data_adj_hold]) # continue # data_adj2 = np.append(data_adj2, np.array([data_adj_hold]), axis=0) # data_adj2 = data_adj2/(5.1*10**14) corr_arr1, reg_slp1, reg_inc1, pvals1 = bam_analysis.build_cov_arr_1st_dem(lw_cld_s1/100, 'allwstats', enso1, lat_c, lon_c, -10, 10, 0.2, 'noname') corr_arr2, reg_slp2, reg_inc2, pvals2 = bam_analysis.build_cov_arr_1st_dem(lw_cld_s2/100, 'allwstats', enso2, lat_c, lon_c, -10, 10, 0.2, 'noname') lon_h = np.arange(0, 180, 6) lon_h1 = np.arange(-180, 0, 6) lon_res = np.append(lon_h, lon_h1) # p = plt.imshow(reg_inc1) # plt.colorbar(p) # plt.savefig('/uufs/chpc.utah.edu/common/home/u1113223/grad_research/my_code/climate_proj/tester.png') #***Also for low cloud occurrence # reg_slp1, reg_inc1, corr_arr1, pvals1 = go_stats.add_cells(cld_occ1_fl, enso1, lat, lon_fl) # reg_slp2, reg_inc2, corr_arr2, pvals2 = go_stats.add_cells(cld_occ2_fl, enso2, lat, lon_fl) # reg_slp1 = oned_corr.flip_co(reg_slp1) # reg_slp2 = oned_corr.flip_co(reg_slp2) p = plt.imshow(reg_slp2, origin='lower') plt.colorbar(p) plt.savefig('/uufs/chpc.utah.edu/common/home/u1113223/grad_research/my_code/climate_proj/testinger.png') plt.close() #***For Pacific low cloud integration. pacific_bounds = [60,-60,-60,122] bounds_arr = jja_global_avg.restrict_domain(lat_c, lon_res, pacific_bounds[0], pacific_bounds[1], pacific_bounds[2], pacific_bounds[3]) res_arr1, lat_r, lon_r = jja_global_avg.region_arr(reg_slp1, lat_c, lon_res, bounds_arr) res_arr2, lat_r, lon_r = jja_global_avg.region_arr(reg_slp2, lat_c, lon_res, bounds_arr) print(lon_r) p = plt.imshow(res_arr2, origin='lower') plt.colorbar(p) plt.savefig('/uufs/chpc.utah.edu/common/home/u1113223/grad_research/my_code/climate_proj/testing.png') plt.close() print(np.nansum(res_arr1)) print(np.nansum(res_arr2)) reg_slp1 = oned_corr.flip_co(reg_slp1) reg_slp2 = oned_corr.flip_co(reg_slp2) # p= plt.imshow(reg_slp1, origin='lower') # plt.colorbar(p) # plt.savefig('/uufs/chpc.utah.edu/common/home/u1113223/grad_research/my_code/climate_proj/testing.png') # plt.close() # p= plt.imshow(reg_slp2, origin='lower') # plt.colorbar(p) # plt.savefig('/uufs/chpc.utah.edu/common/home/u1113223/grad_research/my_code/climate_proj/tester.png') # plt.close() # print(np.sum(reg_slp1, dtype=np.double)) # print(np.sum(reg_slp1, dtype=np.float64)) # print(np.sum(reg_slp2, dtype=np.float64)) # exit() #***This is for calculating the pacific radiative impact # print(test_arr.shape) # print(lat_r) # print(lon_r) # exit() dt_centered(reg_slp1, reg_slp2, season1, season2, pvals1, pvals2)