Commit 8dadaa4d authored by abuddenberg's avatar abuddenberg

Snapshot Mon Sep 25 14:41:49 EDT 2017

parent e3b30962
......@@ -2,12 +2,29 @@ __author__ = 'abuddenberg'
from config import NA_SCENARIOS_FILES, GLOBAL_PRECIP_FILES
from netCDF4 import Dataset
import numpy as np
import subprocess
import os
OUT_DIR = '../data/forweb/'
#This just keeps getting worse and worse
def get_category(pr_data, stipple, lat, lon):
if stipple[lat, lon] == 1.0:
return 3.0
elif pr_data[lat, lon] == 0.0:
return 1.0
elif not(stipple[lat, lon] == 1.0 or pr_data[lat, lon] == 0.0):
return 2.0
else:
raise Exception('Invalid category')
for scenario in NA_SCENARIOS_FILES:
nc = Dataset(NA_SCENARIOS_FILES[scenario], 'r')
# print nc.variables['lon'][9], nc.variables['lat'][23], nc.variables['MAM_percent_change'][23, 9]
print nc.variables['lon'][9], nc.variables['lat'][23], nc.variables['MAM_percent_change'][23, 9]
with open('../data/pr_{s}_1970-1999_2071-2099.txt'.format(s=scenario.lower().replace('.', '')), 'w') as out_file:
with open(OUT_DIR + 'pr_{s}_1970-1999_2071-2099.txt'.format(s=scenario.lower().replace('.', '').replace(' ', '_')), 'w') as out_file:
print nc.variables.keys()
out_file.write('{0:>8} {1:>14} {2:18} {3:11} {4:18} {5:11} {6:18} {7:11} {8:18} {9:11}\n'.format(
......@@ -24,28 +41,111 @@ for scenario in NA_SCENARIOS_FILES:
))
for lon_idx, lon in enumerate(nc.variables['lon'][:]):
for lat_idx, lat in enumerate(nc.variables['lat'][:]):
djf_change = nc.variables['DJF_percent_change']
mam_change = nc.variables['MAM_percent_change']
jja_change = nc.variables['JJA_percent_change']
son_change = nc.variables['SON_percent_change']
out_file.write(
'{lon:>8} {lat: 14.10f} {djf: 18.3f} {djf_s:11.1f} {mam: 18.3f} {mam_s:11.1f} {jja: 18.3f} {jja_s:11.1f} {son: 18.3f} {son_s:11.1f}\n'.format(
lon=lon,
lat=lat,
djf=nc.variables['DJF_percent_change'][lat_idx, lon_idx],
djf_s=nc.variables['DJF_stipple'][lat_idx, lon_idx],
mam=nc.variables['MAM_percent_change'][lat_idx, lon_idx],
mam_s=nc.variables['MAM_stipple'][lat_idx, lon_idx],
jja=nc.variables['JJA_percent_change'][lat_idx, lon_idx],
jja_s=nc.variables['JJA_stipple'][lat_idx, lon_idx],
son=nc.variables['SON_percent_change'][lat_idx, lon_idx],
son_s=nc.variables['SON_stipple'][lat_idx, lon_idx],
djf=djf_change[lat_idx, lon_idx],
djf_s=get_category(djf_change, nc.variables['DJF_stipple'], lat_idx, lon_idx),
mam=mam_change[lat_idx, lon_idx],
mam_s=get_category(mam_change, nc.variables['MAM_stipple'], lat_idx, lon_idx),
jja=jja_change[lat_idx, lon_idx],
jja_s=get_category(jja_change, nc.variables['JJA_stipple'], lat_idx, lon_idx),
son=son_change[lat_idx, lon_idx],
son_s=get_category(son_change, nc.variables['SON_stipple'], lat_idx, lon_idx),
)
)
nc_out = Dataset(OUT_DIR + 'pr_{s}_1970-1999_2071-2099_percent_change.nc'.format(
s=scenario.lower().replace('.', '').replace(' ', '_')), mode='w')
nc_out.createDimension('lat', size=nc.variables['lat'].size)
nc_out.createDimension('lon', size=nc.variables['lon'].size)
nc_out.contact = 'ken.kunkel@noaa.gov'
nc_out.Conventions = "CF-1.6"
latitude = nc_out.createVariable('lat', 'f8', ('lat',))
latitude.units = "degrees_north"
latitude.long_name = "latitude"
latitude.standard_name = "latitude"
latitude.axis = "Y"
latitude[:] = nc.variables['lat'][:]
longitude = nc_out.createVariable('lon', 'f8', ('lon',))
longitude.modulo = 360.
longitude.long_name = "longitude"
longitude.standard_name = "longitude"
longitude.units = "degrees_east"
longitude.axis = "X"
longitude.topology = "circular"
longitude[:] = nc.variables['lon'][:]
seasons = [
('DJF_percent_change', 'DJF_stipple', 'DJF_category', 'Winter'),
('MAM_percent_change', 'MAM_stipple', 'MAM_category', 'Spring'),
('JJA_percent_change', 'JJA_stipple', 'JJA_category', 'Summer'),
('SON_percent_change', 'SON_stipple', 'SON_category', 'Autumn')
]
for pr_var, stipple_var, category_var, long_name in seasons:
data = nc.variables[pr_var][:]
signif = nc.variables[stipple_var][:]
seasonal_change = nc_out.createVariable(pr_var, 'f8', ('lat', 'lon'))
seasonal_change.long_name = '{} percent change'.format(long_name)
seasonal_change.units = 'Percent'
seasonal_change.missing_value = 1.
seasonal_change[:] = data
seasonal_cat = nc_out.createVariable(category_var, 'f8', ('lat', 'lon'))
seasonal_cat.long_name = '{} category'.format(long_name)
seasonal_cat.units = '1'
seasonal_cat.missing_value = 1e20
seasonal_cat.flag_values = '1.0, 2.0, 3.0'
seasonal_cat.flag_meanings = 'model_disagreement within_historical_variations statistically_significant'
category_array = np.zeros(data.shape)
category_array[np.where(signif == 1.0)] = 3.0 #3.0 denotes areas of statistical significance; red
category_array[np.where(data == 0.0)] = 1.0 #1.0 denotes areas of statistical uncertainty; blue
third_cat_idx = np.where(np.logical_not(np.logical_or(signif == 1.0, data == 0.0)))
category_array[third_cat_idx] = 2.0 #2.0 denotes areas where changes are within historical variations; green
seasonal_cat[:] = category_array
nc_out.close()
subprocess.Popen([
'tar', '-czf', 'Figure_2-14_2-15_files.tar.gz',
'pr_rcp26_1970-1999_2071-2099.txt',
'pr_rcp45_1970-1999_2071-2099.txt',
'pr_rcp60_1970-1999_2071-2099.txt',
'pr_rcp85_1970-1999_2071-2099.txt',
'pr_sres_a1b_1970-1999_2071-2099.txt',
'pr_sres_a2_1970-1999_2071-2099.txt',
'pr_sres_b1_1970-1999_2071-2099.txt',
'pr_rcp26_1970-1999_2071-2099_percent_change.nc',
'pr_rcp45_1970-1999_2071-2099_percent_change.nc',
'pr_rcp60_1970-1999_2071-2099_percent_change.nc',
'pr_rcp85_1970-1999_2071-2099_percent_change.nc',
'pr_sres_a1b_1970-1999_2071-2099_percent_change.nc',
'pr_sres_a2_1970-1999_2071-2099_percent_change.nc',
'pr_sres_b1_1970-1999_2071-2099_percent_change.nc',
'pr_sresa1b_1970-1999_2071-2099_percent_change.nc',
'pr_sresa2_1970-1999_2071-2099_percent_change.nc',
'pr_sresb1_1970-1999_2071-2099_percent_change.nc'
], cwd=os.path.join(os.getcwd(), OUT_DIR))
for file, _unused in GLOBAL_PRECIP_FILES:
nc = Dataset(file, 'r')
print nc.variables['lon'][9], nc.variables['lat'][23], nc.variables['annual_percent_change'][23, 9]
nc = Dataset(file, mode='r')
# print nc.variables['lon'][9], nc.variables['lat'][23], nc.variables['annual_percent_change'][23, 9]
scenario = file.split('_')[3]
with open('../data/pr_{s}_global_1970-1999_2071-2099.txt'.format(s=scenario.lower().replace('.', '')), 'w') as out_file:
with open(OUT_DIR + 'pr_{s}_global_1970-1999_2071-2099.txt'.format(s=scenario.lower().replace('.', '')), 'w') as out_file:
print nc.variables.keys()
out_file.write('{0:>9} {1:>14} {2:21} {3:14}\n'.format(
......@@ -62,13 +162,18 @@ for file, _unused in GLOBAL_PRECIP_FILES:
lon=lon,
lat=lat,
ann=nc.variables['annual_percent_change'][lat_idx, lon_idx],
ann_s=nc.variables['annual_stipple'][lat_idx, lon_idx],
ann_s=get_category(nc.variables['annual_percent_change'], nc.variables['annual_stipple'], lat_idx, lon_idx),
)
)
nc_out = Dataset('../data/pr_{s}_1970-1999_2071-2099_global_percent_change.nc'.format(s=scenario.lower().replace('.', '')), 'w')
data = nc.variables['annual_percent_change'][:]
signif = nc.variables['annual_stipple'][:]
nc_out = Dataset(OUT_DIR + 'pr_{s}_1970-1999_2071-2099_global_percent_change.nc'.format(s=scenario.lower().replace('.', '')), 'w')
nc_out.createDimension('lat', size=nc.variables['lat'].size)
nc_out.createDimension('lon', size=nc.variables['lon'].size)
nc_out.contact = 'ken.kunkel@noaa.gov'
nc_out.Conventions = "CF-1.6"
latitude = nc_out.createVariable('lat', 'f8', ('lat',))
latitude.units = "degrees_north"
......@@ -87,11 +192,35 @@ for file, _unused in GLOBAL_PRECIP_FILES:
longitude[:] = nc.variables['lon'][:]
annual_change = nc_out.createVariable('annual_percent_change', 'f8', ('lat', 'lon'))
annual_change.unit = 'Percent'
annual_change.units = 'Percent'
annual_change.missing_value = 1.
annual_change[:] = nc.variables['annual_percent_change'][:]
annual_change[:] = data
annual_cat = nc_out.createVariable('annual_category', 'f8', ('lat', 'lon'))
annual_cat.long_name = 'annual category'
annual_cat.units = '1'
annual_cat.missing_value = 1e20
annual_cat.flag_values = '1.0, 2.0, 3.0'
annual_cat.flag_meanings = 'model_disagreement within_historical_variations statistically_significant'
category_array = np.zeros(data.shape)
category_array[np.where(signif == 1.0)] = 3.0 #3.0 denotes areas of statistical significance; red
category_array[np.where(data == 0.0)] = 1.0 #1.0 denotes areas of statistical uncertainty; blue
third_cat_idx = np.where(np.logical_not(np.logical_or(signif == 1.0, data == 0.0)))
category_array[third_cat_idx] = 2.0 #2.0 denotes areas where changes are within historical variations; green
annual_cat[:] = category_array
nc_out.close()
new_nc = Dataset('../data/pr_{s}_1970-1999_2071-2099_global_percent_change.nc'.format(s=scenario.lower().replace('.', '')), 'r')
print new_nc.variables['lon'][9], new_nc.variables['lat'][23], new_nc.variables['annual_percent_change'][23, 9]
\ No newline at end of file
new_nc = Dataset(OUT_DIR + 'pr_{s}_1970-1999_2071-2099_global_percent_change.nc'.format(s=scenario.lower().replace('.', '')), 'r')
print new_nc.variables['lon'][9], new_nc.variables['lat'][23], new_nc.variables['annual_category'][23, 9]
subprocess.Popen([
'tar', '-czf', 'Figure_2-6_files.tar.gz',
'pr_rcp26_global_1970-1999_2071-2099.txt',
'pr_rcp85_global_1970-1999_2071-2099.txt',
'pr_rcp26_1970-1999_2071-2099_global_percent_change.nc',
'pr_rcp85_1970-1999_2071-2099_global_percent_change.nc'
], cwd=os.path.join(os.getcwd(), OUT_DIR))
\ No newline at end of file
......@@ -87,10 +87,10 @@ for infilename, outfilename in NA_PRECIP_FILES:
data = np.ma.masked_array(data.filled(3.0)) #3.0 denotes areas of statistical significance; red
data.mask = zeros_mask
data = np.ma.masked_array(data.filled(1.0)) #1.0 denotes areas little change; blue
data = np.ma.masked_array(data.filled(1.0)) #1.0 denotes areas of statistical uncertainty; blue
data.mask = third_cat_mask
data = np.ma.masked_array(data.filled(2.0)) #2.0 denotes areas of statistical uncertainty; green
data = np.ma.masked_array(data.filled(2.0)) #2.0 denotes areas where changes are within historical variationst; green
weird = m.pcolor(x_shifted, y_shifted, data)
m.colorbar(weird, location='right', pad="5%")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment