<jupyter_start><jupyter_text>**Load python libraries**<jupyter_code>%matplotlib inline
import datetime
import matplotlib.pyplot as plt
import os.path
import xarray as xr
import pandas as pd
import numpy as np
import requests
import netCDF4
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
from IPython.display import Image, display
from dask.distributed import Client, progress
from glob import glob
import colorsys
from matplotlib.colors import LinearSegmentedColormap # Linear interpolation for color maps
import wrf
from scipy import stats
from statistics import median
latitude=34.5
longitude=359.375-119.375

filename =  '../data/ar_catalog/globalARcatalog_MERRA2_1980-2019_v3.0.nc'
ds = xr.open_dataset(filename, engine='netcdf4')
ds = ds.squeeze()

ds_SB = ds.sel(lat=latitude, lon=longitude)
ds_SB = ds_SB.squeeze()
print('ds size in GB {:0.2f}\n'.format(ds.nbytes / 1e9))
ds.info()
def get_ivt_direction(ds1, ds2):
    da = ds1.shape.dropna(dim='time') ## drop the times where there is no AR present
    idx_times = np.unique(da['time']) ## list of days where AR is present in loc
    times_list = pd.DatetimeIndex(idx_times) ## create a pandas datetime index for easy time use
    ds_new = ds2.sel(time=times_list) # get only AR days

    ## get the ivt direction of the AR object that is over SB
    ivtdir = []
    idx_shape = da.values.astype(int) # get integer index values of ar object
    for i in range(len(times_list)):
        ivtdir.append(ds_new.ivtdir[i].isel(lon=idx_shape[i]).values)
    ivtdir = np.asarray(ivtdir)
    
    return(ivtdir)
## get ivt direction for all AR events in SB
ivtdir = get_ivt_direction(ds_SB, ds)
print(ivtdir.shape)

## get ivt direction for AR events in march only
def is_mar(month):
    return (month == 3)

ds_march = ds_SB.sel(time=is_mar(ds['time.month']))
ivtdir_march = get_ivt_direction(ds_march, ds)
print(ivtdir_march.shape)
## Get direction of specific event ivt direction
ds_event = ds.sel(time='2019-03-05 12:00')
da = ds_SB.shape.dropna(dim='time') ## drop the times where there is no AR present
idx_event = da.sel(time= '2019-03-05 12:00')
ivtdir_event = ds_event.ivtdir.isel(lon=idx_event.astype(int))
print('IVT Direction during March AR Event:', ivtdir_event.values)
results = [ivtdir, ivtdir_march]
n_bins = 100
fig, axs = plt.subplots(1, 1, tight_layout=True)
subtitles = ['All ARs in SB', 'March ARs in SB']

# We can set the number of bins with the `bins` kwarg
 
axs.hist(results[0], bins=n_bins, facecolor='blue')
axs.hist(results[1], bins=n_bins, facecolor='grey')
axs.axvline(ivtdir_event.values, color='r', linestyle='-', linewidth=2)
axs.set_xlim(0, 360)
axs.set_xlabel('Direction of Mean IVT (degree)', fontsize=12)
axs.set_ylabel('Frequency', fontsize=12)
axs.title.set_text('IVT Direction for ARs in Santa Barbara')<jupyter_output><empty_output><jupyter_text>Import MERRA2 0* Isotherm Data for AR Days only<jupyter_code>filename = '/home/sbarc/students/nash/data/MERRA-2/SB_AR_0_degree_isotherm.nc'
ds = xr.open_dataset(filename)<jupyter_output><empty_output><jupyter_text>Check for 0* crossings in the temperature profile between 1000 and 200 hPa.If single zero crossing, altitude is taken as the freezing level.Special case 1: No zero crossing (T < 0degC throughout column)Flag freezing level as missing.Special case 2: Multiple zero crossing due to temperature inversionFlag all locations, store only lowest Z0 value<jupyter_code>ds.T.sel(lat=30.0, lon=-125.0, time='1980-01-09').values - 273.15
ht = ds.H
t = ds.T - 273.15 # convert to Kelvin

