<jupyter_start><jupyter_text>Creates a Pandas DataFrame of Ocetrac attributes<jupyter_code>import xarray as xr
import ocetrac
file_path_list = ('/burg/home/hs3277/marine-heatwaves/data/preprocessed_oisst.nc',
                  '/burg/home/hs3277/marine-heatwaves/data/ocetrac_labels.nc')

ds = xr.open_mfdataset(file_path_list)
ssta = ds.ssta_notrend
labels = ds.blobs
num_events = labels.attrs['final objects tracked']
print(f"There are {num_events} total events tracked.")
ds
%%time
labels.load();<jupyter_output>CPU times: user 0 ns, sys: 2.42 s, total: 2.42 s
Wall time: 2.53 s<jupyter_text>Section of code in Ocetrac `metrics.py`<jupyter_code>import os, glob
import xarray as xr
import numpy as np
import pandas as pd
import warnings
from skimage.measure import label as label_np
from skimage.measure import regionprops 

try:
    from tqdm.auto import tqdm
    tqdm_avail = True
except:
    warnings.warn(
        "Optional dependency `tqdm` not found. This will make progressbars a lot nicer. \
    Install with `conda install -c conda-forge tqdm`"
    )
    tqdm_avail = False
    
def _wrap(labels):
    ''' Impose periodic boundary and wrap labels, then reorder the labels'''
    first_column = labels[..., 0]
    last_column = labels[..., -1]

    stacked = first_column.stack(z=['time','lat'])
    unique_first = np.unique(stacked[stacked.notnull()])

    # This loop iterates over the unique values in the first column, finds the location of those values in 
    # the first columnm and then uses that index to replace the values in the last column with the first column value
    for i in enumerate(unique_first):
        first = np.where(first_column == i[1])
        last = last_column[first[0], first[1]]
        stacked = last.stack(z=['time','lat'])
        bad_labels = np.unique(stacked[stacked.notnull()])
        replace = np.isin(labels, bad_labels)
        labels = labels.where(replace==False, other=i[1])
    
    labels = labels.fillna(0)
    labels_wrapped = np.unique(labels, return_inverse=True)[1].reshape(labels.shape)
    labels_wrapped = xr.DataArray(labels_wrapped, dims=labels.dims, coords=labels.coords)

    return labels_wrapped

def _get_labels(binary_images):
    '''function used to label binary images at each time step using skimage.measure.label'''
    blobs_labels = label_np(binary_images, background=0)
    return blobs_labels
    
def _get_centroids(sub_labels):
    '''This function uses skimage.measure.regionprops to find the centroids of objects assigned 
    to each unique label'''
    props = regionprops(sub_labels.astype('int'))
    centroids = [(float(sub_labels.lat[round(p.centroid[0])].values),
                  float(sub_labels.lon[round(p.centroid[1])].values)) for p in props]
    for i in range(0,len(centroids)):
        if centroids[i][1] >= 359.875:
            centroids[i] = (centroids[i][0], list(centroids[i])[1] - 359.875)
    
    return centroids

def _get_intensity_area(event, ssta, mhw):
    '''Calculates event intensities and area at each time interval using anomaly data and coordinates 
    cooresponding to the event.'''
    
    event_ssta = ssta.where(event>0, drop=True)        
    mhw['intensity_mean'].append(event_ssta.mean(('lat','lon')).values)
    mhw['intensity_max'].append(event_ssta.max(('lat','lon')).values) 
    mhw['intensity_min'].append(event_ssta.min(('lat','lon')).values)
    mhw['intensity_cumulative'].append(np.nansum(event_ssta))
    coords = event.stack(z=('lat','lon'))
    coord_pairs = [(coords.isel(time=t[0]).dropna(dim='z', how='any').z.lat.values, 
                      coords.isel(time=t[0]).dropna(dim='z', how='any').z.lon.values) for t in enumerate(event.time)]

    mhw['coords'].append(coord_pairs)

    # Calculate weighted cell area assuming 0.25º resolution data and 111 km per degree of latitude
    y, x = zip(*coord_pairs)
    dlon = [np.cos(y[c]*np.pi/180)*(111*.25) for c in np.arange(0, len(coord_pairs))]
    dlat = (111*.25) * np.ones(len(dlon))
    cell_area = [np.sum(dlon[c]*dlat[c]) for c in np.arange(0, len(coord_pairs))]
    mhw['area'].append(cell_area)

    return mhw
    

