Skip to content

Commit

Permalink
FDR analysis
Browse files Browse the repository at this point in the history
  • Loading branch information
zmlabe committed Aug 12, 2019
1 parent 9d94e5b commit 690084d
Show file tree
Hide file tree
Showing 4 changed files with 175 additions and 52 deletions.
Binary file not shown.
Binary file not shown.
88 changes: 36 additions & 52 deletions Scripts/plot_Maps_FDR_Comparison.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
Notes
-----
Author : Zachary Labe
Date : 9 August 2019
Date : 12 August 2019
"""

### Import modules
Expand All @@ -28,7 +28,7 @@
titletime = currentmn + '/' + currentdy + '/' + currentyr
print('\n' '----Plotting Monthly Map Comparison- %s----' % titletime)

### Alott time series (300 ensemble members)
### Alott time series (100 ensemble members)
year1 = 1901
year2 = 2000
years = np.arange(year1,year2+1,1)
Expand All @@ -38,25 +38,21 @@
###############################################################################
### Call arguments
varnames = ['U10','Z50','U200','Z500','SLP','P','T2M','RNET']
experi = np.repeat([r'\textbf{$\bf{\Delta}$Cu}',r'\textbf{$\bf{\Delta}$ATL}',
r'\textbf{$\bf{\Delta}$PAC}'],len(varnames))
experi = np.repeat([r'\textbf{$\bf{\Delta}$Cu}',r'\textbf{$\bf{\Delta}$Pd}',
r'\textbf{$\bf{\Delta}$Pi}'],len(varnames))
letters = list(string.ascii_lowercase)
readallinfo = True
period = 'OND'
period = 'JA'

### Define directories
directorydata = '/seley/zlabe/simu/'
directoryfigure = '/home/zlabe/Desktop/STRATOVARI/Comparison/Regional/%s_Maps/' % period
directoryfigure = '/home/zlabe/Desktop/ANTVARI/%s_Maps/' % period

######################
def readDataPeriods(varnames,simulations,period):
### Call function for 4d variable data
lat,lon,lev,varfuture = MO.readExperiAll(varnames,simulations[0],'surface')
lat,lon,lev,varpastq = MO.readExperiAll(varnames,simulations[1],'surface')

### Only 101 ensembles available for the "Current" simulation (PAMIP-1.1)
varfuture = varfuture[:101,:,:,:]
varpast = varpastq[:101,:,:,:]
lat,lon,lev,varpast = MO.readExperiAll(varnames,simulations[1],'surface')

### Create 2d array of latitude and longitude
lon2,lat2 = np.meshgrid(lon,lat)
Expand All @@ -66,50 +62,38 @@ def readDataPeriods(varnames,simulations,period):
varpast[np.where(varpast <= -1e10)] = np.nan

### Calculate over DJF
if period == 'OND':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,-3:,:,:],axis=1)
varpastm = np.nanmean(varpast[:,-3:,:,:],axis=1)
elif period == 'D':
if period == 'JJA':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,-1:,:,:],axis=1)
varpastm = np.nanmean(varpast[:,-1:,:,:],axis=1)
elif period == 'DJF':
varfuturem = np.nanmean(varfuture[:,5:8,:,:],axis=1)
varpastm = np.nanmean(varpast[:,5:8:,:,:],axis=1)
elif period == 'JAS':
print('Calculating over %s months!' % period)
runs = [varfuture,varpast]
var_mo = np.empty((2,varpast.shape[0]-1,varpast.shape[2],varpast.shape[3]))
for i in range(len(runs)):
var_mo[i,:,:,:] = UT.calcDecJanFeb(runs[i],runs[i],lat,lon,'surface',1)
varfuturem = var_mo[0]
varpastm = var_mo[1]
elif period == 'JFM':
varfuturem = np.nanmean(varfuture[:,6:9,:,:],axis=1)
varpastm = np.nanmean(varpast[:,6:9:,:,:],axis=1)
elif period == 'July':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,0:3,:,:],axis=1)
varpastm = np.nanmean(varpast[:,0:3,:,:],axis=1)
elif period == 'JF':
varfuturem = np.nanmean(varfuture[:,6:7,:,:],axis=1)
varpastm = np.nanmean(varpast[:,6:7:,:,:],axis=1)
elif period == 'August':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,0:2,:,:],axis=1)
varpastm = np.nanmean(varpast[:,0:2,:,:],axis=1)
elif period == 'FMA':
varfuturem = np.nanmean(varfuture[:,7:8,:,:],axis=1)
varpastm = np.nanmean(varpast[:,7:8:,:,:],axis=1)
elif period == 'JA':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,1:4,:,:],axis=1)
varpastm = np.nanmean(varpast[:,1:4,:,:],axis=1)
elif period == 'FM':
varfuturem = np.nanmean(varfuture[:,6:8,:,:],axis=1)
varpastm = np.nanmean(varpast[:,6:8:,:,:],axis=1)
elif period == 'S':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,1:3,:,:],axis=1)
varpastm = np.nanmean(varpast[:,1:3,:,:],axis=1)
elif period == 'J':
varfuturem = np.nanmean(varfuture[:,8:9,:,:],axis=1)
varpastm = np.nanmean(varpast[:,8:9:,:,:],axis=1)
elif period == 'AMJ':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,0:1,:,:],axis=1)
varpastm = np.nanmean(varpast[:,0:1,:,:],axis=1)
elif period == 'F':
varfuturem = np.nanmean(varfuture[:,3:6,:,:],axis=1)
varpastm = np.nanmean(varpast[:,3:6:,:,:],axis=1)
elif period == 'OND':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,1:2,:,:],axis=1)
varpastm = np.nanmean(varpast[:,1:2,:,:],axis=1)
elif period == 'M':
print('Calculating over %s months!' % period)
varfuturem = np.nanmean(varfuture[:,2:3,:,:],axis=1)
varpastm = np.nanmean(varpast[:,2:3,:,:],axis=1)
varfuturem = np.nanmean(varfuture[:,-3:,:,:],axis=1)
varpastm = np.nanmean(varpast[:,-3:,:,:],axis=1)
else:
print(ValueError('Selected wrong month period!'))

Expand Down Expand Up @@ -139,13 +123,13 @@ def readDataPeriods(varnames,simulations,period):
pval = np.empty((3,len(varnames),96,144)) # [variables,simulations,lat,lon]
for v in range(len(varnames)):
diffp,climop,pp,lat,lon,lev = readDataPeriods(varnames[v],
['Future','Current'],
['ANT_Fu','ANT_Cu'],
period)
diffcu,climocu,pcu,lat,lon,lev = readDataPeriods(varnames[v],
['BKsea_Fu','Current'],
['ANT_Cu','ANT_Pi'],
period)
diffsit,climosit,psit,lat,lon,lev = readDataPeriods(varnames[v],
['Osea_Fu','Current'],
['ANT_Fu','ANT_Pi'],
period)

vari[:,v,:,:] = np.asarray([diffp,diffcu,diffsit])
Expand Down Expand Up @@ -219,8 +203,8 @@ def readDataPeriods(varnames,simulations,period):

ax1 = plt.subplot(3,len(varnames),i+1)

m = Basemap(projection='ortho',lon_0=0,lat_0=89,resolution='l',
area_thresh=10000.)
m = Basemap(projection='spstere',boundinglat=-40,lon_0=-180,
resolution='l',round =True)

varn, lons_cyclic = addcyclic(var[i], lon)
varn, lons_cyclic = shiftgrid(180., varn, lons_cyclic, start=False)
Expand Down
139 changes: 139 additions & 0 deletions Scripts/read_MonthlyData.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
"""
Script reads in monthly data from PAMIP experiments using SC-WACCM4
for all 100 ensemble members!
Notes
-----
Author : Zachary Labe
Date : 12 August 2019
Usage
-----
readExperiAll(directory,varid,timeperiod,level)
"""

def readExperiAll(varid,timeperiod,level):
"""
Function reads monthly data from PAMIP simulations for 300 members
Parameters
----------
varid : string
variable name to read
timeperiod: string
time of analysis (Future, Current, Past)
level : string
Height of variable (surface or profile)
Returns
-------
lat : 1d numpy array
latitudes
lon : 1d numpy array
longitudes
var : 4d numpy array or 5d numpy array
[year,month,lat,lon] or [year,month,level,lat,lon]
Usage
-----
lat,lon,lev,var = readExperiAll(varid,timeperiod,level)
"""
print('\n>>>>>>>>>> Using readExperiAll function!')

### Import modules
import numpy as np
from netCDF4 import Dataset

###########################################################################
###########################################################################
###########################################################################
### Directories for Antarctic experiments (1-100 members)
if any([timeperiod=='ANT_Fu',timeperiod=='ANT_Cu',timeperiod=='ANT_Pi']):
if timeperiod == 'ANT_Fu':
experi = 'PAMIP-1.8'
directorydata = '/seley/ypeings/simu/'
totaldirectory = directorydata + experi + '/monthly/'
filename = totaldirectory + varid + '_1900-2000.nc'
print('Reading in Antarctic Future Sea Ice!')
elif timeperiod == 'ANT_Cu':
experi = 'PAMIP-1.1'
directorydata = '/seley/ypeings/simu/'
totaldirectory = directorydata + experi + '/monthly/'
filename = totaldirectory + varid + '_1900-2000.nc'
print('Reading in Antarctic Present-Day Sea Ice!')
elif timeperiod == 'ANT_Pi':
experi = 'PAMIP-1.7'
directorydata = '/seley/ypeings/simu/'
totaldirectory = directorydata + experi + '/monthly/'
filename = totaldirectory + varid + '_1900-2000.nc'
print('Reading in Antarctic Pre-Industrial Sea Ice!')
else:
print(ValueError('Selected wrong time period!'))
else:
print(ValueError('Selected wrong experiment name!'))

if varid == 'EGR' and level == 'surface': # integrated from 500-850 hPa
filename = totaldirectory + varid + '_500_850.nc'

### Read in Data
if level == 'surface': # 3d variables
data = Dataset(filename,'r')
lev = 'surface'
lat = data.variables['latitude'][:]
lon = data.variables['longitude'][:]
varq = data.variables['%s' % varid][:]
data.close()
elif level == 'profile': # 4d variables
data = Dataset(filename,'r')
lev = data.variables['level'][:]
lat = data.variables['latitude'][:]
lon = data.variables['longitude'][:]
varq = data.variables['%s' % varid][:]
data.close()
elif level == 'zonmean': # 3d variables (zonal mean!)
varidz = varid + '_' + level
filename = totaldirectory + varidz + '_1701-2000.nc'
data = Dataset(filename,'r')
lev = data.variables['level'][:]
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
varq = data.variables['%s' % varid][:].squeeze()
data.close()
else:
print(ValueError('Selected wrong height - (surface or profile!)!'))
print('Completed: Read data for *%s* : %s!' % (experi[:],varid))

### Reshape to split years and months
months = 12
if level == 'surface': # 3d variables
var = np.reshape(varq,(varq.shape[0]//months,months,
int(lat.shape[0]),int(lon.shape[0])))
elif level == 'profile': # 4d variables
var = np.reshape(varq,(varq.shape[0]//months,months,int(lev.shape[0]),
int(lat.shape[0]),int(lon.shape[0])))
elif level == 'zonmean': # 3d variables (zonal mean!)
var = np.reshape(varq,(varq.shape[0]//months,months,int(lev.shape[0]),
int(lat.shape[0])))
else:
print(ValueError('Selected wrong height - (surface or profile!)!'))
print('Completed: Reshaped %s array!' % (varid))

### Convert units
if varid in ('TEMP','T2M'):
var = var - 273.15 # Kelvin to degrees Celsius
print('Completed: Changed units (K to C)!')
elif varid == 'SWE':
var = var*1000. # Meters to Millimeters
print('Completed: Changed units (m to mm)!')

print('Completed: Read members 1-100!')

print('>>>>>>>>>> Completed: Finished readExperiAll function!')
return lat,lon,lev,var

#### Test function -- no need to use
#varid = 'T2M'
#timeperiod = 'ANT_Fu'
#level = 'surface'
#lat,lon,lev,var = readExperiAll(varid,timeperiod,level)

0 comments on commit 690084d

Please sign in to comment.