interp_results = wrf.interplevel(ht, t, 0.0, missing=np.nan)
interp_results
interp_results.lon
## select SB
latitude=34.5
longitude=-119.375

SB_results = interp_results.sel(lat=latitude, lon=longitude)

## select only March events
def is_mar(month):
    return (month == 3)

ds_march = interp_results.sel(time=is_mar(interp_results['time.month']), lat=latitude, lon=longitude)
ds_march
def preprocess(ds):
    '''keep only the selected lats and lons'''
    return ds.sel(latitude=latitude, longitude=360.-119.5)

filename_pattern =  '/home/sbarc/students/nash/data/CFSv2/netcdf_CFSv2_vertical_levels/201*'
ds = xr.open_mfdataset(filename_pattern, engine='netcdf4', preprocess=preprocess)
print('ds size in GB {:0.2f}\n'.format(ds.nbytes / 1e9))
ds = ds.sel(time=slice('2019-03-06-12', '2019-03-06-18'))
zero_isotherm = ds.hgt_0Cisotherm
print(zero_isotherm.values)
zeroiso_event = np.mean(zero_isotherm.values)
print(np.nanmean(SB_results.values))
# SB_results.values
print(np.nanmean(ds_march.values))
len(SB_results)
results = [SB_results.values, ds_march.values]
n_bins = 100
fig, axs = plt.subplots(1, 1, tight_layout=True)
subtitles = ['All ARs in SB', 'March ARs in SB']

# We can set the number of bins with the `bins` kwarg
 
axs.hist(results[0], bins=n_bins, facecolor='blue')
axs.hist(results[1], bins=n_bins, facecolor='grey')
axs.axvline(zeroiso_event, color='r', linestyle='-', linewidth=2)
axs.axvline(np.nanmean(SB_results.values), color='k', linestyle='--', linewidth=2)
axs.axvline(np.nanmean(ds_march.values), color='grey', linestyle='--', linewidth=2)
# axs.set_xlim(0, 360)
axs.set_xlabel('Height of 0*C Isotherm', fontsize=12)
axs.set_ylabel('Frequency', fontsize=12)
axs.title.set_text('0*C Isotherm Height for ARs in Santa Barbara')
## Compare # ARs in Cannon et al., 2017 to Guan and Waliser 2015
latitude=34.5
longitude=359.375-119.375

filename =  '/home/sbarc/students/nash/data/ar_catalog/globalARcatalog_MERRA2_1980-2019_v2.0.nc'
ds = xr.open_dataset(filename, engine='netcdf4')
ds = ds.squeeze()

## select only March events
def is_ondjfm(month):
    return (month >= 10) | (month <=3)

ds_cannon = ds.sel(time=slice('2014-10-01', '2017-03-31'), lat=slice(20., 60.), lon=slice(360-160, 360-110))
ds_cannon = ds_cannon.sel(time=is_ondjfm(ds_cannon['time.month']))
ds_cannon = ds_cannon.squeeze()
print('ds size in GB {:0.2f}\n'.format(ds.nbytes / 1e9))
ds_cannon

da = ds_cannon.shape.count(dim='time')
da.min()
da.sel(lat=latitude, lon=longitude)
da.max()
## select only months between november and april
def is_ndjfma(month):
    return (month >= 11) | (month <=4)

ds_rutz = ds.sel(time=slice('1988-11-01', '2011-04-30'), lat=latitude, lon=longitude)
ds_rutz = ds_rutz.sel(time=is_ndjfma(ds_rutz['time.month']))
print('ds size in GB {:0.2f}\n'.format(ds.nbytes / 1e9))
ds_rutz
da_rutz = ds_rutz.shape.count(dim='time')
da_rutz<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1332925710791807, 'fraction_numerical': 0.04356465912564965, 'mean_word_length': 5.61577350859454, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '649959', 'license': 'MIT', 'n_tokens_mistral': 2681, 'n_tokens_neox': 2461, 'n_words': 586, 'path': 'analysis/MERRA2 0deg isotherm climatology.ipynb', 'repo_name': 'dlnash/arthunderstorm2019', 'size': '6542'}