def _get_metrics(event, ssta):
    '''
    Creates a Pandas DataFrame of event attributes.
    
    Parameters
    ----------
      event : xarray.DataArray   
              Image set containing only objects corresponding to the event of interest. 
              Dimensions should be ('time', 'lat', 'lon')
              
      ssta  : xarray.DataArray
              Temperature vector [1D numpy array of length T]
    
    Returns
    -------
    
    mhw : pandas.DataFrame
          Marine heat wave event attributes. The keys listed below are 
          are contained within the dataset.
 
        'id'                     Unique label given to the MHW [int]
        'date'                   Dates corresponding to the event [datetime format]
        'coords'                 Latitude and longitude of all points contained in the event [tuple(lat,lon)]
        'centroid'               Center of each object contained in the event [tuple(lat,lon)]
        'duration'               Duration of event [months]
        'intensity_max'          Maximum intensity at each time interval [degC]
        'intensity_mean'         Mean intensity at each time interval [degC]
        'intensity_min'          Minimum intensity at each time interval [degC]
        'intensity_cumulative'   Cumulated intensity over the entire event [degC months]
        'area'                   Area of the event at each time interval [km2]
        
    '''
    
    # Initialize dictionary 
    mhw = {}
    mhw['id'] = [] # event label
    mhw['date'] = [] # datetime format
    mhw['coords'] = [] # (lat, lon)
    mhw['centroid'] = []  # (lat, lon)
    mhw['duration'] = [] # [months]
    mhw['intensity_max'] = [] # [deg C]
    mhw['intensity_mean'] = [] # [deg C]
    mhw['intensity_min'] = [] # [deg C]
    mhw['intensity_cumulative'] = [] # [deg C]
    mhw['area'] = [] # [km2]

    # TO ADD:
    # mhw['rate_onset'] = [] # [deg C / month]
    # mhw['rate_decline'] = [] # [deg C / month]

    mhw['id'].append(int(np.nanmedian(event.values)))
    mhw['date'].append(event.time.values.astype('datetime64[M]'))
    mhw['duration'].append(event.time.shape[0])

    # Turn images into binary
    binary_event = event.where(event>=0, other=0)
    binary_event = binary_event.where(binary_event==0, other=1)
      
    sub_labels = xr.apply_ufunc(_get_labels, binary_event,
                                input_core_dims=[['lat', 'lon']],
                                output_core_dims=[['lat', 'lon']],
                                output_dtypes=[binary_event.dtype],
                                vectorize=True,
                                dask='parallelized')
    
    # Turn background to NaNs
    sub_labels = xr.DataArray(sub_labels, dims=binary_event.dims, coords=binary_event.coords)
    sub_labels = sub_labels.where(sub_labels>0, drop=False, other=np.nan) 

    # The labels are repeated each time step, therefore we relabel them to be consecutive
    for p in range(1, sub_labels.shape[0]):
        sub_labels[p,:,:] = sub_labels[p,:,:].values + sub_labels[p-1,:,:].max().values
    
    sub_labels_wrapped = _wrap(sub_labels)
    
    mhw = _get_intensity_area(event, ssta, mhw)
    
    centroid = []
    for s in np.arange(0, sub_labels_wrapped.shape[0]):
        lx = sub_labels_wrapped.isel(time=s)
        east = lx.where(lx.lon < 180, drop=True)
        east['lon'] = np.arange(360.125, 540.125, .25)
        append_east = xr.concat([lx.where(lx.lon >= 180, drop=True), east], dim="lon")
        centroid.append(_get_centroids(append_east))
    mhw['centroid'].append(centroid)
    
    mhw = pd.DataFrame(dict([(name, pd.Series(data)) for name,data in mhw.items()]))
#     mhw.to_csv('df_'+str(mhw['id'].values[0]).zfill(4)+'.csv', index=False)
    return mhw



num_events = int(np.nanmax(labels))
num_events
df.append(df)<jupyter_output><empty_output><jupyter_text>How do we parrelize this instead of looping over event ID<jupyter_code>dataframes = []  
for i in range(1,num_events+1):
    print(i)
    event = labels.where(labels==i, drop=True).reindex({"lon": labels.lon.values})
    df = _get_metrics(event, ssta)
    dataframes.append(df)
# def to_dataframe(blobs, ssta):

#     directory = os.getcwd()+'/ocetrac_csv/'
#     if not os.path.exists(directory):
#         os.makedirs(directory)
#     os.chdir(directory)
    
#     num_events = int(np.nanmax(blobs)) # number of events from Ocetrac
#     print(f'CSV files are saved here: {directory}')
#     for i in tqdm(range(1,num_events+1)):
#         event = blobs.where(blobs==i, drop=True).reindex({"lon": blobs.lon.values})
#         _save_csv(event, ssta)
    
#     csvfiles = glob.glob(os.path.join(directory, '*.csv'))
#     dataframes = []  
#     for csvfile in csvfiles:
#         df = pd.read_csv(csvfile)
#         dataframes.append(df)
    
#     df = pd.concat(dataframes, ignore_index=True)
    
#     for f in glob.glob(os.path.join(directory, '*.csv')):
#         if os.path.exists(f):
#             os.remove(f)
#     os.rmdir(directory)

#     return df
# df = ocetrac.to_dataframe(labels, ssta)
dff = pd.concat(dataframes, ignore_index=True)
dff<jupyter_output><empty_output><jupyter_text>Save metrics as CSV<jupyter_code>dff.to_csv('/burg/home/hs3277/marine-heatwaves/data/ocetrac_metrics.csv', index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13333333333333333, 'fraction_numerical': 0.01216931216931217, 'mean_word_length': 2.989447024060785, 'pattern_counts': {'":': 2, '<': 17, '<?xml version=': 0, '>': 21, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '454778', 'license': 'MIT', 'n_tokens_mistral': 3270, 'n_tokens_neox': 3127, 'n_words': 717, 'path': 'ocetrac_metrics.ipynb', 'repo_name': 'ocetrac/marine-heatwaves', 'size': '9450'}