search

OpenDrift / opendrift

Open source framework for ocean trajectory modelling
https://opendrift.github.io
GNU General Public License v2.0
250 stars 120 forks source link

Telemac + own Shapefile + PlastDrift #1439

Open adegroodt opened 2 weeks ago

adegroodt commented 2 weeks ago

Hello OpenDrift community!

I am having issues trying to use the OpenTelemac output files and a shapefile created for a specific estuary. I've tried: using a global landmask - works fine, however the resolution is way too low for my area which is an estuary and includes sandbanks... So I want to add my own shapefile with the variables from opentelemac. Let me break it into what I have and what problems arise...

So I am using the gallery example of Telemac to open the Selafin file (https://opendrift.github.io/gallery/example_long_selafin.html), after importing all the required packages and setting log=0, this is how the code starts:

1. Calling filename path, defining projection (North Portugal), and adding it as reader

proj ="+proj=utm +zone=29 +north +datum=WGS84 +units=m +no_defs"     ##north Portugal -  Minho estuary
start_time= datetime(2013,1,1,0,0)
reader_selafin = reader_telemac_selafin.Reader(filename=filename, proj4 = proj, start_time=start_time)

output=

Reader: /Volumes/ExtremeSSD/BlueWWater/r2d_test_minho_11_dias.slf
Projection: 
  +proj=utm +zone=29 +north +datum=WGS84 +units=m
Coverage: [degrees]
  xmin: -71185.218750   xmax: -34285.222656
  ymin: 236247.703125   ymax: 265617.718750
  Corners (lon, lat):
    (-14.13,   2.39)  (-13.80,   2.39)
    (-14.13,   2.13)  (-13.80,   2.13)
Vertical levels [m]: 
  Not specified
Available time range:
  start: 2013-01-01 00:00:00   end: 2013-01-11 23:00:12   step: None
Variables:
  x_sea_water_velocity
  y_sea_water_velocity
===========================
 0:00:42.8  open dataset
 0:00:00.0  build index

--> Issue A) : the coordinates in the file is x/y, however the output of the reader is [degrees]. I've tried changing the projection to +proj=merc, it still stays in degrees, while I think the coverage should be in [m]?

--> Issue B): the variables x and y velocities are well read but not the other one such as water depth, sea surface... So I guess the issue is that in the source code of reader telemac it does not recognise the variable of the selafin file because they are not called the way opendrift would recognise them? Is it a solution to change the names of the variables?

2. Openin shapefile and adding it as reader

SHAPEFILE

from opendrift.readers import reader_shape

fn_shp ='./minho_MSLshp/min_margin.shp'
reader_shp = reader_shape.Reader.from_shpfiles(fn_shp, proj4_str=proj, invert=True)

Output=

Reader: shape Projection: +proj=utm +zone=29 +north +datum=WGS84 +units=m Coverage: [degrees] xmin: inf xmax: inf ymin: inf ymax: inf Corners (lon, lat): ( inf, inf) ( inf, inf) ( inf, inf) ( inf, inf) Vertical levels [m]: None Available time range: start: None end: None step: None Variables: land_binary_mask

--> Issue C): similar to issue A

Adding the readers and disable auto landmask

o.add_reader([reader_shp, reader_selafin])
o.set_config('general:use_auto_landmask', False)

Convert x and y to Lat/Lon

import pyproj
datum73_proj = pyproj.Proj(
    proj='utm',
    zone=29,
    ellps='GRS80',  # Using GRS80 as a close approximation, as DATUM73 is not directly available
    units='m',
    no_defs=True
)
latitude = 44.235,
longitude = -8.870893

# Convert to UTM (X: easting, Y: northing)
x, y = datum73_proj(longitude, latitude)

print(f"Easting (X): {x}, Northing (Y): {y}")
p= Proj(proj, preserve_units=False)
lon, lat = p(x,y,inverse=True)

z = -np.random.rand(200)*50

Seeding elements

o.seed_elements(lon=lon, lat=lat, z=z, radius=0, number= 200, time= start_time)
print(o)
Model:  PlastDrift     (OpenDrift version 1.10.7)
    0 active PlastElement particles  (0 deactivated, 200 scheduled)
-------------------
Environment variables:
  -----
  land_binary_mask
     1) shape
  -----
  x_sea_water_velocity
  y_sea_water_velocity
     1) /Volumes/ExtremeSSD/BlueWWater/r2d_test_minho_11_dias.slf
  -----
