Commit 4e9a2817 authored by Carl Schreck's avatar Carl Schreck

mapping prism differences

parent 7a29d214
#!/bin/bash --login
#if [ -z "$1" ]; then
# QUEUE=`pick_queue`
#else
# QUEUE=$1
#fi
QUEUE=allnodes
echo $QUEUE
PY_DIR=`pwd`
PY_SCRIPT=calc_prism_diff
LOG_DIR=$PY_DIR/log
mkdir -p $LOG_DIR
for month in {1..12}
do
# for var_name in tmin tmax tavg prcp
for var_name in pgt0
do
JOB_NAME=$PY_SCRIPT"_"$var_name"_"$month
echo $JOB_NAME `date`
LOG_FILE=$LOG_DIR/$JOB_NAME.log
ERR_FILE=$LOG_DIR/$JOB_NAME.err
PY_OPTION="--var_name "$var_name" --month "$month
rm $LOG_FILE $ERR_FILE
bsub \
-J $JOB_NAME \
-o $LOG_FILE \
-e $ERR_FILE \
-q $QUEUE \
-n 1 -R "span[hosts=1]" -W 48:00 \
-sp 75 \
~/template/run_py.sh $PY_DIR $PY_SCRIPT "$PY_OPTION"
done
done
""" For a given month, calculate the mean bias/difference at each point."""
__author__ = "Carl Schreck"
__email__ = "cjschrec@ncsu.edu"
__copyright__ = "Copyright 2020, North Carolina State University"
__license__ = "BSD-3.0"
import xarray as xr
import sys
import numpy as np
import argparse
import os
import cjs # '~carl/lib/python'
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='pgt0')
parser.add_argument('--month', default=9)
args = parser.parse_args()
var_name = args.var_name
month = int(args.month)
cjs.tstamp(f'{var_name} -- {month}')
years = range(1981,2018+1)
ncg_var_to_prism_var = {
'tmin': 'tmin',
'tmax': 'tmax',
'tavg': 'tmean',
'prcp': 'ppt',
}
if var_name == 'pgt0': # precip greater than zero
ncg_var = 'prcp'
else:
ncg_var = var_name
prism_var = ncg_var_to_prism_var[ncg_var]
for year in years:
ncg_path = f'{os.environ["DATA_DIR"]}/nclimgrid/beta/{ncg_var}/' \
f'{ncg_var}-{year}{month:02}-grd-scaled.nc'
prism_path = f'{os.environ["DATA_DIR"]}/prism/netcdf/' \
f'PRISM_{prism_var}_{year}{month:02}.nc'
ncg_ds = xr.open_dataset(ncg_path)
prism_ds = xr.open_dataset(prism_path)
sub_prism = prism_ds.sel(lat=ncg_ds.lat, lon=ncg_ds.lon, method='nearest')
sub_prism['lat'] = ncg_ds.lat
sub_prism['lon'] = ncg_ds.lon
if var_name == 'pgt0': # precip greater than zero
ncg_ds[ncg_var] = xr.where(ncg_ds[ncg_var] == 0, 0, 1)
sub_prism[prism_var] = xr.where(sub_prism[prism_var] == 0, 0, 1)
if year == years[0]:
cjs.tstamp(year)
bias = ((ncg_ds[ncg_var] - sub_prism[prism_var]).mean(dim='time')
/ len(years))
absdiff = (abs(ncg_ds[ncg_var] - sub_prism[prism_var]).mean(dim='time')
/ len(years))
else:
cjs.tstamp(year)
bias += ((ncg_ds[ncg_var] - sub_prism[prism_var]).mean(dim='time')
/ len(years))
absdiff += (abs(ncg_ds[ncg_var] - sub_prism[prism_var]).mean(dim='time')
/ len(years))
cjs.tstamp('Writing')
bias.attrs['long_name'] = 'Bias ncg - prism'
absdiff.attrs['long_name'] = 'Absolute difference ncg - prism'
if var_name == 'pgt0':
bias.attrs['units'] = 'frequency'
absdiff.attrs['units'] = 'frequency'
else:
bias.attrs['units'] = ncg_ds[ncg_var].attrs['units']
absdiff.attrs['units'] = ncg_ds[ncg_var].attrs['units']
ds = xr.Dataset(data_vars={'bias': bias, 'absdiff': absdiff})
out_path = f'{os.environ["DATA_DIR"]}/nclimgrid/beta/compare/' \
f'{var_name}-prism-{month:02}.nc'
ds.to_netcdf(out_path)
cjs.tstamp('Thank you, come again.')
""" Draw a map of comparison between nclimgrid and prism."""
__author__ = "Carl Schreck"
__email__ = "cjschrec@ncsu.edu"
__copyright__ = "Copyright 2020, North Carolina State University"
__license__ = "BSD-3.0"
import xarray as xr
import numpy as np
import os
import glob
from matplotlib import pyplot as plt
import cartopy.crs as ccrs
import cjs # '~carl/lib/python'
cjs.tstamp('Here we go!')
# These are some parameters that could be useful to have up top
var_name = 'tavg'
in_path = f'{os.environ["DATA_DIR"]}/nclimgrid/beta/compare/{var_name}*.nc'
var_name_to_contour_level = {
'tmin': 0.5,
'tmax': 0.5,
'tavg': 0.5,
'prcp': 0.5,
'pgt0': 0.5,
}
ds = xr.open_mfdataset(glob.glob(in_path), concat_dim='month', combine='nested')
ds = ds.mean(dim='month', keep_attrs=True)
ax_bias = plt.subplot(211,
projection=ccrs.AlbersEqualArea(
central_longitude=360-95,
standard_parallels=(40,60)))
plt.gca().outline_patch.set_visible(False)
ax_bias.set_extent([-118, -75, 23.5, 50.5])
shp_name = 'ne_10m_admin_1_states_provinces_lakes'
cjs.add_map(ax_bias, shp_name, edgecolor='none', facecolor='silver', zorder=1,
usa_only=True)
cjs.add_map(ax_bias, shp_name, edgecolor=(0.2, 0.2, 0.2), facecolor='none',
zorder=3, usa_only=True)
levels = list(np.arange(-5,6)*0.5)
ds['bias'].plot(ax=ax_bias, transform=ccrs.PlateCarree(),
cmap='RdBu_r', levels=levels, zorder=2, extend='both',
vmin=-np.inf, vmax=np.inf,
)
ax_bias.set_title('a) Bias', loc='left')
ax_bias.set_title(f'mean: {ds["bias"].mean().values:0.2f}', loc='right')
cjs.tstamp(f'mean bias: {ds["bias"].mean().values:0.2f}')
ax_absdiff = plt.subplot(212,
projection=ccrs.AlbersEqualArea(
central_longitude=360-95,
standard_parallels=(40,60)))
plt.gca().outline_patch.set_visible(False)
ax_absdiff.set_extent([-118, -75, 23.5, 50.5])
shp_name = 'ne_10m_admin_1_states_provinces_lakes'
cjs.add_map(ax_absdiff, shp_name, edgecolor='none', facecolor='silver',
zorder=1, usa_only=True)
cjs.add_map(ax_absdiff, shp_name, edgecolor=(0.2, 0.2, 0.2), facecolor='none',
zorder=3, usa_only=True)
levels = list(np.arange(1,8)*0.5)
ds['absdiff'].plot(ax=ax_absdiff, transform=ccrs.PlateCarree(),
levels=levels, zorder=2, extend='both',
vmin=-np.inf, vmax=np.inf,
)
ax_absdiff.set_title('b) Abs. diff.', loc='left')
ax_absdiff.set_title(f'mean: {ds["absdiff"].mean().values:0.2f}', loc='right')
cjs.tstamp(f'mean absdiff: {ds["absdiff"].mean().values:0.2f}')
plt.suptitle(var_name)
plt.savefig(f'figures/{var_name}.png', dpi=200)
plt.show()
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