OceanParcels / Parcels

Main code for Parcels (Probably A Really Computationally Efficient Lagrangian Simulator)
https://www.oceanparcels.org
MIT License
294 stars 135 forks source link

NestedField from FieldSet.from_nemo #1429

Closed pablohs27 closed 3 months ago

pablohs27 commented 1 year ago

Hello I'm testing OceanParcels' NestedField function as I have two fields with different resolutions from a FieldSet.from_nemo field. Until now, when I run the code with the structure suggested in the function, it does not recognize the second field when loading the nestings in U and V, but it does not give me any errors either.

This is the structure I am following:

# Field 1
mfpathu = 'E:/MAESTRIA/DATOS/U/'
mfpathv = 'E:/MAESTRIA/DATOS/V/'
mfpathd = 'E:/MAESTRIA/DATOS/Domain/1_domain_cfg.nc'
mfufiles = (mfpathu + '1_CARIB_1h_20090101_20090131_grid_U.nc')
mfvfiles = (mfpathv + '1_CARIB_1h_20090101_20090131_grid_V.nc')

mffilenames = {'U': {'lon': mfpathd,'lat': mfpathd,'data': mfufiles},
             'V': {'lon': mfpathd,'lat': mfpathd,'data': mfvfiles}}
mfvariables = {'U': 'uo', 'V': 'vo'}
mfdimensions = {'lat': 'gphif','lon': 'glamf','time': 'time_counter'}

mffieldset = FieldSet.from_nemo(mffilenames, mfvariables, mfdimensions, allow_time_extrapolation=False)

# Field 2
mgpathu = 'E:/MAESTRIA/DATOS/MallaGrande/'
mgpathv = 'E:/MAESTRIA/DATOS/MallaGrande/'
mgpathd = 'E:/MAESTRIA/DATOS/MallaGrande/domain_cfg.nc'
mgufiles = (mgpathu + 'CARIB_1h_20090101_20090131_grid_U.nc')
mgvfiles = (mgpathv + 'CARIB_1h_20090101_20090131_grid_V.nc')

mgfilenames = {'U': {'lon': mgpathd,'lat': mgpathd,'data': mgufiles},
             'V': {'lon': mgpathd,'lat': mgpathd,'data': mgvfiles}}
mgvariables = {'U': 'uo', 'V': 'vo'}
mgdimensions = {'lat': 'gphif','lon': 'glamf','time': 'time_counter'}

mgfieldset = FieldSet.from_nemo(mgfilenames, mgvariables, mgdimensions, allow_time_extrapolation=False)

U = NestedField('U', [mffieldset.U, mgfieldset.U])
V = NestedField('V', [mffieldset.V, mgfieldset.V])

fieldset = FieldSet(U, V)

Why could it not read the second field? I've noticed that when I reverse the fields in NestedField for U and V, it still doesn't recognize the second input field. Greetings

erikvansebille commented 1 year ago

Dear @pablohs27, thanks for asking this question. And also good that you posted some code. But in order to help, I'd need to understand better what you mean with 'read the second field'. What behaviour do you see, and what did you expect? Please be as specific as possible, otherwise it'll be very difficult for me to find out what's going wrong

pablohs27 commented 1 year ago

Thank you Dr. Erick for your response. I'm specifically referring to the following part of the code:

U = NestedField('U', [mffieldset.U, mgfieldset.U])
V = NestedField('V', [mffieldset.V, mgfieldset.V])

When I load the velocities in U and V for both fields. "mffieldset" is my field 1 (with a higher resolution) and "mgfieldset" is my field 2.

The code works but when I plot the trajectories of the particles, only the trajectories of field 1 entered in the NestedField function appear. I expect my particles to go from my field 1 to my field 2. But field 2 is not recognized.

erikvansebille commented 1 year ago

Thanks for this extra information, but it's not really what I meant. To be able to help you, I would need a minimal reproducible example. So provide the full code, the full output etc.

pablohs27 commented 1 year ago

This the full code:

from parcels import FieldSet, ParticleSet, JITParticle, ParticleFile, plotTrajectoriesFile, ErrorCode, AdvectionRK4, Field, NestedField
from argparse import ArgumentParser
import numpy as np
import pytest
from glob import glob
from datetime import timedelta as delta 
from os import path
import time
from netCDF4 import Dataset
import xarray as xr
import cartopy
import matplotlib.pyplot as plt

# Field 1
mfpathu = 'E:/MAESTRIA/DATOS/U/'
mfpathv = 'E:/MAESTRIA/DATOS/V/'
mfpathd = 'E:/MAESTRIA/DATOS/Domain/1_domain_cfg.nc'
mfufiles = (mfpathu + '1_CARIB_1h_20090101_20090131_grid_U.nc')
mfvfiles = (mfpathv + '1_CARIB_1h_20090101_20090131_grid_V.nc')

mffilenames = {'U': {'lon': mfpathd,'lat': mfpathd,'data': mfufiles},
             'V': {'lon': mfpathd,'lat': mfpathd,'data': mfvfiles}}
mfvariables = {'U': 'uo', 'V': 'vo'}
mfdimensions = {'lat': 'gphif','lon': 'glamf','time': 'time_counter'}

mffieldset = FieldSet.from_nemo(mffilenames, mfvariables, mfdimensions)

# Field 2
mgpathu = 'E:/MAESTRIA/DATOS/MallaGrande/'
mgpathv = 'E:/MAESTRIA/DATOS/MallaGrande/'
mgpathd = 'E:/MAESTRIA/DATOS/MallaGrande/domain_cfg.nc'
mgufiles = (mgpathu + 'CARIB_1h_20090101_20090131_grid_U.nc')
mgvfiles = (mgpathv + 'CARIB_1h_20090101_20090131_grid_V.nc')

mgfilenames = {'U': {'lon': mgpathd,'lat': mgpathd,'data': mgufiles},
             'V': {'lon': mgpathd,'lat': mgpathd,'data': mgvfiles}}
mgvariables = {'U': 'uo', 'V': 'vo'}
mgdimensions = {'lat': 'gphif','lon': 'glamf','time': 'time_counter'}

mgfieldset = FieldSet.from_nemo(mgfilenames, mgvariables, mgdimensions)

U = NestedField('U', [mffieldset.U, mgfieldset.U])
V = NestedField('V', [mffieldset.V, mgfieldset.V])

fieldset = FieldSet(U, V)

pset = ParticleSet(fieldset, pclass=JITParticle, lon=[-93.33], lat=[18.51])

output_file = pset.ParticleFile(
    name="27NFTT.zarr", outputdt=7200)

pset.execute(AdvectionRK4, runtime=1728000, dt=600, output_file=output_file)

ds = xr.open_zarr("27NFTT.zarr")
plt.plot(ds.lon.T, ds.lat.T, ".-")
plt.show()

And this is the output when plot trajectories:

image

Although the code runs and doesn't give me any errors. The straight line shows the edge of field 1, so the trajectory does not pass into field 2

erikvansebille commented 1 year ago

OK, this helps a bit but I still don't know how to replicate this. I don't know what your velocity feels are like, where they come from, what their domains are, how the grids are defined etc.

I'm sorry, but you need to provide me with a version that I can run myself; otherwise I can't help.