<jupyter_start><jupyter_text>OpenDrift Documentation[https://github.com/opendrift/opendrift/wiki](https://github.com/opendrift/opendrift/wiki)*** Installation`git clone https://github.com/OpenDrift/opendrift.git``python setup.py develop --user`***<jupyter_code>import numpy as np
import xarray as xr
import netCDF4 as nc
import matplotlib.pyplot as plt
import os
from glob import glob
from datetime import datetime, timedelta
from dateutil.parser import parse
from itertools import repeat
from salishsea_tools import viz_tools
from rotate_fields_python import rotate_fields
from tqdm import tqdm_notebook as tqdm
from IPython.display import HTML

from opendrift.readers import reader_netCDF_CF_generic
from opendrift.models.oceandrift import OceanDrift

%matplotlib inline
plt.rcParams['font.size'] = 12
plt.rcParams['animation.html'] = 'html5'<jupyter_output><empty_output><jupyter_text>*** Prepare forcing filesForcing must be gridded on lon/lat coordinates. The easiest way is to reshape everything onto a single spatial dimension and use `opendrift.readers.reader_netCDF_CF_unstructured`.<jupyter_code># Define paths
paths = {
    'NEMO': '/results/SalishSea/hindcast.201812',
    'NEMO_cutoff': '/results2/SalishSea/hindcast.201812_annex',
    'HRDPS': '/results/forcing/atmospheric/GEM2.5/operational',
    'out': '/data/bmoorema/results/opendrift',
    'date_cutoff': '2016 Nov 21',
}
mask = xr.open_dataset('/data/bmoorema/MEOPAR/grid/mesh_mask201702.nc')
coords = xr.open_dataset('/data/bmoorema/MEOPAR/grid/coordinates_seagrid_SalishSea201702.nc', decode_times=False)
daterange = [parse(d) for d in ['2017 Nov 10', '2017 Nov 12']]<jupyter_output><empty_output><jupyter_text>Prefix constructor for SalishSeaCast<jupyter_code>def make_prefix(date, paths, res='h'):
    """Construct path prefix for local SalishSeaCast results given date object and paths dict
    e.g., /results/SalishSea/hindcast.201812/ddmmmyy/SalishSea_1h_yyyymmdd_yyyymmdd
    """

    path = paths['NEMO']
    if 'date_cutoff' in paths and date >= parse(paths['date_cutoff']):
        path = paths['NEMO_cutoff']
    fn = '_'.join([f'SalishSea_1{res}', *repeat(date.strftime('%Y%m%d'), 2)])
    prefix = os.path.join(path, date.strftime('%d%b%y').lower(), fn)
    
    return prefix<jupyter_output><empty_output><jupyter_text>Calculate approximate lon/lat coordinate arrays to give similar spacing to NEMO grid_X, grid_Y<jupyter_code>def calculate_coordinate_arrays(coords):
    """Calculate approximate lon and lat indices to give similar spacing as NEMO grid_X, grid_Y
    """

    lonlat = []
    for var1, var2, axis, fac in zip(['glamt', 'gphit'], ['e1t', 'e2t'], ['x', 'y'], [np.cos(49 * np.pi / 180), 1]):
        deg2m = fac * 111e3
        n = coords.dims[axis] - 1
        coord, spacing = coords[var1].values, mask[var2].values.mean()
        coord_0 = (coord.max() + coord.min()) / 2 - (n / 2) * spacing / deg2m
        lonlat.append(coord_0 + np.arange(n) * spacing / deg2m)

    return lonlat<jupyter_output><empty_output><jupyter_text>Process SalishSeaCast<jupyter_code>def process_NEMO(date, paths):
    """
    """
    
    # Process input
    prefix = make_prefix(date, paths)
    NEMO_out = os.path.join(paths['out'], 'forcing', os.path.split(prefix)[-1] + '_opendrift.nc')
    if os.path.exists(NEMO_out):
        print(f'{NEMO_out} already exists! Skipping ...')
        return
    else:
        print(f'Processing {NEMO_out}')
        raw = []
        for k, key in zip(['U', 'V'], ['vozocrtx', 'vomecrty']):
            with xr.open_dataset(prefix + f'_grid_{k}.nc') as data:
                raw.append(data[key][:, 0, ...].values)
                if k is 'U': time = data.time_counter.values.astype('datetime64[s]').astype(datetime)
        u, v = viz_tools.unstagger(*raw)
        lon, lat = calculate_coordinate_arrays(coords)
        
        # Write to netCDF out
        with nc.Dataset(NEMO_out, 'w') as ds:
            
            # Create time dimension
            time_units = 'seconds since 1900-01-01 00:00:00'
            ds.createDimension('time', None)
            tempvar = ds.createVariable('time', float, 'time', zlib=True)
            tempvar[:] = nc.date2num(time, units=time_units, calendar='gregorian')
            tempvar.standard_name = 'time'
            tempvar.calendar = 'gregorian'
            tempvar.units = time_units
            
            # Create lon/lat
            dims = ('latitude', 'longitude')
            for var, name, n in zip([lat, lon], dims, u[0, ...].shape):
                ds.createDimension(name, n)
                tempvar = ds.createVariable(name, float, name, zlib=True, fill_value=1e20)
                tempvar[:] = var
                tempvar.standard_name = name
            
            # Create velocities
            for var, name in zip([u, v], ['x_sea_water_velocity', 'y_sea_water_velocity']):
                tempvar = ds.createVariable(name, float, ('time',) + dims, zlib=True, fill_value=1e20)
                tempvar[:] = var
                tempvar.standard_name = name
def process_mask(paths, mask):
    """
    """
    
    # Process input
    mask_out = os.path.join(paths['out'], 'forcing', 'meshmask_opendrift.nc')
    if os.path.exists(mask_out):
        print(f'{mask_out} already exists! Skipping ...')
        return
    else:
        print(f'Processing {mask_out}')
        tmask = mask.tmask[0, 0, 1:, 1:].values
        lon, lat = calculate_coordinate_arrays(coords)
        
        # Write to netCDF out
        with nc.Dataset(mask_out, 'w') as ds:
            
            # Create time dimension
            time_units = 'seconds since 1900-01-01 00:00:00'
            ds.createDimension('time', None)
            tempvar = ds.createVariable('time', float, 'time', zlib=True)
            tempvar[:] = nc.date2num(datetime(2019, 1, 1, 0, 0, 0), units=time_units, calendar='gregorian')
            tempvar.standard_name = 'time'
            tempvar.calendar = 'gregorian'
            tempvar.units = time_units
            
            # Create lon/lat
            dims = ('latitude', 'longitude')
            for var, name, n in zip([lat, lon], dims, tmask.shape):
                ds.createDimension(name, n)
                tempvar = ds.createVariable(name, float, name, zlib=True, fill_value=1e20)
                tempvar[:] = var
                tempvar.standard_name = name
        
            # Create landmask
            tempvar = ds.createVariable('land_binary_mask', float, dims, zlib=True, fill_value=1e20)
            tempvar[:] = tmask
            tempvar.standard_name = 'land_binary_mask'<jupyter_output><empty_output><jupyter_text>Process HRDPS<jupyter_code>def process_HRDPS(date, paths):
    """
    """

    datestr = date.strftime('ops_y%Ym%md%d.nc')
    HRDPS_out = os.path.join(paths['out'], 'forcing', datestr)
    if os.path.exists(HRDPS_out):
        print(f'{HRDPS_out} already exists! Skipping ...')
        return
    else:
        print(f'Processing {HRDPS_out}')
        with xr.open_dataset(os.path.join(paths['HRDPS'], datestr)) as data:
            time = data.time_counter
            time.attrs['standard_name'] = 'time'
            xr.Dataset(
                {
                    'longitude': ('flat', data.nav_lon.values.reshape(-1)),
                    'latitude': ('flat', data.nav_lat.values.reshape(-1)),
                    'u_wind': (['time_counter', 'flat'], data.u_wind.values.reshape(time.size, -1), {'standard_name': 'x_wind'}),
                    'v_wind': (['time_counter', 'flat'], data.v_wind.values.reshape(time.size, -1), {'standard_name': 'y_wind'}),
                },
                coords={'time_counter': time},
            ).to_netcdf(HRDPS_out)<jupyter_output><empty_output><jupyter_text>Loop processing<jupyter_code># Date range
process_mask(paths, mask)
for day in tqdm(range(np.diff(daterange)[0].days + 1)):
    date = daterange[0] + timedelta(days=day)
    process_NEMO(date, paths)
    #process_HRDPS(date, paths)
ds = xr.open_dataset('/data/bmoorema/results/opendrift/forcing/meshmask_opendrift.nc')
ds
ds = xr.open_dataset('/data/bmoorema/results/opendrift/forcing/SalishSea_1h_20171110_20171110_opendrift.nc')
ds<jupyter_output><empty_output><jupyter_text>***<jupyter_code>o = OceanDrift(loglevel=0)
datestr = daterange[0].strftime('%Y%m')
filenames = sorted(glob(os.path.join(paths['out'], 'forcing', f'SalishSea_1h_{datestr}*')))
#o.add_readers_from_list(filenames, lazy=False)
reader_landmask = reader_netCDF_CF_generic.Reader('/data/bmoorema/results/opendrift/forcing/meshmask_opendrift.nc')
reader_velocity = reader_netCDF_CF_generic.Reader(filenames[0])
o.add_reader([reader_landmask, reader_velocity])
o.seed_elements(-123.6, 49, radius=100, time=daterange[0], number=100)
o.set_config('general:coastline_action', 'previous')
o.set_config('drift:scheme', 'runge-kutta4')
o.set_config('general:use_basemap_landmask', False)
o.run(time_step=90, time_step_output=3600, duration=timedelta(days=1))
o
anim = o.animation()
HTML(anim.to_html5_video())<jupyter_output>14:18:40 INFO: Animation.save using <class 'matplotlib.animation.FFMpegWriter'>
14:18:40 INFO: figure size (inches) has been adjusted from 9.966152505214811 x 10.0 to 9.944444444444445 x 10.0
14:18:40 INFO: MovieWriter.run: running command: ['ffmpeg', '-f', 'rawvideo', '-vcodec', 'rawvideo', '-s', '716x720', '-pix_fmt', 'rgba', '-r', '20.0', '-loglevel', 'quiet', '-i', 'pipe:', '-vcodec', 'h264', '-pix_fmt', 'yuv420p', '-y', '/tmp/tmpp6djj609/temp.m4v']
{'fraction_non_alphanumeric': 0.14239038785834737, 'fraction_numerical': 0.02824620573355818, 'mean_word_length': 3.2723998199009454, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 34, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '41651', 'license': 'Apache-2.0', 'n_tokens_mistral': 3474, 'n_tokens_neox': 3277, 'n_words': 632, 'path': 'notebooks/OpenDrift/opendrift_forcing.ipynb', 'repo_name': 'SalishSeaCast/analysis-ben', 'size': '9488'}