issues running regridded ROMS

[jwongala]

Hello!

I needed to regrid a ROMS (Kaneohe) from a curvilinear to a rectilinear grid to nest it into a larger ROMS of the Hawaiian islands. I realized a few days ago that when I released particles from the area inside of the regridded ROMS, particles aren't moving. Right now, I am just trying to see if I can use only the regridded ROMS to advect particles but I am having trouble assigning the fieldset to the pset (attached error screenshot). The ROMS was regridded using pyferret. I also attached a screenshot of the ROMS info and pasted below the parcels code I am using. Any help or advice is appreciated! Thank you!

# Test for only KANEOHE ROMS RUNS

from parcels import FieldSet, Field, ParticleSet, JITParticle, ParticleFile, plotTrajectoriesFile
from parcels import AdvectionRK4
import numpy as np
from datetime import timedelta as delta
%matplotlib inline
import xarray as xr

file1_u = 'kbay_on_rect.nc'

file1_v = 'kbay_on_rect.nc'

u1 = xr.open_mfdataset(file1_u)
u2 = u1.rename({'OUT_U': 'u','XAX1': 'longitude','YAX': 'latitude', 'DEPTH1_1':'depth', 'AX009':'time'})

v1 = xr.open_mfdataset(file1_v)
v2 = v1.rename({'OUT_V': 'v','XAX1': 'longitude','YAX': 'latitude', 'DEPTH1_1':'depth', 'AX009':'time'})

