##################################################################################### #Code to create grid plots for simple objects. Since hpc-lw-direct machine failed #in early August of 2021, this code is not complete and needs work. # #-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/unmatched' 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_obs = working_dir+'latlon_vday'+str(validday)+'_ERO'+ero_category+'_simpleobj_obs'+'.csv' latlon_mod = working_dir+'latlon_vday'+str(validday)+'_ERO'+ero_category+'_simpleobj_mod'+'.csv' #unmatched initial unmatched_lon_mod=[] unmatched_lat_mod=[] unmatched_lon_obs=[] unmatched_lat_obs=[] unmatched_date_obs=[] unmatched_date_mod=[] unmatched_valid_obs=[] unmatched_valid_mod=[] unmatched_area_obs=[] unmatched_area_mod=[] #unmatched jja unmatched_lon_mod_jja=[] unmatched_lat_mod_jja=[] unmatched_lon_obs_jja=[] unmatched_lat_obs_jja=[] unmatched_date_obs_jja=[] unmatched_date_mod_jja=[] unmatched_valid_mod_jja=[] unmatched_valid_obs_jja=[] unmatched_area_obs_jja=[] unmatched_area_mod_jja=[] #unmatched son unmatched_lon_mod_son=[] unmatched_lat_mod_son=[] unmatched_lon_obs_son=[] unmatched_lat_obs_son=[] unmatched_date_obs_son=[] unmatched_date_mod_son=[] unmatched_valid_mod_son=[] unmatched_valid_obs_son=[] unmatched_area_obs_son=[] unmatched_area_mod_son=[] #unmatched djf unmatched_lon_mod_djf=[] unmatched_lat_mod_djf=[] unmatched_lon_obs_djf=[] unmatched_lat_obs_djf=[] unmatched_date_obs_djf=[] unmatched_date_mod_djf=[] unmatched_valid_obs_djf=[] unmatched_valid_mod_djf=[] unmatched_area_obs_djf=[] unmatched_area_mod_djf=[] #unmatched mam unmatched_lon_mod_mam=[] unmatched_lat_mod_mam=[] unmatched_lon_obs_mam=[] unmatched_lat_obs_mam=[] unmatched_date_obs_mam=[] unmatched_date_mod_mam=[] unmatched_valid_obs_mam=[] unmatched_valid_mod_mam=[] unmatched_area_obs_mam=[] unmatched_area_mod_mam=[] ########################################Read in observation simple objects####################################################################################### target_obs=open(latlon_obs,'r') alllines_obs=target_obs.readlines() for line in alllines_obs: values=line.strip().split(',') if (float(values[2]) == 0): unmatched_date_obs=np.append(unmatched_date_obs,float(values[0][0:8])) unmatched_valid_obs=np.append(unmatched_valid_obs,values[1]) unmatched_lat_obs=np.append(unmatched_lat_obs,float(values[3])) unmatched_lon_obs=np.append(unmatched_lon_obs,float(values[4])) unmatched_area_obs=np.append(unmatched_area_obs,float(values[5])) #jja if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: if (float(values[2]) == 0): unmatched_date_obs_jja=np.append(unmatched_date_obs_jja,float(values[0][0:8])) unmatched_valid_obs_jja=np.append(unmatched_valid_obs_jja,values[1]) unmatched_lat_obs_jja=np.append(unmatched_lat_obs_jja,float(values[3])) unmatched_lon_obs_jja=np.append(unmatched_lon_obs_jja,float(values[4])) unmatched_area_obs_jja=np.append(unmatched_area_obs_jja,float(values[5])) #son if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: if (float(values[2]) == 0): unmatched_date_obs_son=np.append(unmatched_date_obs_son,float(values[0][0:8])) unmatched_valid_obs_son=np.append(unmatched_valid_obs_son,values[1]) unmatched_lat_obs_son=np.append(unmatched_lat_obs_son,float(values[3])) unmatched_lon_obs_son=np.append(unmatched_lon_obs_son,float(values[4])) unmatched_area_obs_son=np.append(unmatched_area_obs_son,float(values[5])) #djf if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: if (float(values[2]) == 0): unmatched_date_obs_djf=np.append(unmatched_date_obs_djf,float(values[0][0:8])) unmatched_valid_obs_djf=np.append(unmatched_valid_obs_djf,values[1]) unmatched_lat_obs_djf=np.append(unmatched_lat_obs_djf,float(values[3])) unmatched_lon_obs_djf=np.append(unmatched_lon_obs_djf,float(values[4])) unmatched_area_obs_djf=np.append(unmatched_area_obs_djf,float(values[5])) #mam if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: if (float(values[2]) == 0): unmatched_date_obs_mam=np.append(unmatched_date_obs_mam,float(values[0][0:8])) unmatched_valid_obs_mam=np.append(unmatched_valid_obs_mam,values[1]) unmatched_lat_obs_mam=np.append(unmatched_lat_obs_mam,float(values[3])) unmatched_lon_obs_mam=np.append(unmatched_lon_obs_mam,float(values[4])) unmatched_area_obs_mam=np.append(unmatched_area_obs_mam,float(values[5])) ##########################################################read in model (ERO) simple objects################################################################### target_mod=open(latlon_mod,'r') alllines_mod=target_mod.readlines() for line in alllines_mod: values=line.strip().split(',') if (float(values[2]) == 0): unmatched_date_mod=np.append(unmatched_date_mod,float(values[0][0:8])) unmatched_valid_mod=np.append(unmatched_valid_mod,values[1]) unmatched_lat_mod=np.append(unmatched_lat_mod,float(values[3])) unmatched_lon_mod=np.append(unmatched_lon_mod,float(values[4])) unmatched_area_mod=np.append(unmatched_area_mod,float(values[5])) #jja if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: if (float(values[2]) == 0): unmatched_date_mod_jja=np.append(unmatched_date_mod_jja,float(values[0][0:8])) unmatched_valid_mod_jja=np.append(unmatched_valid_mod_jja,values[1]) unmatched_lat_mod_jja=np.append(unmatched_lat_mod_jja,float(values[3])) unmatched_lon_mod_jja=np.append(unmatched_lon_mod_jja,float(values[4])) unmatched_area_mod_jja=np.append(unmatched_area_mod_jja,float(values[5])) #son if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: if (float(values[2]) == 0): unmatched_date_mod_son=np.append(unmatched_date_mod_son,float(values[0][0:8])) unmatched_valid_mod_son=np.append(unmatched_valid_mod_son,values[1]) unmatched_lat_mod_son=np.append(unmatched_lat_mod_son,float(values[3])) unmatched_lon_mod_son=np.append(unmatched_lon_mod_son,float(values[4])) unmatched_area_mod_son=np.append(unmatched_area_mod_son,float(values[5])) #djf if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: if (float(values[2]) == 0): unmatched_date_mod_djf=np.append(unmatched_date_mod_djf,float(values[0][0:8])) unmatched_valid_mod_djf=np.append(unmatched_valid_mod_djf,values[1]) unmatched_lat_mod_djf=np.append(unmatched_lat_mod_djf,float(values[3])) unmatched_lon_mod_djf=np.append(unmatched_lon_mod_djf,float(values[4])) unmatched_area_mod_djf=np.append(unmatched_area_mod_djf,float(values[5])) #mam if np.float(values[0][4:6])==6 or np.float(values[0][4:6])==7 or np.float(values[0][4:6])==8: if (float(values[2]) == 0): unmatched_date_mod_mam=np.append(unmatched_date_mod_mam,float(values[0][0:8])) unmatched_valid_mod_mam=np.append(unmatched_valid_mod_mam,values[1]) unmatched_lat_mod_mam=np.append(unmatched_lat_mod_mam,float(values[3])) unmatched_lon_mod_mam=np.append(unmatched_lon_mod_mam,float(values[4])) unmatched_area_mod_mam=np.append(unmatched_area_mod_mam,float(values[5])) ################################################################################################################################################################ 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 unmatched obs scatter plots######################################################################################### #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(unmatched_lon_obs,unmatched_lat_obs, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched PP Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_year_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_year_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #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(unmatched_lon_obs_jja,unmatched_lat_obs_jja, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched PP Summer Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_jja_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_jja_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #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(unmatched_lon_obs_son,unmatched_lat_obs_son, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched PP Fall Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_son_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_son_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #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(unmatched_lon_obs_djf,unmatched_lat_obs_djf, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched PP Winter Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_djf_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_djf_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #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(unmatched_lon_obs_mam,unmatched_lat_obs_mam, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched PP Spring Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_mam_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_mam_obs_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) ########################################################make unmatched mod scatter plots######################################################################### #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(unmatched_lon_mod,unmatched_lat_mod, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched ERO Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_year_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_year_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #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(unmatched_lon_mod_jja,unmatched_lat_mod_jja, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched ERO Summer Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_jja_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_jja_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #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(unmatched_lon_mod_son,unmatched_lat_mod_son, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched ERO Fall Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_son_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_son_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #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(unmatched_lon_mod_djf,unmatched_lat_mod_djf, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched ERO Winter Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_djf_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_djf_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #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(unmatched_lon_mod_mam,unmatched_lat_mod_mam, color='violet',marker='o',s=5,zorder=2,transform=ccrs.PlateCarree()) plt.title('All Unmatched ERO Spring Day '+str(validday)+' Analysis Centroid Location for '+ero_category,fontsize=14) plt.savefig(working_dir+'centroid_scatter_unmatched_mam_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_scatter_unmatched_mam_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #################################calculate and plot frequency map################################################################################################# #np.nanmax if ero_category == 'SLGT': vmax_all_obs=200 vmax_all_mod=20 vmax_season_obs=30 vmax_season_mod=10 if ero_category =='MDT': vmax_all_obs=35 vmax_all_mod=6 vmax_season_obs=25 vmax_season_mod=3 if ero_category=='HIGH': vmax_all_obs=6 vmax_all_mod=6 vmax_season_obs=3 vmax_season_mod=3 #all seasons [lon_nat_obs,lat_nat_obs] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res_freq)+1), np.linspace(50,22,int((50-22)/grid_res_freq)+1)) [lon_nat_mod,lat_nat_mod] = np.meshgrid(np.linspace(-126,-64,int((126-64)/grid_res_freq)+1), np.linspace(50,22,int((50-22)/grid_res_freq)+1)) unmatched_grid_count_obs = np.zeros(lat_nat_obs.shape) unmatched_grid_count_mod = np.zeros(lat_nat_mod.shape) for x in range(0,lon_nat_obs.shape[1]-1): for y in range(0,lat_nat_obs.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_obs[y,x] = np.nansum((unmatched_lat_obs<=lat_nat_obs[y,x]) & (unmatched_lat_obs>lat_nat_obs[y+1,x]) &\ (unmatched_lon_obs>lon_nat_obs[y,x]) & (unmatched_lon_obs<=lon_nat_obs[y,x+1])) for x in range(0,lon_nat_mod.shape[1]-1): for y in range(0,lat_nat_mod.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_mod[y,x] = np.nansum((unmatched_lat_mod<=lat_nat_mod[y,x]) & (unmatched_lat_mod>lat_nat_mod[y+1,x]) &\ (unmatched_lon_mod>lon_nat_mod[y,x]) & (unmatched_lon_mod<=lon_nat_mod[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_obs, lat_nat_obs, unmatched_grid_count_obs, cmap = 'Reds',vmin=0, vmax=vmax_all_obs, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Day '+str(validday)+' Unmatched PP Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_year_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_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_mod, lat_nat_mod, unmatched_grid_count_mod, cmap = 'Reds',vmin=0, vmax=vmax_all_mod, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Day '+str(validday)+' Unmatched ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_year_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_year_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #jja [lon_nat_obs_jja,lat_nat_obs_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)) [lon_nat_mod_jja,lat_nat_mod_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)) unmatched_grid_count_obs_jja = np.zeros(lat_nat_obs_jja.shape) unmatched_grid_count_mod_jja = np.zeros(lat_nat_mod_jja.shape) for x in range(0,lon_nat_obs_jja.shape[1]-1): for y in range(0,lat_nat_obs_jja.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_obs_jja[y,x] = np.nansum((unmatched_lat_obs_jja<=lat_nat_obs_jja[y,x]) & (unmatched_lat_obs_jja>lat_nat_obs_jja[y+1,x]) &\ (unmatched_lon_obs_jja>lon_nat_obs_jja[y,x]) & (unmatched_lon_obs_jja<=lon_nat_obs_jja[y,x+1])) for x in range(0,lon_nat_mod_jja.shape[1]-1): for y in range(0,lat_nat_mod_jja.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_mod_jja[y,x] = np.nansum((unmatched_lat_mod_jja<=lat_nat_mod_jja[y,x]) & (unmatched_lat_mod_jja>lat_nat_mod_jja[y+1,x]) &\ (unmatched_lon_mod_jja>lon_nat_mod_jja[y,x]) & (unmatched_lon_mod_jja<=lon_nat_mod_jja[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_obs_jja, lat_nat_obs_jja, unmatched_grid_count_obs_jja, cmap = 'Reds',vmin=0, vmax=vmax_season_obs, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Summer Day '+str(validday)+' Unmatched PP Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_jja_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_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=ax.pcolormesh(lon_nat_mod_jja, lat_nat_mod_jja, unmatched_grid_count_mod_jja, cmap = 'Reds',vmin=0, vmax=vmax_season_mod, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Summer Day '+str(validday)+' Unmatched ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_jja_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_jja_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #son [lon_nat_obs_son,lat_nat_obs_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)) [lon_nat_mod_son,lat_nat_mod_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)) unmatched_grid_count_obs_son = np.zeros(lat_nat_obs_son.shape) unmatched_grid_count_mod_son = np.zeros(lat_nat_mod_son.shape) for x in range(0,lon_nat_obs_son.shape[1]-1): for y in range(0,lat_nat_obs_son.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_obs_son[y,x] = np.nansum((unmatched_lat_obs_son<=lat_nat_obs_son[y,x]) & (unmatched_lat_obs_son>lat_nat_obs_son[y+1,x]) &\ (unmatched_lon_obs_son>lon_nat_obs_son[y,x]) & (unmatched_lon_obs_son<=lon_nat_obs_son[y,x+1])) for x in range(0,lon_nat_mod_son.shape[1]-1): for y in range(0,lat_nat_mod_son.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_mod_son[y,x] = np.nansum((unmatched_lat_mod_son<=lat_nat_mod_son[y,x]) & (unmatched_lat_mod_son>lat_nat_mod_son[y+1,x]) &\ (unmatched_lon_mod_son>lon_nat_mod_son[y,x]) & (unmatched_lon_mod_son<=lon_nat_mod_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=ax.pcolormesh(lon_nat_obs_son, lat_nat_obs_son, unmatched_grid_count_obs_son, cmap = 'Reds',vmin=0, vmax=vmax_season_obs, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Fall Day '+str(validday)+' Unmatched PP Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_son_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_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=ax.pcolormesh(lon_nat_mod_son, lat_nat_mod_son, unmatched_grid_count_mod_son, cmap = 'Reds',vmin=0, vmax=vmax_season_mod, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Fall Day '+str(validday)+' Unmatched ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_son_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_son_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #djf [lon_nat_obs_djf,lat_nat_obs_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)) [lon_nat_mod_djf,lat_nat_mod_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)) unmatched_grid_count_obs_djf = np.zeros(lat_nat_obs_djf.shape) unmatched_grid_count_mod_djf = np.zeros(lat_nat_mod_djf.shape) for x in range(0,lon_nat_obs_djf.shape[1]-1): for y in range(0,lat_nat_obs_djf.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_obs_djf[y,x] = np.nansum((unmatched_lat_obs_djf<=lat_nat_obs_djf[y,x]) & (unmatched_lat_obs_djf>lat_nat_obs_djf[y+1,x]) &\ (unmatched_lon_obs_djf>lon_nat_obs_djf[y,x]) & (unmatched_lon_obs_djf<=lon_nat_obs_djf[y,x+1])) for x in range(0,lon_nat_mod_djf.shape[1]-1): for y in range(0,lat_nat_mod_djf.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_mod_djf[y,x] = np.nansum((unmatched_lat_mod_djf<=lat_nat_mod_djf[y,x]) & (unmatched_lat_mod_djf>lat_nat_mod_djf[y+1,x]) &\ (unmatched_lon_mod_djf>lon_nat_mod_djf[y,x]) & (unmatched_lon_mod_djf<=lon_nat_mod_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=ax.pcolormesh(lon_nat_obs_djf, lat_nat_obs_djf, unmatched_grid_count_obs_djf, cmap = 'Reds',vmin=0, vmax=vmax_season_obs, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Fall Day '+str(validday)+' Unmatched PP Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_djf_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_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=ax.pcolormesh(lon_nat_mod_djf, lat_nat_mod_djf, unmatched_grid_count_mod_djf, cmap = 'Reds',vmin=0, vmax=vmax_season_mod, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Fall Day '+str(validday)+' Unmatched ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_djf_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_djf_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) #mam [lon_nat_obs_mam,lat_nat_obs_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)) [lon_nat_mod_mam,lat_nat_mod_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)) unmatched_grid_count_obs_mam = np.zeros(lat_nat_obs_mam.shape) unmatched_grid_count_mod_mam = np.zeros(lat_nat_mod_mam.shape) for x in range(0,lon_nat_obs_mam.shape[1]-1): for y in range(0,lat_nat_obs_mam.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_obs_mam[y,x] = np.nansum((unmatched_lat_obs_mam<=lat_nat_obs_mam[y,x]) & (unmatched_lat_obs_mam>lat_nat_obs_mam[y+1,x]) &\ (unmatched_lon_obs_mam>lon_nat_obs_mam[y,x]) & (unmatched_lon_obs_mam<=lon_nat_obs_mam[y,x+1])) for x in range(0,lon_nat_mod_mam.shape[1]-1): for y in range(0,lat_nat_mod_mam.shape[0]-1): # print([lon_nat[y,x],lat_nat[y,x]]) unmatched_grid_count_mod_mam[y,x] = np.nansum((unmatched_lat_mod_mam<=lat_nat_mod_mam[y,x]) & (unmatched_lat_mod_mam>lat_nat_mod_mam[y+1,x]) &\ (unmatched_lon_mod_mam>lon_nat_mod_mam[y,x]) & (unmatched_lon_mod_mam<=lon_nat_mod_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=ax.pcolormesh(lon_nat_obs_mam, lat_nat_obs_mam, unmatched_grid_count_obs_mam, cmap = 'Reds',vmin=0, vmax=vmax_season_obs, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Spring Day '+str(validday)+' Unmatched PP Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_mam_obs_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_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=ax.pcolormesh(lon_nat_mod_mam, lat_nat_mod_mam, unmatched_grid_count_mod_mam, cmap = 'Reds',vmin=0, vmax=vmax_season_mod, transform=ccrs.PlateCarree(), zorder = 1) plt.title('Spring Day '+str(validday)+' Unmatched ERO Centroid Location Frequency for '+ero_category,fontsize=14) plt.colorbar(cs, ax=ax) plt.savefig(working_dir+'centroid_freq_unmatched_mam_mod_'+ero_category+'_day'+str(validday)+'.png',bbox_inches='tight',dpi=150) plt.close() subprocess.call('scp '+working_dir+'/centroid_freq_unmatched_mam_mod_'+ero_category+'_day'+str(validday)+'.png '+FIG_DIR,shell=True) print('year obs') print(np.nanmax(unmatched_grid_count_obs)) print('year mod') print(np.nanmax(unmatched_grid_count_mod)) print('jja obs') print(np.nanmax(unmatched_grid_count_obs_jja)) print('jja mod') print(np.nanmax(unmatched_grid_count_mod_jja)) print('son obs') print(np.nanmax(unmatched_grid_count_obs_son)) print('son mod') print(np.nanmax(unmatched_grid_count_mod_son)) print('djf obs') print(np.nanmax(unmatched_grid_count_obs_djf)) print('djf mod') print(np.nanmax(unmatched_grid_count_mod_djf)) print('mam obs') print(np.nanmax(unmatched_grid_count_obs_mam)) print('mam mod') print(np.nanmax(unmatched_grid_count_mod_mam)) #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)])