Commit 05862c1c authored by Carl Schreck's avatar Carl Schreck

Automated Nightly Commit - Wed Jan 22 13:40:59 EST 2020

parent 79339e26
""" Draw maps of total nClimGrid data."""
__author__ = "Carl Schreck"
__email__ = "cjschrec@ncsu.edu"
__copyright__ = "Copyright 2019, North Carolina State University"
__license__ = "BSD-3.0"
import xarray as xr
import pandas as pd
import numpy as np
import cartopy.crs as ccrs
import cartopy.io.shapereader as shapereader
import calendar
import os
import cjs # '~carl/lib/python'
import matplotlib.pyplot as plt
cjs.tstamp('Here we go!')
# These are some parameters that could be useful to have up top
var_name = 'case'
if(var_name == 'tmin'):
year = 1983
month = 12
day = 25
levels = list(np.arange(-50, 70, 10))
units = 'Degrees Fahrenheit'
cmap = 'BuPu_r'
file_type = 'scaled'
elif (var_name == 'tmax'):
year = 2006
month = 7
day = 17
levels = list(np.arange(50, 130, 10))
units = 'Degrees Fahrenheit'
cmap = 'OrRd'
file_type = 'scaled'
elif (var_name == 'tavg'):
year = 1989
month = 2
day = 3
levels = list(np.arange(-30, 80, 10))
units = 'Degrees Fahrenheit'
cmap = 'RdBu_r'
file_type = 'scaled'
elif(var_name == 'prcp'):
year = 2015
month = 12
day = 14
levels = [1, 2, 5, 10, 20, 30, 40, 50, 60, 70]
units = 'Millimeters'
cmap = 'GnBu'
file_type = 'scaled'
else:
var_name = 'prcp'
year = 2019
month = 7
day = 14
# levels = list(np.arange(50, 101, 2))
# units = 'Degrees Fahrenheit'
# cmap = 'OrRd'
levels = [1, 2, 5, 10, 20, 30, 40, 50, 60, 70]
units = 'Millimeters'
cmap = 'GnBu'
file_type = 'prelim'
cjs.tstamp('Reading')
file_name = f'{var_name}-{year}{month:02}-grd-{file_type}'
in_path = f'{os.environ["DATA_DIR"]}/nclimgrid/beta/{var_name}/{file_name}.nc'
ds = xr.open_dataset(in_path)
data = ds[var_name].sel(time=f'{year}-{month:02}-{day:02}')
if(var_name != 'prcp'):
data = (data * 9 / 5) + 32
cjs.tstamp('Draw map')
ax = plt.axes(projection=ccrs.AlbersEqualArea(central_longitude=360-95,
standard_parallels=(40,60)))
plt.gca().outline_patch.set_visible(False)
ax.set_extent([-118, -75, 23.5, 50.5])
state_shp_path = f'{os.environ["DATA_DIR"]}' \
f'/geography/ne_10m_admin_1_states_provinces_lakes' \
f'/ne_10m_admin_1_states_provinces_lakes.shp'
for shape in shapereader.Reader(state_shp_path).records():
if shape.attributes['adm0_a3'] == 'USA':
ax.add_geometries(shape.geometry, ccrs.PlateCarree(),
edgecolor='none', facecolor='silver', zorder=1)
ax.add_geometries(shape.geometry, ccrs.PlateCarree(),
edgecolor=(0.2, 0.2, 0.2), linewidth=0.5,
facecolor='none', zorder=3)
cjs.tstamp('Draw data')
plot = data.plot(ax=ax, transform=ccrs.PlateCarree(),
cmap=cmap, levels=levels, zorder=2, extend='both',
vmin=-np.inf, vmax=np.inf,
cbar_kwargs={'orientation': 'horizontal',
'label': units,
'pad': 0.0})
plot.colorbar.set_ticklabels(levels)
plot.colorbar.set_ticks(levels)
plt.title(f'nClimGrid-{file_type} {var_name.upper()}: {day}'
f' {calendar.month_name[month]}'
f' {year}')
plt.savefig(f'figures/nclim-{file_type}-{var_name}-'
f'{year}-{month:02}-{day:02}.png', dpi=200)
plt.show()
cjs.tstamp('Thank you, come again.')
......@@ -14,17 +14,14 @@ import cartopy.io.shapereader as shapereader
import matplotlib as mpl
import cjs # '~carl/lib/python'
from datetime import datetime
import numpy as np
# These are some parameters that could be useful to have up top
year = 2019
month = 10
<<<<<<< HEAD
first_day = 24
last_day = 31
=======
first_day = 6
last_day = 10
>>>>>>> e147b0908d6602919d2fb90d875416404f3f207b
time_range = slice(f'{year}-{month:02}-{first_day:02}',
f'{year}-{month:02}-{last_day:02}')
......@@ -62,6 +59,9 @@ cjs.tstamp('Anomalize')
anom_data = 100. * total_data / clim_data
# Take care of infinities when the observed prcp is nonzero
anom_data = anom_data.where((total_data == 0) | (clim_data != 0), 600)
# and when observed is zero
anom_data = anom_data.where((total_data != 0) | (clim_data != 0), 0)
anom_data = anom_data.where(total_ds['prcp'][0].notnull())
cjs.tstamp('Draw map')
......@@ -84,7 +84,7 @@ for shape in shapereader.Reader(state_shp_path).records():
cjs.tstamp('Draw data')
levels = [5, 10, 25, 50, 75, 100, 125, 150, 200, 300, 500]
plot = anom_data.plot(levels=levels, ax=ax, transform=ccrs.PlateCarree(),
cmap='BrBG', zorder=2,
cmap='BrBG', zorder=2, add_labels=False,
cbar_kwargs={'orientation': 'horizontal',
'label': 'Percent',
'pad': 0.0})
......@@ -106,21 +106,13 @@ today_string = (f'{datetime.now().strftime("Created: %a %b %d %Y")}\n'
plt.text(-0.1, -0.25, today_string,
fontsize=7, ha='left', transform=ax.transAxes)
cjs.tstamp('Add legend')
grey_patch = mpl.patches.Rectangle((0.0, 0.15), 0.1, 0.09, facecolor='silver',
edgecolor='black', linewidth=0.5,
transform=ax.transAxes, zorder=5)
ax.add_patch(grey_patch)
grey_string = ('Climatological and observed\n'
'precipitation were both 0.00"')
plt.text(0.00, 0.06, grey_string, fontsize=7, ha='left', transform=ax.transAxes)
cjs.tstamp('Save figure')
plt.tight_layout()
plt.savefig(f'figures/prcp-pon-{year}-{month:02}{first_day:02}'
f'-{month:02}{last_day:02}.png')
f'-{month:02}{last_day:02}.png', dpi=200)
plt.show()
plt.close()
cjs.tstamp('Thank you, come again.')
""" Read in PRISM ascii files and map them."""
__author__ = "Carl Schreck"
__email__ = "cjschrec@ncsu.edu"
__copyright__ = "Copyright 2019, North Carolina State University"
__license__ = "BSD-3.0"
import xarray as xr
import pandas as pd
import numpy as np
import cartopy.crs as ccrs
import cartopy.io.shapereader as shapereader
import calendar
import os
import cjs # '~carl/lib/python'
import matplotlib.pyplot as plt
cjs.tstamp('Here we go!')
# These are some parameters that could be useful to have up top
var_name = 'ppt'
if(var_name == 'tmin'):
year = 1983
month = 12
day = 25
levels = list(np.arange(-50, 70, 10))
units = 'Degrees Fahrenheit'
cmap = 'BuPu_r'
elif (var_name == 'tmax'):
year = 2006
month = 7
day = 17
levels = list(np.arange(50, 130, 10))
units = 'Degrees Fahrenheit'
cmap = 'OrRd'
elif (var_name == 'tmean'):
year = 1989
month = 2
day = 3
levels = list(np.arange(-30, 80, 10))
units = 'Degrees Fahrenheit'
cmap = 'RdBu_r'
elif(var_name == 'ppt'):
year = 2015
month = 12
day = 14
levels = [1, 2, 5, 10, 20, 30, 40, 50, 60, 70]
units = 'mm'
cmap = 'GnBu'
var_name_to_long_name = {
'tmin': 'Minimum Temperature',
'tmax': 'Maximum Temperature',
'tmean': 'Average Temperature',
'ppt': 'Precipitation',
}
var_name_to_units = {
'tmin': 'degree_Celsius',
'tmax': 'degree_Celsius',
'tmean': 'degree_Celsius',
'ppt': 'Millimeters',
}
cjs.tstamp('Reading')
file_name = f'PRISM_{var_name}_{year}{month:02}{day:02}'
in_path = f'{os.environ["DATA_DIR"]}/prism/{file_name}.txt'
in_data = pd.read_csv(in_path, sep=' ', header=None, skiprows=6,
na_values=-9999, skipinitialspace=True)
in_data = in_data[::-1]
if(var_name != 'ppt'):
in_data = (in_data * 9 / 5) + 32
xllcorner = -125.020833333333
yllcorner = 24.062499999979
cellsize = 0.041666666667
ncols = 1405
nrows = 621
lon = xr.DataArray((xllcorner + (np.arange(ncols) * cellsize)),
dims='lon', name='lon',
attrs={'long_name': 'Longitude', 'units': 'degrees_north'})
lon.assign_coords({'lon': lon})
lat = xr.DataArray((yllcorner + (np.arange(nrows) * cellsize)),
dims='lat', name='lat',
attrs={'long_name': 'Latitude', 'units': 'degrees_east'})
lat.assign_coords({'lat': lat})
data = xr.DataArray(in_data, coords=[lat, lon], name=var_name,
attrs={'long_name': var_name_to_long_name[var_name],
'units': var_name_to_units[var_name]})
cjs.tstamp('Draw map')
ax = plt.axes(projection=ccrs.AlbersEqualArea(central_longitude=360-95,
standard_parallels=(40,60)))
plt.gca().outline_patch.set_visible(False)
ax.set_extent([-118, -75, 23.5, 50.5])
state_shp_path = f'{os.environ["DATA_DIR"]}' \
f'/geography/ne_10m_admin_1_states_provinces_lakes' \
f'/ne_10m_admin_1_states_provinces_lakes.shp'
for shape in shapereader.Reader(state_shp_path).records():
if shape.attributes['adm0_a3'] == 'USA':
ax.add_geometries(shape.geometry, ccrs.PlateCarree(),
edgecolor='none', facecolor='silver', zorder=1)
ax.add_geometries(shape.geometry, ccrs.PlateCarree(),
edgecolor=(0.2, 0.2, 0.2), linewidth=0.5,
facecolor='none', zorder=3)
cjs.tstamp('Draw data')
plot = data.plot(ax=ax, transform=ccrs.PlateCarree(),
cmap=cmap, levels=levels, zorder=2, extend='both',
vmin=-np.inf, vmax=np.inf,
cbar_kwargs={'orientation': 'horizontal',
'label': units,
'pad': 0.0})
plot.colorbar.set_ticklabels(levels)
plot.colorbar.set_ticks(levels)
plt.title(f'PRISM {var_name.upper()}: {day} {calendar.month_name[month]} {year}')
plt.savefig(f'figures/prism-{var_name}-{year}-{month:02}-{day:02}.png', dpi=200)
plt.show()
cjs.tstamp('Thank you, come again.')
......@@ -7,17 +7,70 @@ __license__ = "BSD-3.0"
import xarray as xr
import pandas as pd
import numpy as np
import os
import cjs # '~carl/lib/python'
import argparse
import calendar
cjs.tstamp('Here we go!')
# These are some parameters that could be useful to have up top
parser = argparse.ArgumentParser()
parser.add_argument('--var_name', default='tmin')
parser.add_argument('--year', default=2019)
parser.add_argument('--month', default=9)
args = parser.parse_args()
var_name = args.var_name
year = args.year
month = args.month
cjs.tstamp('Here we go!')
file_name = 'PRISM_tmin_20190919.txt'
in_path = f'{os.environ["DATA_DIR"]}/prism/{file_name}'
var_name_to_long_name = {
'tmin': 'Minimum Temperature',
'tmax': 'Maximum Temperature',
'tave': 'Average Temperature',
'prcp': 'Precipitation',
}
var_name_to_units = {
'tmin': 'degree_Celsius',
'tmax': 'degree_Celsius',
'tave': 'degree_Celsius',
'prcp': 'mm',
}
last_day = calendar.monthrange(year, month)[1]
# Create time array
for day in np.arange(1, last_day+1):
file_name = f'PRISM_{var_name}_{year}{month:02}{day:02}'
in_path = f'{os.environ["DATA_DIR"]}/prism/{file_name}.txt'
out_path = f'{os.environ["DATA_DIR"]}/prism/netcdf/{file_name}.nc'
in_data = pd.read_csv(in_path, sep=' ', header=None, skiprows=6,
na_values=-9999, skipinitialspace=True)
in_data = in_data[::-1]
xllcorner = -125.020833333333
yllcorner = 24.062499999979
cellsize = 0.041666666667
ncols = 1405
nrows = 621
lon = xr.DataArray((xllcorner + (np.arange(ncols) * cellsize)),
dims='lon', name='lon',
attrs={'long_name': 'Longitude', 'units': 'degrees_north'})
lon.assign_coords({'lon': lon})
lat = xr.DataArray((yllcorner + (np.arange(nrows) * cellsize)),
dims='lat', name='lat',
attrs={'long_name': 'Latitude', 'units': 'degrees_east'})
lat.assign_coords({'lat': lat})
data = xr.DataArray(in_data, coords=[lat, lon], name=var_name,
attrs={'long_name': var_name_to_long_name[var_name],
'units': var_name_to_units[var_name]})
data.to_netcdf(out_path)
in_data = pd.read_csv(in_path, sep=' ', header=None, skiprows=6,
na_values=-9999, skipinitialspace=True)
cjs.tstamp('Thank you, come again.')
......@@ -19,14 +19,9 @@ from datetime import datetime
# These are some parameters that could be useful to have up top
# mpl.use('Agg')
year = 2019
month = 10
<<<<<<< HEAD
first_day = 24
last_day = 31
=======
first_day = 11
last_day = 17
>>>>>>> e147b0908d6602919d2fb90d875416404f3f207b
month = 8
first_day = 22
last_day = 30
time_range = slice(f'{year}-{month:02}-{first_day:02}',
f'{year}-{month:02}-{last_day:02}')
......@@ -112,8 +107,9 @@ cjs.tstamp('Save figure')
plt.tight_layout()
plt.savefig(f'figures/tave-anom-{year}-{month:02}{first_day:02}'
f'-{month:02}{last_day:02}.png')
f'-{month:02}{last_day:02}.png', dpi=200)
plt.show()
plt.close()
cjs.tstamp('Thank you, come again.')
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