dimensions = {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'}

U=Field.from_xarray(u2.u, 'u', dimensions, allow_time_extrapolation=True)
V=Field.from_xarray(v2.v, 'v', dimensions, allow_time_extrapolation=True)

fieldset = FieldSet(U,V)

pset = ParticleSet.from_list(fieldset=fieldset, pclass=JITParticle, lon=-157.7806, lat=21.44058)

[erikvansebille]

Thanks for reporting, @jwongala. What do your u2 and v2 look like (which attributes do they have)? And your fieldset?

[erikvansebille]

Dear @jwongala; do you still experience this Issue? I am closing it for now but feel free to reopen if you have more information

[jwongala]

Hi @erikvansebille ! I have come back to running these simulations and am having the same issue. Going back to your questions, this is what u2 and v2 look like (attached). When I try to look at the fieldset all I get as output is <parcels.fieldset.FieldSet at 0x1a9511090>. Is there another way I can view the fieldset output?

[erikvansebille]

Is it because you use lower-case u and v in your Field creation? What if you change to

U=Field.from_xarray(u2.u, 'U', dimensions, allow_time_extrapolation=True)
V=Field.from_xarray(v2.v, 'V', dimensions, allow_time_extrapolation=True)

[jwongala]

Thank you! The model works now.

Best, Jenn

[jwongala]

Hi @erikvansebille I was able to get the fine res ROMS (Kaneohe) running on its own, but now that I tried to nest it in the coarse ROMS (MHI), particles are only advected in the domain of the MHI and only released in the Kaneohe ROMS. Attached is the code for the model run. Please let me know if I should send it another way. I tried to paste it in but the formatting got hectic.

The only thing that I think might be an issue is that the Kaneohe ROMS is regridded so that it now has a rectilinear grid instead of a curvilinear grid. Thank you for your help!

Ch2_KANEOHE_NESTED_Bounce_TOOTS_Parcels.pdf

[jwongala]

Below is the code as well

import numpy as np import numpy.ma as ma from netCDF4 import Dataset import xarray as xr import pandas as pd from scipy import interpolate

from parcels import FieldSet, ParticleSet, JITParticle, ScipyParticle, AdvectionRK4, DiffusionUniformKh, Variable, Field, GeographicPolar, Geographic, ErrorCode, plotTrajectoriesFile, NestedField from datetime import timedelta as delta

file_name="test_toots_KANEOHE_NESTED_2019.zarr"

path2 = '/Users/wongalaj/TOOTS/TOOTS_parcels/test_toots_KANEOHE_NESTED_2019.zarr' # give path to where .nc file is

feather_path = 'test_toots_KANEOHE_NESTED_2019.feather' # assign name of feather file

startDate = '2019-01-01' # YYYY-MM-DD endDate = '2019-01-03' # YYYY-MM-DD

run_days = 4 # day 2018: 240, 2019: 365, 2020: 366, 2021: 365

simDepthMHI = 1

simDepthKANE = 1

kh = 10 # This is the eddy diffusivity in m2/s

pld = 3 # in days

file1_u = '/Users/wongalaj/TOOTS/TOOTS_parcels/ROMS_KANEOHE_Regrid/2019/u_vel/*.nc'

file1_v = '/Users/wongalaj/TOOTS/TOOTS_parcels/ROMS_KANEOHE_Regrid/2019/v_vel/*.nc'

u1 = xr.open_mfdataset(file1_u) u2 = u1.rename({'OUT_U_TIME': 'u','XAX': 'longitude','YAX': 'latitude', 'DEPTH':'depth', 'TIME':'time'})

v1 = xr.open_mfdataset(file1_v) v2 = v1.rename({'OUT_V_TIME': 'v','XAX': 'longitude','YAX': 'latitude', 'DEPTH':'depth', 'TIME':'time'})

dimensions = {'lon': 'longitude', 'lat': 'latitude', 'depth': 'depth', 'time': 'time'}

U1=Field.from_xarray(u2.u, 'U', dimensions, allow_time_extrapolation=True) # , interp_method = {'u': 'freeslip'} V1=Field.from_xarray(v2.v, 'V', dimensions, allow_time_extrapolation=True) # , interp_method = {'v': 'freeslip'}

Ufineres = U1 Vfineres = V1

file2 = 'https://pae-paha.pacioos.hawaii.edu/erddap/griddap/roms_hiig' # MHI ROMS (BASE)

ds2 = xr.open_dataset(file2) # this puts the opendap data into a xarray dataset

myDat2 = ds2.isel({'depth': simDepthMHI}).sel({'time': slice(startDate,endDate)}) # subset based on time and depth layer

variables = {'U': 'u', 'V': 'v'} dimensions = {'lon': 'longitude', 'lat': 'latitude', 'time': 'time'}

U=Field.from_xarray(myDat2.u, 'U', dimensions, allow_time_extrapolation=True) V=Field.from_xarray(myDat2.v, 'V', dimensions, allow_time_extrapolation=True)

Ucoarse = U Vcoarse = V

Ucoarse.show()

UNested = NestedField('U', [Ufineres, Ucoarse]) VNested = NestedField('V', [Vfineres, Vcoarse])

U = UNested V = VNested

fieldset = FieldSet(U, V)

def make_landmask1(fielddata): """Returns landmask where land = 1 and ocean = 0 fielddata is a netcdf file. """ datafile = Dataset(fielddata)

landmask = datafile.variables['OUT_U_TIME'][0, 0]
landmask = np.ma.masked_invalid(landmask)
landmask = landmask.mask.astype('int')

return landmask

def make_landmask2(fielddata): """Returns landmask where land = 1 and ocean = 0 fielddata is a netcdf file. """ datafile = Dataset(fielddata)

landmask = datafile.variables['u'][0, 0]
landmask = np.ma.masked_invalid(landmask)
landmask = landmask.mask.astype('int')

return landmask

file_path1 = '/Users/wongalaj/TOOTS/TOOTS_parcels/ROMS_KANEOHE_Regrid/2018/u_vel/out_u_timerect_2018_01.nc'

file_path2 = '/Users/wongalaj/TOOTS/TOOTS_parcels/ROMS_MHI/roms_hiig_2018_05_05_685a_c38d_b3fe.nc'

landmask_fineres = make_landmask1(file_path1) # KANEOHE

landmask_coarse = make_landmask2(file_path2) # MHI

def get_coastal_nodes(landmask): """Function that detects the coastal nodes, i.e. the ocean nodes directly next to land. Computes the Laplacian of landmask.

- landmask: the land mask built using `make_landmask`, where land cell = 1
            and ocean cell = 0.

Output: 2D array array containing the coastal nodes, the coastal nodes are
        equal to one, and the rest is zero.
"""
mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
mask_lap -= 4*landmask
coastal = np.ma.masked_array(landmask, mask_lap > 0)
coastal = coastal.mask.astype('int')

return coastal

def get_shore_nodes(landmask): """Function that detects the shore nodes, i.e. the land nodes directly next to the ocean. Computes the Laplacian of landmask.

- landmask: the land mask built using `make_landmask`, where land cell = 1
            and ocean cell = 0.

Output: 2D array array containing the shore nodes, the shore nodes are
        equal to one, and the rest is zero.
"""
mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
mask_lap -= 4*landmask
shore = np.ma.masked_array(landmask, mask_lap < 0)
shore = shore.mask.astype('int')

return shore

def get_coastal_nodes_diagonal(landmask): """Function that detects the coastal nodes, i.e. the ocean nodes where one of the 8 nearest nodes is land. Computes the Laplacian of landmask and the Laplacian of the 45 degree rotated landmask.

- landmask: the land mask built using `make_landmask`, where land cell = 1
            and ocean cell = 0.

Output: 2D array array containing the coastal nodes, the coastal nodes are
        equal to one, and the rest is zero.
"""
mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
mask_lap += np.roll(landmask, (-1,1), axis=(0,1)) + np.roll(landmask, (1, 1), axis=(0,1))
mask_lap += np.roll(landmask, (-1,-1), axis=(0,1)) + np.roll(landmask, (1, -1), axis=(0,1))
mask_lap -= 8*landmask
coastal = np.ma.masked_array(landmask, mask_lap > 0)
coastal = coastal.mask.astype('int')

return coastal

def get_shore_nodes_diagonal(landmask): """Function that detects the shore nodes, i.e. the land nodes where one of the 8 nearest nodes is ocean. Computes the Laplacian of landmask and the Laplacian of the 45 degree rotated landmask.

- landmask: the land mask built using `make_landmask`, where land cell = 1
            and ocean cell = 0.

Output: 2D array array containing the shore nodes, the shore nodes are
        equal to one, and the rest is zero.
"""
mask_lap = np.roll(landmask, -1, axis=0) + np.roll(landmask, 1, axis=0)
mask_lap += np.roll(landmask, -1, axis=1) + np.roll(landmask, 1, axis=1)
mask_lap += np.roll(landmask, (-1,1), axis=(0,1)) + np.roll(landmask, (1, 1), axis=(0,1))
mask_lap += np.roll(landmask, (-1,-1), axis=(0,1)) + np.roll(landmask, (1, -1), axis=(0,1))
mask_lap -= 8*landmask
shore = np.ma.masked_array(landmask, mask_lap < 0)
shore = shore.mask.astype('int')

return shore

coastal_fineres = get_coastal_nodes_diagonal(landmask_fineres) shore_fineres = get_shore_nodes_diagonal(landmask_fineres)

coastal_coarse = get_coastal_nodes_diagonal(landmask_coarse) shore_coarse = get_shore_nodes_diagonal(landmask_coarse)

def create_displacement_field(landmask, double_cell=False): """Function that creates a displacement field 1 m/s away from the shore.

- landmask: the land mask dUilt using `make_landmask`.
- double_cell: Boolean for determining if you want a double cell.
  Default set to False.

Output: two 2D arrays, one for each camponent of the velocity.
"""
shore = get_shore_nodes(landmask)
shore_d = get_shore_nodes_diagonal(landmask) # bordering ocean directly and diagonally
shore_c = shore_d - shore                    # corner nodes that only border ocean diagonally

Ly = np.roll(landmask, -1, axis=0) - np.roll(landmask, 1, axis=0) # Simple derivative
Lx = np.roll(landmask, -1, axis=1) - np.roll(landmask, 1, axis=1)

Ly_c = np.roll(landmask, -1, axis=0) - np.roll(landmask, 1, axis=0)
Ly_c += np.roll(landmask, (-1,-1), axis=(0,1)) + np.roll(landmask, (-1,1), axis=(0,1)) # Include y-component of diagonal neighbours
Ly_c += - np.roll(landmask, (1,-1), axis=(0,1)) - np.roll(landmask, (1,1), axis=(0,1))

Lx_c = np.roll(landmask, -1, axis=1) - np.roll(landmask, 1, axis=1)
Lx_c += np.roll(landmask, (-1,-1), axis=(1,0)) + np.roll(landmask, (-1,1), axis=(1,0)) # Include x-component of diagonal neighbours
Lx_c += - np.roll(landmask, (1,-1), axis=(1,0)) - np.roll(landmask, (1,1), axis=(1,0))

v_x = -Lx*(shore)
v_y = -Ly*(shore)

v_x_c = -Lx_c*(shore_c)
v_y_c = -Ly_c*(shore_c)

v_x = v_x + v_x_c
v_y = v_y + v_y_c

magnitude = np.sqrt(v_y**2 + v_x**2)
# the coastal nodes between land create a problem. Magnitude there is zero
# I force it to be 1 to avoid problems when normalizing.
ny, nx = np.where(magnitude == 0)
magnitude[ny, nx] = 1

v_x = v_x/magnitude
v_y = v_y/magnitude

return v_x, v_y

v_x_f, v_y_f = create_displacement_field(landmask_fineres) v_x_c, v_y_c = create_displacement_field(landmask_coarse)

def distance_to_shore(landmask, dx=1): """Function that computes the distance to the shore. It is based in the the get_coastal_nodes algorithm.

- landmask: the land mask dUilt using `make_landmask` function.
- dx: the grid cell dimension. This is a crude approxsimation of the real
distance (be careful).

Output: 2D array containing the distances from shore.
"""
ci = get_coastal_nodes(landmask) # direct neighbours
dist = ci*dx                     # 1 dx away

ci_d = get_coastal_nodes_diagonal(landmask) # diagonal neighbours
dist_d = (ci_d - ci)*np.sqrt(2*dx**2)       # sqrt(2) dx away

return dist+dist_d

d_2_s_f = distance_to_shore(landmask_fineres) d_2_s_c = distance_to_shore(landmask_coarse)

u_displacement_f = v_x_f v_displacement_f = v_y_f

u_displacement_c = v_x_c v_displacement_c = v_y_c

fieldset.add_field(Field('dispUF', data=u_displacement_f, lon=fieldset.U[0].grid.lon, lat=fieldset.U[0].grid.lat, mesh='spherical'))

fieldset.add_field(Field('dispVF', data=v_displacement_f, lon=fieldset.V[0].grid.lon, lat=fieldset.V[0].grid.lat, mesh='spherical'))

fieldset.add_field(Field('dispUC', data=u_displacement_c, lon=fieldset.U[1].grid.lon, lat=fieldset.U[1].grid.lat, mesh='spherical'))

fieldset.add_field(Field('dispVC', data=v_displacement_c, lon=fieldset.V[1].grid.lon, lat=fieldset.V[1].grid.lat, mesh='spherical'))

fieldset.dispUF.units = GeographicPolar() fieldset.dispUC.units = GeographicPolar() fieldset.dispVF.units = Geographic() fieldset.dispVC.units = Geographic()

fieldset.add_field(Field('landmask_fine', landmask_fineres, lon=fieldset.U[0].grid.lon, lat=fieldset.U[0].grid.lat, mesh='spherical'))

fieldset.add_field(Field('landmask_coarse', landmask_coarse, lon=fieldset.U[1].grid.lon, lat=fieldset.U[1].grid.lat, mesh='spherical'))

fieldset.add_field(Field('distance2shore_fine', d_2_s_f, lon=fieldset.U[0].grid.lon, lat=fieldset.U[0].grid.lat, mesh='spherical'))

fieldset.add_field(Field('distance2shore_coarse', d_2_s_c, lon=fieldset.U[1].grid.lon, lat=fieldset.U[1].grid.lat, mesh='spherical'))

fieldset.add_constant_field('Kh_zonal', kh, mesh='spherical') fieldset.add_constant_field('Kh_meridional', kh, mesh='spherical')

fieldset.add_constant('pld', (pld*86400))

class DisplacementParticle(JITParticle): dU = Variable('dU') dV = Variable('dV') d2s = Variable('d2s', initial=1e3) age = Variable('age', dtype=np.float32, initial=0.) releaseSite = Variable('releaseSite', dtype=np.int32) IslandReleaseSite = Variable('IslandReleaseSite', dtype=np.int32)

def set_displacement(particle, fieldset, time):

try saying in between max and min

if  particle.lon <= fieldset.max_lon and particle.lon >= fieldset.min_lon and particle.lat <= fieldset.max_lat and particle.lat >= fieldset.min_lat:   
    particle.d2s = fieldset.distance2shore_fine[time, particle.depth,
                           particle.lat, particle.lon]
    if  particle.d2s < 0.5:
        dispUab = fieldset.dispUF[time, particle.depth, particle.lat,
                           particle.lon]
        dispVab = fieldset.dispVF[time, particle.depth, particle.lat,
                           particle.lon]

        particle.dU = dispUab
        particle.dV = dispVab
    else:
        particle.dU = 0.
        particle.dV = 0.
else:
    particle.d2s = fieldset.distance2shore_coarse[time, particle.depth,
                           particle.lat, particle.lon]

    if  particle.d2s < 0.5:
        dispUab = fieldset.dispUC[time, particle.depth, particle.lat,
                           particle.lon]
        dispVab = fieldset.dispVC[time, particle.depth, particle.lat,
                           particle.lon]

        particle.dU = dispUab
        particle.dV = dispVab
    else:
        particle.dU = 0.
        particle.dV = 0.

def displace(particle, fieldset, time):
if particle.d2s < 0.5: particle.lon += particle.dUparticle.dt particle.lat += particle.dVparticle.dt

def DeleteParticle(particle, fieldset, time): print('deleted particle') particle.delete()

def Ageing(particle, fieldset, time): particle.age += particle.dt if particle.age >= fieldset.pld: particle.delete()

fieldset.add_constant('min_lon', -158.02)

fieldset.add_constant('max_lon', -157.6207)

fieldset.add_constant('min_lat', 21.34312)

fieldset.add_constant('max_lat', 21.71657)

fieldset.add_constant('min_lon', -163.83070)

fieldset.add_constant('max_lon', -152.51930)

fieldset.add_constant('min_lat', 17.01843)

fieldset.add_constant('max_lat', 23.98238)

nrepeat = 1 # how many times do you want locations to repeat (original = 500)

source_loc = pd.read_csv("all_release_locs_oocysts_final.csv") # read in file I make

toots_norm = source_loc.oocysts_normalized #assign column of oocysts_normalized to toots-norm

par_release_prop = nrepeat * toots_norm # multiple nrepeat to the toots_norm data to get the number of particles that need to be added in addition to the 500 minimum

release_par_final = par_release_prop + nrepeat # add the additional particles to nrepeat (500) so all locations release a minimum of 500 particles (max of 1000)

lon_fin = np.array([]) lat_fin = np.array([]) site_fin = np.array([]) islandsite_fin = np.array([])

for i in range(0,117):

tmp = release_par_final[i] # get out number of particles to be released at that specific location using [i]

lon1 = np.repeat(source_loc.lon[i], tmp) # subset out the variables I need to be repeated in the pset and replicate it by number of particles needed
lat1 = np.repeat(source_loc.lat[i], tmp)
site1 = np.repeat(source_loc.ID[i], tmp)
islandsite1 = np.repeat(source_loc.island_release[i], tmp)

lon_fin = np.append(lon_fin, lon1) # bind all of the initial conditions data to one final vector for the model psete input
lat_fin = np.append(lat_fin, lat1)
site_fin = np.append(site_fin, site1)
islandsite_fin = np.append(islandsite_fin, islandsite1)

habilon = lon_fin habilat = lat_fin habisite = site_fin islandsite = islandsite_fin

release_days = 1 # how often to release particles (original == 15 days)

release_int = release_days*86400 # 15 days converted to seconds

pset = ParticleSet.from_list(fieldset=fieldset, #parameter found in ParticleSet pclass=DisplacementParticle, #parameter found in ParticleSet lon=habilon, #parameter found in ParticleSet lat=habilat, #parameter found in ParticleSet releaseSite=habisite, IslandReleaseSite=islandsite, repeatdt=release_int) #parameter found in ParticleSet

kernels = pset.Kernel(displace) + pset.Kernel(AdvectionRK4) + pset.Kernel(DiffusionUniformKh) + pset.Kernel(set_displacement) + pset.Kernel(Ageing)

output_file = pset.ParticleFile(name=file_name, outputdt=delta(hours=4))

pset.execute(kernels, runtime=delta(days=run_days), dt=delta(minutes=15), recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}, output_file=output_file)

pset.repeatdt = None

pset.execute(kernels, runtime=delta(days=pld+1), dt=delta(minutes=15), recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}, output_file=output_file)

pset

plt = plotTrajectoriesFile('test_toots_KANEOHE_NESTED_2019.zarr')

[erikvansebille]

Hi @jwongala, I'm afraid I can't help you here. It will be very difficult for me to debug your code without access to your hydrodynamic files. I had a quick look at your code, but didn't notice anything odd or suspicious.

[jwongala]

Hello,Okay thank you very much! I will update this issues if anything new happens. I think I will just do the run using only the Kaneohe ROMS like I originally planned. On Dec 11, 2023, at 12:03 AM, Erik van Sebille @.***> wrote: Hi @jwongala, I'm afraid I can't help you here. It will be very difficult for me to debug your code without access to your hydrodynamic files. I had a quick look at your code, but didn't notice anything odd or suspicious.

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>