More clarity on landmask issue

[brucelaughlin]

Apologies, I was unclear in my previous post. I will try to explain more clearly.

Glorys has a few ocean points that are masked by the GSHHG land mask, and the GSHHG land mask has some ocean points which are land points in Glorys. For example, San Francisco bay is finely resolved in the GSHHG land mask, but, in Glorys, the bay is modeled by 10-20 grid cells. So, using the GSHHG land mask, Opendrift can send drifters into parts of the SF Bay which are land points in Glorys. I want to avoid this.

As a related example, this is from an Opendrift run using ROMS data in which the SF Bay didn't exist (ie was completely masked as land). Since I forgot to turn off the GSHHG land mask, Opendrift let particles drift into a masked part of the forcing data (the bay), apparently taking velocity values from the nearest non-masked points ( otside of the bay)

Disabling the GSHHG land mask solved the problem in this case, but I don't have mask variables in my Glorys data, so disabling the GSHHG land mask with Glorys data leads to "ValueError: Readers must be added for the following required variables: ['land_binary_mask']"

I hope this explanation of the problem is more clear. I'm wondering if there's a way to make Opendrift respect the masking (ie difference between land and ocean points) in forcing data, when there isn't a mask variable stored in the forcing data.

Thanks,
Bruce

[knutfrode]

Ok, yes this was clarifying.
OpenDrift will extrapolate current (and other) fields towards (vector) coastline to prevent particles from being trapped whenever model has no current, but still outside of actual landmask, which is independent and not perfectly consistent.

Normally this is ok/desired, but I agree in your case that you might not want too much extrapolation, and thus better to use ocean model landmask.
Since this seem not to be present in your dataset, one can e.g. define land whenever current is NaN. There is a built-in mechanism ("environment mapping") to calculate derived variables from others (e.g. land_binary_mask = 1 whenever current is infinite, and 0 elsewhere):

https://opendrift.github.io/_modules/opendrift/readers/basereader/variables.html#Variables
https://opendrift.github.io/_modules/opendrift/readers/basereader/variables.html#land_binary_mask_from_ocean_depth

However, this is not well documented, and thus you might instead want to use the workaround described in this new example:

https://opendrift.github.io/gallery/example_custom_xarray.html#sphx-glr-gallery-example-custom-xarray-py

[brucelaughlin]

Thank you so much, what you're saying makes sense, and I've gotten the example to work!

However, there seems to be another issue, UNRELATED TO THE MASK(!!), but related to usingxr.open_mfdataset vsxr.open_dataset, which I tried in order to pass multiple files to the reader (I just used these three lines, and the error appears):

his_files = "/home/blaughli/sample_files/test_forcing/aa_Mercator/onshore_uo/onshore*.nc"
ds = xr.open_mfdataset(his_files)
reader_current_job = Reader(ds)

This fails with the error output posted below (first problem is with my line reader_current_job = Reader(ds) and the final error message is RecursionError: maximum recursion depth exceeded). Note that there's only one file in my his_files directory - I was just testing with the wildcard.

[brucelaughlin]

Ok, since I couldn't get this approach to work with multiple input files, I've tried another strategy: adding a land binary mask and depths to my forcing file. Here's the output from ncks -m:

[brucelaughlin]

But, these two new variables in my netcdf file, land_binary_mask and sea_floor_depth_below_sea_level, are not being processed when I use the generic netcdf reader:

14:54:18 DEBUG opendrift.readers.reader_netCDF_CF_generic:368: Skipped variables without standard_name: ['land_binary_mask', 'sea_floor_depth_below_sea_level']

I'm also confused by this line:
14:54:18 DEBUG opendrift.readers.basereader:186: Variable mapping: ['sea_floor_depth_below_sea_level'] -> ['land_binary_mask'] is not activated
Is there an option to use the depth field as a land mask? That seems reasonable, though I shouldn't need it if I can get my new forcing file to work.

Please advise, thank you!

[brucelaughlin]

Ok, it's working now! I had to add standard names to the attributes of the variables I added to my model forcing. Now the generic netcdf reader seems to know about my added fields:

land_binary_mask = ncfile.createVariable("land_binary_mask","b",("latitude","longitude",))
land_binary_mask[:] = mask_rho
land_binary_mask.standard_name = "land_binary_mask"

sea_floor_depth_below_sea_level = ncfile.createVariable("sea_floor_depth_below_sea_level","f",("latitude","longitude",))
sea_floor_depth_below_sea_level[:] = model_grid_depths
sea_floor_depth_below_sea_level.standard_name = "sea_floor_depth_below_sea_level"