Commit 201693f2 authored by Carl Schreck's avatar Carl Schreck

Working version of prcp_steps.py

parent ed03b03f
......@@ -25,52 +25,93 @@ scaled_path = f'{base_path}{var_name}-201809-grd-scaled.nc'
extent = [-85.0, -75.0, 32.0, 38.0]
central_longitude = np.mean(extent[0:2])
levels = [0.5, 1, 5, 10, 25, 50, 75, 100, 150, 200, 250]
precip_levels = [0.5, 1, 5, 10, 25, 50, 75, 100, 150, 200]
occur_levels = np.arange(0.1, 1.0, 0.1)
cmap = 'GnBu'
tif_name = 'MSR_50M'
cjs.tstamp('Here we go!')
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(6.4, 9.6),
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(6.4, 8.0),
subplot_kw={
'projection': ccrs.AlbersEqualArea(
central_longitude=central_longitude),
'anchor': 'N'
})
cjs.tstamp('Drawing anusplin occurrence')
ax = axs[0,0]
helpers.draw_station_data(ax, station_path, levels=0.1, cmap=cmap,
data_column=6)
anu_occur_path = f'{base_path}BPRC_20180914_D6_A2.pnt.gz'
anu_occur_data = helpers.read_unscaled(anu_occur_path, var_name, data_column=5)
plot = anu_occur_data.plot(ax=ax, transform=ccrs.PlateCarree(),
add_colorbar=False,
cmap=cmap, levels=occur_levels, zorder=1,
vmin=-np.inf, vmax=np.inf, extend='both',
)
ax.set_title(f'(a) Anusplin occurrence', loc='left')
ax.set_extent(extent)
ax.outline_patch.set_visible(True)
cjs.tstamp('Drawing spheremap occurrence')
ax = axs[0,1]
helpers.draw_station_data(ax, station_path, levels=0.1, cmap=cmap,
data_column=6)
sph_occur_path = f'{base_path}BPRC_20180914_D6_S1.pnt.gz'
sph_occur_data = helpers.read_unscaled(sph_occur_path, var_name, data_column=2)
plot = sph_occur_data.plot(ax=ax, transform=ccrs.PlateCarree(),
add_colorbar=False,
cmap=cmap, levels=occur_levels, zorder=1,
vmin=-np.inf, vmax=np.inf, extend='both',
)
ax.set_title(f'(b) Spheremap occurrence', loc='left')
ax.set_extent(extent)
ax.outline_patch.set_visible(True)
fig.colorbar(plot, ax=axs[0, :], shrink=0.6, orientation='horizontal',
anchor=(0.5, 1.0),
label='Probability of precipitation')
cjs.tstamp('Drawing anusplin precip')
ax = axes[0,0]
helpers.draw_station_data(ax, station_path, levels=levels, cmap=cmap,
ax = axs[1,0]
helpers.draw_station_data(ax, station_path, levels=precip_levels, cmap=cmap,
data_column=6)
anusplin_path = f'{base_path}{var_name.upper()}_20180914_D6_A2.pnt.gz'
anusplin_data = helpers.read_unscaled(anusplin_path, var_name, data_column=5)
plot = anusplin_data.plot(ax=ax, transform=ccrs.PlateCarree(),
add_colorbar=True,
cmap=cmap, levels=levels, zorder=1,
add_colorbar=False,
cmap=cmap, levels=precip_levels, zorder=1,
vmin=-np.inf, vmax=np.inf, extend='both',
cbar_kwargs={'label': 'mm',}
)
ax.set_title(f'(a) Anusplin Precipitation', loc='left')
ax.set_title(f'(c) Anusplin Precipitation', loc='left')
ax.set_extent(extent)
ax.outline_patch.set_visible(True)
cjs.tstamp('Drawing scaled data')
ax = axes[1,1]
helpers.draw_station_data(ax, station_path, levels=levels, cmap=cmap,
ax = axs[1,1]
helpers.draw_station_data(ax, station_path, levels=precip_levels, cmap=cmap,
data_column=6)
scaled_ds = xr.open_dataset(scaled_path)
scaled_data = scaled_ds[var_name].sel(time=time).drop('time')
plot = scaled_data.plot(ax=ax, transform=ccrs.PlateCarree(),
add_colorbar=True,
cmap=cmap, levels=levels, zorder=1,
add_colorbar=False,
cmap=cmap, levels=precip_levels, zorder=1,
vmin=-np.inf, vmax=np.inf, extend='both',
cbar_kwargs={'label': 'mm',}
)
ax.set_title(f'(f) Scaled data', loc='left')
ax.set_title(f'(d) Scaled data', loc='left')
ax.set_extent(extent)
ax.outline_patch.set_visible(True)
str_levels = [str(lev) for lev in precip_levels]
cb = fig.colorbar(plot, ax=axs[1, :], shrink=0.6, orientation='horizontal',
anchor=(0.5, 1.0),
label='Precipitation (mm)')
cb.ax.set_xticklabels(str_levels)
cjs.tstamp('Drawing figure')
......
......@@ -33,7 +33,7 @@ tif_name = 'MSR_50M'
cjs.tstamp('Here we go!')
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(6.4, 9.6),
fig, axs = plt.subplots(nrows=3, ncols=1, figsize=(6.4, 9.6),
subplot_kw={
'projection': ccrs.AlbersEqualArea(
central_longitude=central_longitude),
......@@ -41,7 +41,7 @@ fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(6.4, 9.6),
})
i = 0
for ax, var_name in zip(axes[:2], var_names):
for ax, var_name in zip(axs[:2], var_names):
cjs.tstamp('Drawing scaled data')
station_path = f'{base_path}{var_name.upper()}19890203map.txt'
unscaled_path = f'{base_path}{var_name.upper()}_19890203_D6_A2.pnt.gz'
......@@ -65,7 +65,7 @@ for ax, var_name in zip(axes[:2], var_names):
# cjs.add_shaded_relief(ax, tif_name=tif_name)
cjs.tstamp('Adding elevation plot')
ax = axes[2]
ax = axs[2]
elev_path = f'{os.environ["DATA_DIR"]}/geography/elev.americas.5-min.nc'
elev_ds = xr.open_dataset(elev_path)
levels = np.arange(800, 3800, 100)
......
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