##################################################################################### #Code to create gridded plots for paired objects. # #-Run runMODE.py first to create the netCDF files (modify mode_config for MODE sensitivity studies) #-Run createCSV.py to create a csv of paired objects #-Run createCSV_simpleobj.py to create a csv of simple objects #-Run plotting.py to create non-grid plots for paired objects #-Run plotting_grid.py to create grid plots for paired objects #-Run plotting_grid_simpleobj.py to create grid plots for simple objects (needs work) # #Created in the summer/fall of 2020 by CLC, editted by MJE on 20210601. #Additional modifications by AJ throughout summer of 2021. #Further documented by MJE. 20210921. # ###################################################################################### #!/usr/bin/python import pygrib import sys import subprocess import string import gzip import time import os #from netCDF4 import Dataset import numpy as np import datetime import csv import math import haversine import glob import cartopy.crs as ccrs # import projections import cartopy.feature as cf # import features import numpy as np import matplotlib as mpl import matplotlib.pyplot as pyplot mpl.use('Agg') #So plots can be saved in cron import matplotlib.pyplot as plt from windrose import WindroseAxes import matplotlib.cm as cm import cartopy_maps import importlib from scipy.io import netcdf from netCDF4 import Dataset from scipy import ndimage from scipy import stats from scipy import interpolate from mpl_toolkits.axes_grid1 import make_axes_locatable from scipy.ndimage.filters import gaussian_filter #Load predefined functions import METConFigGenV100 import METLoadEnsemble import METPlotEnsemble importlib.reload(METConFigGenV100) importlib.reload(METLoadEnsemble) importlib.reload(METPlotEnsemble) #######################LIST OF VARIABLES######################################################## FIG_DIR = 'hpc@vm-lnx-rzdm01:/home/people/hpc/ftp/erickson/lapenta/matched' working_dir = '/export/hpc-lw-dtbdev5/merickson/code/python/lapenta/work/' validday = 1 ero_category = 'SLGT' #for displacement maps grid_res = 3 #for frequency maps grid_res_freq= 2 ################################################################################################ #Define the output file latlon = working_dir+'latlon_vday'+str(validday)+'_ERO'+ero_category+'.csv' #initial lon_mod=[] lat_mod=[] lon_obs=[] lat_obs=[] date=[] valid=[] magnitude=[] angle=[] area_obs=[] area_mod=[] ydiff_lat=[] xdiff_lon=[] #jja lon_mod_jja=[] lat_mod_jja=[] lon_obs_jja=[] lat_obs_jja=[] date_jja=[] valid_jja=[] magnitude_jja=[] angle_jja=[] area_obs_jja=[] area_mod_jja=[] ydiff_lat_jja=[] xdiff_lon_jja=[] #son lon_mod_son=[] lat_mod_son=[] lon_obs_son=[] lat_obs_son=[] date_son=[] valid_son=[] magnitude_son=[] angle_son=[] area_obs_son=[] area_mod_son=[] ydiff_lat_son=[] xdiff_lon_son=[] #djf lon_mod_djf=[] lat_mod_djf=[] lon_obs_djf=[] lat_obs_djf=[] date_djf=[] valid_djf=[] magnitude_djf=[] angle_djf=[] area_obs_djf=[] area_mod_djf=[] ydiff_lat_djf=[] xdiff_lon_djf=[] #mam lon_mod_mam=[] lat_mod_mam=[] lon_obs_mam=[] lat_obs_mam=[] date_mam=[] valid_mam=[] magnitude_mam=[] angle_mam=[] area_obs_mam=[] area_mod_mam=[] ydiff_lat_mam=[] xdiff_lon_mam=[] target=open(latlon,'r') alllines=target.readlines() for line in alllines: values=line.strip().split(',') date=np.append(date,float(values[0][0:8])) valid=np.append(valid,values[1]) lat_mod=np.append(lat_mod,float(values[2])) lon_mod=np.append(lon_mod,float(values[3])) area_mod=np.append(area_mod,float(values[4])) lat_obs=np.append(lat_obs,float(values[5])) lon_obs=np.append(lon_obs,float(values[6])) area_obs=np.append(area_obs,float(values[7])) magnitude=np.append(magnitude,float(values[8])) angle=np.append(angle,float(values[9])) ydiff_lat=np.append(ydiff_lat,float(values[10])) xdiff_lon=np.append(xdiff_lon,float(values[11])) #June,July,August if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: #print('jja') #print(values[0][0:8]) date_jja=np.append(date_jja,float(values[0][0:8])) valid_jja=np.append(valid_jja,values[1]) lat_mod_jja=np.append(lat_mod_jja,float(values[2])) lon_mod_jja=np.append(lon_mod_jja,float(values[3])) area_mod_jja=np.append(area_mod_jja,float(values[4])) lat_obs_jja=np.append(lat_obs_jja,float(values[5])) lon_obs_jja=np.append(lon_obs_jja,float(values[6])) area_obs_jja=np.append(area_obs_jja,float(values[7])) magnitude_jja=np.append(magnitude_jja,float(values[8])) angle_jja=np.append(angle_jja,float(values[9])) ydiff_lat_jja=np.append(ydiff_lat_jja,float(values[10])) xdiff_lon_jja=np.append(xdiff_lon_jja,float(values[11])) #September,October,November elif np.float(values[0][4:6])==9 or np.float(values[0][4:6])==10 or np.float(values[0][4:6])==11: #print('son') #print(values[0][0:8]) date_son=np.append(date_son,float(values[0][0:8])) valid_son=np.append(valid_son,values[1]) lat_mod_son=np.append(lat_mod_son,float(values[2])) lon_mod_son=np.append(lon_mod_son,float(values[3])) area_mod_son=np.append(area_mod_son,float(values[4])) lat_obs_son=np.append(lat_obs_son,float(values[5])) lon_obs_son=np.append(lon_obs_son,float(values[6])) area_obs_son=np.append(area_obs_son,float(values[7])) magnitude_son=np.append(magnitude_son,float(values[8])) angle_son=np.append(angle_son,float(values[9])) ydiff_lat_son=np.append(ydiff_lat_son,float(values[10])) xdiff_lon_son=np.append(xdiff_lon_son,float(values[11])) #December,Jan,Feb elif np.float(values[0][4:6])==12 or np.float(values[0][4:6])==1 or np.float(values[0][4:6])==2: #print('djf') #print(values[0][0:8]) date_djf=np.append(date_djf,float(values[0][0:8])) valid_djf=np.append(valid_djf,values[1]) lat_mod_djf=np.append(lat_mod_djf,float(values[2])) lon_mod_djf=np.append(lon_mod_djf,float(values[3])) area_mod_djf=np.append(area_mod_djf,float(values[4])) lat_obs_djf=np.append(lat_obs_djf,float(values[5])) lon_obs_djf=np.append(lon_obs_djf,float(values[6])) area_obs_djf=np.append(area_obs_djf,float(values[7])) magnitude_djf=np.append(magnitude_djf,float(values[8])) angle_djf=np.append(angle_djf,float(values[9])) ydiff_lat_djf=np.append(ydiff_lat_djf,float(values[10])) xdiff_lon_djf=np.append(xdiff_lon_djf,float(values[11])) #March,April,May elif np.float(values[0][4:6])==3 or np.float(values[0][4:6])==4 or np.float(values[0][4:6])==5: #print('mam') #print(values[0][0:8]) date_mam=np.append(date_mam,float(values[0][0:8])) valid_mam=np.append(valid_mam,values[1]) lat_mod_mam=np.append(lat_mod_mam,float(values[2])) lon_mod_mam=np.append(lon_mod_mam,float(values[3])) area_mod_mam=np.append(area_mod_mam,float(values[4])) lat_obs_mam=np.append(lat_obs_mam,float(values[5])) lon_obs_mam=np.append(lon_obs_mam,float(values[6])) area_obs_mam=np.append(area_obs_mam,float(values[7])) magnitude_mam=np.append(magnitude_mam,float(values[8])) angle_mam=np.append(angle_mam,float(values[9])) ydiff_lat_mam=np.append(ydiff_lat_mam,float(values[10])) xdiff_lon_mam=np.append(xdiff_lon_mam,float(values[11])) if ero_category=='SLGT': bin = [0,50,100,150,200,250,300,350,400,450,500] elif ero_category=='MDT': bin = [0,25,50,75,100,125,150,175,200,250,300] elif ero_category=='HIGH': bin = [0,20,40,60,80,100,120,140,160,170,180] ########################################make scatter plots######################################################################################################## if len(angle) > 0: #scatterplot for all seasons latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) plt.scatter(lon_obs,lat_obs, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_matched_year_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_matched_year_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) if len(angle_jja) > 0: #scatterplot jja latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) plt.scatter(lon_obs_jja,lat_obs_jja, color='blue',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('Summer Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_matched_jja_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_matched_jja_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) if len(angle_son) > 0: #scatterplot son latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) plt.scatter(lon_obs_son,lat_obs_son, color='orange',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('Fall Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_matched_son_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_matched_son_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) if len(angle_djf) > 0: #scatterplot djf latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) plt.scatter(lon_obs_djf,lat_obs_djf, color='red',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('Winter Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_matched_djf_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_matched_djf_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) if len(angle_mam) > 0: #scatterplot mam latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) plt.scatter(lon_obs_mam,lat_obs_mam, color='green',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('Spring Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_matched_mam_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_matched_mam_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #################################calculate and plot frequency map################################################################################################# if ero_category == 'SLGT': vmax_all=30 vmax_season=20 if ero_category =='MDT': vmax_all=15 vmax_season=5 if ero_category=='HIGH': vmax_all=4 vmax_season=2 #all seasons [lon_nat,lat_nat] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res_freq)+1), np.linspace(50,22,int((50-22)/grid_res_freq)+1)) grid_count_obs = np.zeros(lat_nat.shape) grid_count_mod = np.zeros(lat_nat.shape) if len(angle) > 0: for x in range(0,lon_nat.shape[1]-1): for y in range(0,lat_nat.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) grid_count_obs[y,x] = np.nansum((lat_obs<=lat_nat[y,x]) & (lat_obs>lat_nat[y+1,x]) & (lon_obs>lon_nat[y,x]) & (lon_obs<=lon_nat[y,x+1])) grid_count_mod[y,x] = np.nansum((lat_mod<=lat_nat[y,x]) & (lat_mod>lat_nat[y+1,x]) & (lon_mod>lon_nat[y,x]) & (lon_mod<=lon_nat[y,x+1])) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs=ax.pcolormesh(lon_nat, lat_nat, grid_count_obs, cmap = 'Reds',vmin=0, vmax=vmax_all, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Day '+str(validday)+' Analysis Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_matched_year_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_year_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs=ax.pcolormesh(lon_nat, lat_nat, grid_count_mod, cmap = 'Reds',vmin=0, vmax=vmax_all, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Day '+str(validday)+' ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_matched_year_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_year_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #jja [lon_nat_jja,lat_nat_jja] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res_freq)+1), np.linspace(50,22,int((50-22)/grid_res_freq)+1)) grid_count_obs_jja = np.zeros(lat_nat_jja.shape) grid_count_mod_jja = np.zeros(lat_nat_jja.shape) if len(angle_jja) > 0: for x in range(0,lon_nat_jja.shape[1]-1): for y in range(0,lat_nat_jja.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) grid_count_obs_jja[y,x] = np.nansum((lat_obs_jja<=lat_nat_jja[y,x]) & (lat_obs_jja>lat_nat_jja[y+1,x]) & (lon_obs_jja>lon_nat_jja[y,x]) &\ (lon_obs_jja<=lon_nat_jja[y,x+1])) grid_count_mod_jja[y,x] = np.nansum((lat_mod_jja<=lat_nat_jja[y,x]) & (lat_mod_jja>lat_nat_jja[y+1,x]) & (lon_mod_jja>lon_nat_jja[y,x]) &\ (lon_mod_jja<=lon_nat_jja[y,x+1])) #print(np.unique(grid_count_jja)) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs_jja=ax.pcolormesh(lon_nat_jja, lat_nat_jja, grid_count_obs_jja, cmap = 'Reds',vmin=0,vmax=vmax_season, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Summer Day '+str(validday)+' Analysis Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs_jja,ax=ax) plt.savefig(working_dir+'centroid_freq_matched_jja_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_jja_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs_jja=ax.pcolormesh(lon_nat_jja, lat_nat_jja, grid_count_mod_jja, cmap = 'Reds',vmin=0,vmax=vmax_season, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Summer Day '+str(validday)+' ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs_jja,ax=ax) plt.savefig(working_dir+'centroid_freq_matched_jja_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_jja_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #son [lon_nat_son,lat_nat_son] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res_freq)+1), np.linspace(50,22,int((50-22)/grid_res_freq)+1)) grid_count_obs_son = np.zeros(lat_nat_son.shape) grid_count_mod_son = np.zeros(lat_nat_son.shape) if len(angle_son) > 0: for x in range(0,lon_nat_son.shape[1]-1): for y in range(0,lat_nat_son.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) grid_count_obs_son[y,x] = np.nansum((lat_obs_son<=lat_nat_son[y,x]) & (lat_obs_son>lat_nat_son[y+1,x]) & (lon_obs_son>lon_nat_son[y,x])\ & (lon_obs_son<=lon_nat_son[y,x+1])) grid_count_mod_son[y,x] = np.nansum((lat_mod_son<=lat_nat_son[y,x]) & (lat_mod_son>lat_nat_son[y+1,x]) & (lon_mod_son>lon_nat_son[y,x])\ & (lon_mod_son<=lon_nat_son[y,x+1])) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs_son=ax.pcolormesh(lon_nat_son, lat_nat_son, grid_count_obs_son, cmap = 'Reds',vmin=0,vmax=vmax_season, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Fall Day '+str(validday)+' Analysis Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs_son,ax=ax) plt.savefig(working_dir+'centroid_freq_matched_son_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_son_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs_son=ax.pcolormesh(lon_nat_son, lat_nat_son, grid_count_mod_son, cmap = 'Reds',vmin=0,vmax=vmax_season, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Fall Day '+str(validday)+' ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs_son,ax=ax) plt.savefig(working_dir+'centroid_freq_matched_son_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_son_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #djf [lon_nat_djf,lat_nat_djf] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res_freq)+1), np.linspace(50,22,int((50-22)/grid_res_freq)+1)) grid_count_obs_djf = np.zeros(lat_nat_djf.shape) grid_count_mod_djf = np.zeros(lat_nat_djf.shape) if len(angle_djf) > 0: for x in range(0,lon_nat_djf.shape[1]-1): for y in range(0,lat_nat_djf.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) grid_count_obs_djf[y,x] = np.nansum((lat_obs_djf<=lat_nat_djf[y,x]) & (lat_obs_djf>lat_nat_djf[y+1,x]) & (lon_obs_djf>lon_nat_djf[y,x]) &\ (lon_obs_djf<=lon_nat_djf[y,x+1])) grid_count_mod_djf[y,x] = np.nansum((lat_mod_djf<=lat_nat_djf[y,x]) & (lat_mod_djf>lat_nat_djf[y+1,x]) & (lon_mod_djf>lon_nat_djf[y,x]) &\ (lon_mod_djf<=lon_nat_djf[y,x+1])) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs_djf=ax.pcolormesh(lon_nat_djf, lat_nat_djf, grid_count_obs_djf, cmap = 'Reds',vmin=0,vmax=vmax_season, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Winter Day '+str(validday)+' Analysis Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs_djf,ax=ax) plt.savefig(working_dir+'centroid_freq_matched_djf_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_djf_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs_djf=ax.pcolormesh(lon_nat_djf, lat_nat_djf, grid_count_mod_djf, cmap = 'Reds',vmin=0,vmax=vmax_season, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Winter Day '+str(validday)+' ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs_djf,ax=ax) plt.savefig(working_dir+'centroid_freq_matched_djf_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_djf_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #mam [lon_nat_mam,lat_nat_mam] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res_freq)+1), np.linspace(50,22,int((50-22)/grid_res_freq)+1)) grid_count_obs_mam = np.zeros(lat_nat_mam.shape) grid_count_mod_mam = np.zeros(lat_nat_mam.shape) if len(angle_mam) > 0: for x in range(0,lon_nat_mam.shape[1]-1): for y in range(0,lat_nat_mam.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) grid_count_obs_mam[y,x] = np.nansum((lat_obs_mam<=lat_nat_mam[y,x]) & (lat_obs_mam>lat_nat_mam[y+1,x]) & (lon_obs_mam>lon_nat_mam[y,x]) &\ (lon_obs_mam<=lon_nat_mam[y,x+1])) grid_count_mod_mam[y,x] = np.nansum((lat_mod_mam<=lat_nat_mam[y,x]) & (lat_mod_mam>lat_nat_mam[y+1,x]) & (lon_mod_mam>lon_nat_mam[y,x]) &\ (lon_mod_mam<=lon_nat_mam[y,x+1])) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs_mam=ax.pcolormesh(lon_nat_mam, lat_nat_mam, grid_count_obs_mam, cmap = 'Reds',vmin=0,vmax=vmax_season, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Spring Day '+str(validday)+' Analysis Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs_mam,ax=ax) plt.savefig(working_dir+'centroid_freq_matched_mam_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_mam_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) cs_mam=ax.pcolormesh(lon_nat_mam, lat_nat_mam, grid_count_mod_mam, cmap = 'Reds',vmin=0,vmax=vmax_season, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Spring Day '+str(validday)+' ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs_mam,ax=ax) plt.savefig(working_dir+'centroid_freq_matched_mam_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_matched_mam_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) ############################################magnitude/angle average calculations and plotting##################################################################### if ero_category == 'SLGT': num_points = 5 if ero_category == 'MDT': num_points = 3 if ero_category == 'HIGH': num_points = 1 #all seasons [lon_nat,lat_nat] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res)+1), np.linspace(50,22,int((50-22)/grid_res)+1)) xdiff_lon_ave = np.zeros(lat_nat.shape) ydiff_lat_ave = np.zeros(lat_nat.shape) grid_count_me = np.zeros(lat_nat.shape) #example #cows = np.array([2,4,5]) #log_cows = cows > 2 #cows2 = cows[log_cows] #print(log_cows) #print(cows2) #cow if len(angle) > 0: for x in range(0,lon_nat.shape[1]-1): for y in range(0,lat_nat.shape[0]-1): grid_count_me[y,x] =np.nansum((lat_obs<=lat_nat[y,x]) & (lat_obs>lat_nat[y+1,x]) & (lon_obs>lon_nat[y,x]) & (lon_obs<=lon_nat[y,x+1])) if grid_count_me[y,x] >= num_points: log_array = (lat_obs<=lat_nat[y,x]) & (lat_obs>lat_nat[y+1,x]) & (lon_obs>lon_nat[y,x]) & (lon_obs<=lon_nat[y,x+1]) xdiff_lon_ave[y,x] = np.nanmedian(xdiff_lon[log_array]) ydiff_lat_ave[y,x] = np.nanmedian(ydiff_lat[log_array]) else: pass #finds the magnitude #print(np.sqrt(xdiff_lon_ave[2,1]**2+ydiff_lat_ave[2,1]**2)) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) ma=plt.quiver(lon_nat+(grid_res/2), lat_nat-(grid_res/2), xdiff_lon_ave, ydiff_lat_ave, transform=ccrs.PlateCarree(), zorder = 1, units='inches',angles='xy',\ scale=4,scale_units='inches',width=0.02, linewidths=0.2) plt.title('Day '+str(validday)+' Paired Object Displacement for '+ero_category,fontsize=14) cs = plt.quiverkey(ma, 0.8, 0.28, 2.0,'2 Degrees', labelpos='E',coordinates='figure',transform=ccrs.PlateCarree(), color= 'Red') plt.savefig(working_dir+'arrow_map_displ_year_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/arrow_map_displ_year_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #jja [lon_nat_jja,lat_nat_jja] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res)+1), np.linspace(50,22,int((50-22)/grid_res)+1)) xdiff_lon_ave_jja = np.zeros(lat_nat_jja.shape) ydiff_lat_ave_jja = np.zeros(lat_nat_jja.shape) grid_count_me_jja = np.zeros(lat_nat_jja.shape) if len(angle_jja) > 0: for x in range(0,lon_nat_jja.shape[1]-1): for y in range(0,lat_nat_jja.shape[0]-1): grid_count_me_jja[y,x] =np.nansum((lat_obs<=lat_nat[y,x]) & (lat_obs>lat_nat[y+1,x]) & (lon_obs>lon_nat[y,x]) & (lon_obs<=lon_nat[y,x+1])) if grid_count_me_jja[y,x]>= num_points: log_array_jja = (lat_obs_jja<=lat_nat_jja[y,x]) & (lat_obs_jja>lat_nat_jja[y+1,x]) & (lon_obs_jja>lon_nat_jja[y,x])\ & (lon_obs_jja<=lon_nat_jja[y,x+1]) xdiff_lon_ave_jja[y,x] = np.nanmedian(xdiff_lon_jja[log_array_jja]) ydiff_lat_ave_jja[y,x] = np.nanmedian(ydiff_lat_jja[log_array_jja]) else: pass #finds the magnitude #print(np.sqrt(xdiff_lon_ave[2,1]**2+ydiff_lat_ave[2,1]**2)) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) ma=plt.quiver(lon_nat_jja+(grid_res/2), lat_nat_jja-(grid_res/2), xdiff_lon_ave_jja, ydiff_lat_ave_jja, transform=ccrs.PlateCarree(), zorder = 1,\ units='inches',angles='xy', scale=4,scale_units='inches',width=0.02, linewidths=0.2) plt.title('Summer Day '+str(validday)+' Paired Object Displacement for '+ero_category,fontsize=14) cs = plt.quiverkey(ma, 0.8, 0.28, 2.0,'2 Degrees', labelpos='E',coordinates='figure',transform=ccrs.PlateCarree(), color= 'Red') plt.savefig(working_dir+'arrow_map_displ_jja_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/arrow_map_displ_jja_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #son [lon_nat_son,lat_nat_son] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res)+1), np.linspace(50,22,int((50-22)/grid_res)+1)) xdiff_lon_ave_son = np.zeros(lat_nat_son.shape) ydiff_lat_ave_son = np.zeros(lat_nat_son.shape) grid_count_me_son = np.zeros(lat_nat_son.shape) if len(angle_son) > 0: for x in range(0,lon_nat_son.shape[1]-1): for y in range(0,lat_nat_son.shape[0]-1): grid_count_me_son[y,x] =np.nansum((lat_obs<=lat_nat[y,x]) & (lat_obs>lat_nat[y+1,x]) & (lon_obs>lon_nat[y,x]) & (lon_obs<=lon_nat[y,x+1])) if grid_count_me_son[y,x]>=num_points: log_array_son = (lat_obs_son<=lat_nat_son[y,x]) & (lat_obs_son>lat_nat_son[y+1,x]) & (lon_obs_son>lon_nat_son[y,x])\ & (lon_obs_son<=lon_nat_son[y,x+1]) xdiff_lon_ave_son[y,x] = np.nanmedian(xdiff_lon_son[log_array_son]) ydiff_lat_ave_son[y,x] = np.nanmedian(ydiff_lat_son[log_array_son]) else: pass #finds the magnitude #print(np.sqrt(xdiff_lon_ave[2,1]**2+ydiff_lat_ave[2,1]**2)) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) ma=plt.quiver(lon_nat_son+(grid_res/2), lat_nat_son-(grid_res/2), xdiff_lon_ave_son, ydiff_lat_ave_son, transform=ccrs.PlateCarree(), zorder = 1,\ units='inches',angles='xy', scale=4,scale_units='inches',width=0.02, linewidths=0.2) plt.title('Fall Day '+str(validday)+' Paired Object Displacement for '+ero_category,fontsize=14) cs = plt.quiverkey(ma, 0.8, 0.28, 2.0,'2 Degrees', labelpos='E',coordinates='figure',transform=ccrs.PlateCarree(), color= 'Red') plt.savefig(working_dir+'arrow_map_displ_son_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/arrow_map_displ_son_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #djf [lon_nat_djf,lat_nat_djf] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res)+1), np.linspace(50,22,int((50-22)/grid_res)+1)) xdiff_lon_ave_djf = np.zeros(lat_nat_djf.shape) ydiff_lat_ave_djf = np.zeros(lat_nat_djf.shape) grid_count_me_djf = np.zeros(lat_nat_djf.shape) if len(angle_djf) > 0: for x in range(0,lon_nat_djf.shape[1]-1): for y in range(0,lat_nat_djf.shape[0]-1): grid_count_me_djf[y,x] =np.nansum((lat_obs<=lat_nat[y,x]) & (lat_obs>lat_nat[y+1,x]) & (lon_obs>lon_nat[y,x]) & (lon_obs<=lon_nat[y,x+1])) if grid_count_me_djf[y,x]>=num_points: log_array_djf = (lat_obs_djf<=lat_nat_djf[y,x]) & (lat_obs_djf>lat_nat_djf[y+1,x]) & (lon_obs_djf>lon_nat_djf[y,x])\ & (lon_obs_djf<=lon_nat_djf[y,x+1]) xdiff_lon_ave_djf[y,x] = np.nanmedian(xdiff_lon_djf[log_array_djf]) ydiff_lat_ave_djf[y,x] = np.nanmedian(ydiff_lat_djf[log_array_djf]) else: pass #finds the magnitude #print(np.sqrt(xdiff_lon_ave[2,1]**2+ydiff_lat_ave[2,1]**2)) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) ma=plt.quiver(lon_nat_djf+(grid_res/2), lat_nat_djf-(grid_res/2), xdiff_lon_ave_djf, ydiff_lat_ave_djf, transform=ccrs.PlateCarree(), zorder = 1,\ units='inches',angles='xy', scale=4,scale_units='inches',width=0.02, linewidths=0.2) plt.title('Winter Day '+str(validday)+' Paired Object Displacement for '+ero_category,fontsize=14) cs = plt.quiverkey(ma, 0.8, 0.28, 2.0,'2 Degrees', labelpos='E',coordinates='figure',transform=ccrs.PlateCarree(), color= 'Red') plt.savefig(working_dir+'arrow_map_displ_djf_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/arrow_map_displ_djf_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #mam [lon_nat_mam,lat_nat_mam] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res)+1), np.linspace(50,22,int((50-22)/grid_res)+1)) xdiff_lon_ave_mam = np.zeros(lat_nat_mam.shape) ydiff_lat_ave_mam = np.zeros(lat_nat_mam.shape) grid_count_me_mam = np.zeros(lat_nat_mam.shape) if len(angle_mam) > 0: for x in range(0,lon_nat_mam.shape[1]-1): for y in range(0,lat_nat_mam.shape[0]-1): grid_count_me_mam[y,x] =np.nansum((lat_obs<=lat_nat[y,x]) & (lat_obs>lat_nat[y+1,x]) & (lon_obs>lon_nat[y,x]) & (lon_obs<=lon_nat[y,x+1])) if grid_count_me_mam[y,x]>=num_points: log_array_mam = (lat_obs_mam<=lat_nat_mam[y,x]) & (lat_obs_mam>lat_nat_mam[y+1,x]) & (lon_obs_mam>lon_nat_mam[y,x])\ & (lon_obs_mam<=lon_nat_mam[y,x+1]) xdiff_lon_ave_mam[y,x] = np.nanmedian(xdiff_lon_mam[log_array_mam]) ydiff_lat_ave_mam[y,x] = np.nanmedian(ydiff_lat_mam[log_array_mam]) else: pass #finds the magnitude #print(np.sqrt(xdiff_lon_ave[2,1]**2+ydiff_lat_ave[2,1]**2)) latlon_dims = [-129.8,25.0,-65.0,49.8] [fig,ax] = cartopy_maps.plot_map([latlon_dims[1],latlon_dims[3]], [latlon_dims[0], latlon_dims[2]]) ma=plt.quiver(lon_nat_mam+(grid_res/2), lat_nat_mam-(grid_res/2), xdiff_lon_ave_mam, ydiff_lat_ave_mam, transform=ccrs.PlateCarree(), zorder = 1,\ units='inches',angles='xy', scale=4,scale_units='inches',width=0.02, linewidths=0.2) plt.title('Spring Day '+str(validday)+' Paired Object Displacement for '+ero_category,fontsize=14) cs = plt.quiverkey(ma, 0.8, 0.28, 2.0,'2 Degrees', labelpos='E',coordinates='figure',transform=ccrs.PlateCarree(), color= 'Red') plt.savefig(working_dir+'arrow_map_displ_mam_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/arrow_map_displ_mam_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #with open('readdata.csv', mode='a') as lines: #readdata_lines=csv.writer(lines, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL) #readdata_lines.writerow([str(date),str(valid),str(lat_obs),str(lon_obs),str(magnitude),str(angle)])