Readers not added for the following variables:
  ocean_mixed_layer_thickness
  ocean_vertical_diffusivity
  sea_floor_depth_below_sea_level
  sea_surface_height
  sea_surface_wave_significant_height
  sea_surface_wave_stokes_drift_x_velocity
  sea_surface_wave_stokes_drift_y_velocity
  x_wind
  y_wind

o.run(time_step=time_step/10, duration=timedelta(hours=24))

16:09:27 DEBUG   opendrift.models.basemodel:2400: 

Software and hardware:
  OpenDrift version 1.10.7
  Platform: Darwin, 23.0.0
  8.0 GB memory
  8 processors (arm)
  NumPy version 1.24.2
  SciPy version 1.11.3
  Matplotlib version 3.7.1
  NetCDF4 version 1.6.4
  Xarray version 2023.3.0
  Python version 3.10.10 (v3.10.10:aad5f6a891, Feb  7 2023, 08:47:40) [Clang 13.0.0 (clang-1300.0.29.30)]

16:09:27 DEBUG   opendrift.models.basemodel:2417: No output file is specified, neglecting export_buffer_length
16:09:27 INFO    opendrift.models.basemodel:2446: Fallback values will be used for the following variables which have no readers: 
16:09:27 INFO    opendrift.models.basemodel:2449:   sea_surface_height: 0.000000
16:09:27 INFO    opendrift.models.basemodel:2449:   sea_surface_wave_stokes_drift_x_velocity: 0.000000
16:09:27 INFO    opendrift.models.basemodel:2449:   sea_surface_wave_stokes_drift_y_velocity: 0.000000
16:09:27 INFO    opendrift.models.basemodel:2449:   sea_surface_wave_significant_height: 0.000000
16:09:27 INFO    opendrift.models.basemodel:2449:   x_wind: 0.000000
16:09:27 INFO    opendrift.models.basemodel:2449:   y_wind: 0.000000
16:09:27 INFO    opendrift.models.basemodel:2449:   ocean_vertical_diffusivity: 0.020000
16:09:27 INFO    opendrift.models.basemodel:2449:   ocean_mixed_layer_thickness: 50.000000
16:09:27 INFO    opendrift.models.basemodel:2449:   sea_floor_depth_below_sea_level: 10000.000000
16:09:27 DEBUG   opendrift.models.basemodel:2573: Preparing readers for simulation coverage ([-11.043663040142619, 42.67824385359481, -6.698122009295858, 45.79175736710832]) and time (2013-01-01 00:00:00 to 2013-01-02 00:00:00)
16:09:27 DEBUG   opendrift.models.basemodel:2577:   Preparing shape
16:09:27 DEBUG   opendrift.readers.basereader:188: Nothing more to prepare for shape
16:09:27 DEBUG   opendrift.models.basemodel:2577:   Preparing /Volumes/ExtremeSSD/BlueWWater/r2d_test_minho_11_dias.slf
16:09:27 DEBUG   opendrift.readers.basereader:188: Nothing more to prepare for /Volumes/ExtremeSSD/BlueWWater/r2d_test_minho_11_dias.slf
16:09:27 INFO    opendrift.models.basemodel:2608: Adding a dynamical landmask with max. priority based on assumed maximum speed of 2 m/s. Adding a customised landmask may be faster...

*** 16:09:27 DEBUG   opendrift.readers.basereader:176: Variable mapping: ['sea_floor_depth_below_sea_level'] -> ['land_binary_mask'] is not activated**
*** 16:09:32 DEBUG   opendrift.models.basemodel:786: Added reader global_landmask
*** 16:09:32 INFO    opendrift.models.basemodel:1556: Using existing reader for land_binary_mask
*** 16:09:32 DEBUG   opendrift.models.basemodel:990: Discarding reader (too far north): shape
*** 16:09:32 DEBUG   opendrift.models.basemodel:990: Discarding reader (too far south): /Volumes/ExtremeSSD/BlueWWater/r2d_test_minho_11_dias.slf

--> Issue D): Amongst the other errors that come by trying to change some parameters - the ones with the start *** are the ones that keep coming back and I have no idea why one is "too south" and the other "too north". I've tried changing the latitude/longitude but it still does no change anything....

Help would be greatly appreciated !! Thank you :))