Problem with changing the land mask

[zdulou]

Hello, I work in the Port of Bordeaux in numerical modelisation. I simulate the transport of plastic with openDrift in the Garonne and the Gironde estuary.

I run the simulation with a selafin file (.slf) made by open TELEMAC wich give me all the data of the water. The simulation are good on the ocean part.

I have a problem. The default land mask of opendrift doesn't include the river like Garonne and Dordogne. I want to change the land mask.

I have 2 options :
Part I : Use the land mask from the slf file

I want to use the land mask of the slf file wich include the river but the variable land_binary_mask doesn't exist in this file.

I change the land mask with this command :

o = OceanDrift(loglevel=0)
selafin = reader_telemac_selafin.Reader(filename=filename, proj4=proj, start_time=start_time)
o.add_reader(selafin, variables='land_binary_mask')
\\] o.add_reader(selafin)
o.set_config('general:use_auto_landmask', False)

and I have this error

error_part1.txt

The file doesn't have a variable call land_binary_mask, when I do print(selafin) it also doesn't appear in the variables.

I don't know how to modify this file.

Part II : Use the land mask from a shapeFile (.shp)
I have a python code that give me a shp file by using the slf file. So I have the same geometry with a shp file. Now I want to use the land mask of this shp file. I create another reader with this command :
shpfilename = 'gxl.shp'
reader_natural = reader_shape.Reader.from_shpfiles(shpfilename)
\\o.add_reader([reader_natural, selafin])
o.add_reader([reader_natural], variables='land_binary_mask')
o.set_config('general:use_auto_landmask', False) # Disabling the automatic GSHHG landmask
o.set_config('general:coastline_action', 'stranding')

Now the code run but I have a problem with the SRC the coordinate system.

I have this error :

error_part2.txt

It seems like OpenDrift mix 2 coordinate system the Lambert 93 and the WGS 84. The last ligne of the error is interesting : "Argcheck: all 10 particles (-1.14--0.34E, 45.11-45.51N) are outside domain of shape (258379.17-469018.50E, 6389451.50-6610523.50N)". The domain shape is inacurrate. There the value of the Lambert 93 system with the Wgs 84 coordinate system.

the selafin reader have an argument call proj4, how does it work ? I tried a lot of things. By default I have : proj = "+proj=lcc +lat_1=49.50000000000001 +lat_0=49.50000000000001 +lon_0=0
+k_0=0.999877341 +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356515
+units=m +no_defs"

I can have my shapefile in Lambert93 or in WGS 84

It's my first time posting here, tell me if I did something wrong or if you need more informations. I red a lot of forum, my code is a mix of a lot of things I found on internet and OpenDrift documentation.

Thanks

[knutfrode]

Hi,

I have not been involved in development of the Selafin-reader, but it seems that variable names there are hardcoded.
https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/reader_telemac_selafin.py

Thus you must probably add the missing variable (land_binary_mask), or update the reader to detect variables with standard_name attribute.
This could be done as in the generic unstructured reader (and ideally a specific Selafin-reader should not be necessary at all):
https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/reader_netCDF_CF_unstructured.py#L159

For the shapereader, it seems that you have to provide a string proj4_str upon creation, probably you can use the the one you have.
WGS84 is there default, so I would expect it to work if the shapefile used such coordinates already.