### Obtain Surface Pressure for specified dates ### ### Create composites of the specified dates ### ### Created by Michael Fischer on 5/10/16 ### ### Edited by Dylan Card on 5/24/16 ### ### Import Libraries ### import numpy as np import netCDF4 as NC4 import datetime as DT import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap import matplotlib.colors as colors import matplotlib as mpl ### Establish dates to composite ### warm_enhance=np.genfromtxt('/home11/ugrad/2014/dc187197/MHC/warm_enhance.csv',delimiter=',') warm=np.genfromtxt('/home11/ugrad/2014/dc187197/MHC/warm_pure.csv',delimiter=',') #warm=np.genfromtxt('/home11/ugrad/2014/dc187197/MHC/warm.csv',delimiter=',') cold=np.genfromtxt('/home11/ugrad/2014/dc187197/MHC/cold.csv',delimiter=',') ### Set up array lengths ### for var in (warm,cold): year=var[0] month=var[1] day=var[2] hour=var[3] print "############ Starting Process Sea Level Pressure ###################" print ### Load a sample year to see dimensions ### datafile = 'http://ramadda.atmos.albany.edu:8080/repository/opendap/Top/UAlbany+Repository/CFSR/2005/pmsl.2005.0p5.anl.nc/entry.das' fin = NC4.Dataset(datafile) lat = np.array(fin.variables['lat']) lon = np.array(fin.variables['lon']) ### Create a limited domain ### def find_closest(array,value): idx = (np.abs(array-value)).argmin() return idx north = 50 south = 35 west = -90 east = -60 lat_north = find_closest(lat[:], north) lat_south = find_closest(lat[:], south) lon_west = find_closest(lon[:], west) lon_east = find_closest(lon[:], east) lat_lim = lat[lat_south:lat_north+1] lon_lim = lon[lon_west:lon_east+1] pmsl = np.zeros((np.size(year),np.size(lat_lim),np.size(lon_lim))) #create the surface pressur array for sn in xrange(len(year)): #This loops over every storm (sn=storm number) year_int = int(year[sn]) month_int = int(month[sn]) day_int = int(day[sn]) hour_int = int(hour[sn]) year_str = str(year_int)#This is the string of the year used to retrieve the data print sn print month_int print day_int print year_int print '---------------------------------------------' pmsl_datafile = 'http://ramadda.atmos.albany.edu:8080/repository/opendap/Top/UAlbany+Repository/CFSR/'+year_str+'/pmsl.'+year_str+'.0p5.anl.nc/entry.das' pmsl_fin = NC4.Dataset(pmsl_datafile) time = pmsl_fin.variables['time'] dates = NC4.num2date(time[:],time.units).tolist() #Get all dates for the current year time_ind = dates.index(DT.datetime(year_int,month_int,day_int,hour_int,0)) #This is the time index #retrieve the surface pressure pmsl[sn,:,:] = pmsl_fin.variables['pmsl'][time_ind,lat_south:lat_north+1,lon_west:lon_east+1] #Close the data files to prevent memory leaks pmsl_fin.close() ### Create the composites ### pmsl = np.divide(pmsl,100) #pressure converted to hPa pmsl = list(pmsl) comp_pmsl = np.mean(pmsl,axis=0) comp_pmsl = np.array(comp_pmsl) ### Plot Composite Pressure ### #Plot Warm Cases if len(var[0])==len(warm[0]): plt.figure(figsize=(11,7)) b = Basemap(projection='cyl', lat_ts=42 ,llcrnrlon=west, urcrnrlon=east,llcrnrlat=south,urcrnrlat=north,resolution='h') #Recall spatial bounds were defined earlier x , y = b(lon_lim,lat_lim) b.drawcoastlines(color='k') b.drawcountries(color='k') b.drawstates(color='k') b.drawcounties(color='k') b.drawparallels(np.arange(south,north,10),labels=[1,0,0,1],fontsize=20) b.drawmeridians(np.arange(west,east,10),labels=[1,0,0,1],fontsize=20) bounds = [1004,1008,1012,1016,1020,1024,1028] #Pick bounds for the colorbar cmap=plt.cm.jet norm = colors.BoundaryNorm(bounds, cmap.N) cont = plt.contourf(x,y,comp_pmsl,bounds) CS = plt.contour(x,y,comp_pmsl,bounds,colors='k') #plt.title('Composite Mean Sea Level Pressure- Cold Cases (hPa) (N=12)',fontsize=14) #To remove trailing zeros class nf(float): def __repr__(self): str = '%.1f' % (self.__float__(),) if str[-1] == '0': return '%.0f' % self.__float__() else: return '%.1f' % self.__float__() # Recast levels to new class CS.levels = [nf(val) for val in CS.levels] # Label levels with specially formatted floats if plt.rcParams["text.usetex"]: fmt = r'%r \%%' else: fmt = '%r' plt.clabel(CS, CS.levels, inline=1, fmt=fmt, fontsize=15) plt.savefig('mslp_warm.png') #Plot Cold Cases if len(var[0])==len(cold[0]): plt.figure(figsize=(11,7)) b = Basemap(projection='cyl', lat_ts=42 ,llcrnrlon=west, urcrnrlon=east,llcrnrlat=south,urcrnrlat=north,resolution='h') #Recall spatial bounds were defined earlier x , y = b(lon_lim,lat_lim) b.drawcoastlines(color='k') b.drawcountries(color='k') b.drawstates(color='k') b.drawcounties(color='k') b.drawparallels(np.arange(south,north,10),labels=[1,0,0,1],fontsize=20) b.drawmeridians(np.arange(west,east,10),labels=[1,0,0,1],fontsize=20) bounds = np.arange(1004,1029,4) #Pick bounds for the colorbar cmap=plt.cm.jet norm = colors.BoundaryNorm(bounds, cmap.N) cont = plt.contourf(x,y,comp_pmsl) CS = plt.contour(x,y,comp_pmsl,colors='k') #plt.title('Composite Mean Sea Level Pressure- Cold Cases (hPa) (N=12)',fontsize=14) #To remove trailing zeros class nf(float): def __repr__(self): str = '%.1f' % (self.__float__(),) if str[-1] == '0': return '%.0f' % self.__float__() else: return '%.1f' % self.__float__() # Recast levels to new class CS.levels = [nf(val) for val in CS.levels] # Label levels with specially formatted floats if plt.rcParams["text.usetex"]: fmt = r'%r \%%' else: fmt = '%r' plt.clabel(CS, CS.levels, inline=1, fmt=fmt, fontsize=15) plt.savefig('mslp_cold.png') #Plot Enhanced Warm Cases if len(var[0])==len(warm_enhance[0]): plt.figure(figsize=(11,7)) b = Basemap(projection='cyl', lat_ts=42 ,llcrnrlon=west, urcrnrlon=east,llcrnrlat=south,urcrnrlat=north,resolution='h') #Recall spatial bounds were defined earlier x , y = b(lon_lim,lat_lim) b.drawcoastlines(color='k') b.drawcountries(color='k') b.drawstates(color='k') b.drawcounties(color='k') b.drawparallels(np.arange(south,north,10),labels=[1,0,0,1]) b.drawmeridians(np.arange(west,east,10),labels=[1,0,0,1]) bounds = np.arange(960,1050,1) #Pick bounds for the colorbar cmap=plt.cm.jet norm = colors.BoundaryNorm(bounds, cmap.N) cont = plt.contourf(x,y,comp_pmsl) CS = plt.contour(x,y,comp_pmsl,colors='k') plt.title('Composite Mean Sea Level Pressure- Enhanced Warm Cases (hPa) (N=17)',fontsize=14) #To remove trailing zeros class nf(float): def __repr__(self): str = '%.1f' % (self.__float__(),) if str[-1] == '0': return '%.0f' % self.__float__() else: return '%.1f' % self.__float__() # Recast levels to new class CS.levels = [nf(val) for val in CS.levels] # Label levels with specially formatted floats if plt.rcParams["text.usetex"]: fmt = r'%r \%%' else: fmt = '%r' plt.clabel(CS, CS.levels, inline=1, fmt=fmt, fontsize=10) plt.savefig('mslp_warm_e.png') #Plot Pure Warm Cases if len(var[0])==len(warm[0]): plt.figure(figsize=(11,7)) b = Basemap(projection='cyl', lat_ts=42 ,llcrnrlon=west, urcrnrlon=east,llcrnrlat=south,urcrnrlat=north,resolution='h') #Recall spatial bounds were defined earlier x , y = b(lon_lim,lat_lim) b.drawcoastlines(color='k') b.drawcountries(color='k') b.drawstates(color='k') b.drawcounties(color='k') b.drawparallels(np.arange(south,north,10),labels=[1,0,0,1]) b.drawmeridians(np.arange(west,east,10),labels=[1,0,0,1]) bounds = np.arange(1004,1029,4) #Pick bounds for the colorbar cmap=plt.cm.jet norm = colors.BoundaryNorm(bounds, cmap.N) cont = plt.contourf(x,y,comp_pmsl) CS = plt.contour(x,y,comp_pmsl,colors='k') #plt.title('Composite Mean Sea Level Pressure- Cold Cases (hPa) (N=12)',fontsize=14) #To remove trailing zeros class nf(float): def __repr__(self): str = '%.1f' % (self.__float__(),) if str[-1] == '0': return '%.0f' % self.__float__() else: return '%.1f' % self.__float__() # Recast levels to new class CS.levels = [nf(val) for val in CS.levels] # Label levels with specially formatted floats if plt.rcParams["text.usetex"]: fmt = r'%r \%%' else: fmt = '%r' plt.clabel(CS, CS.levels, inline=1, fmt=fmt, fontsize=15) plt.savefig('mslp_warm_p.png') ### End of Script ###