search

OpenDrift / opendrift

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

'OceanDrift' object has no attribute 'environment' #593

Closed GuiCruz closed 3 years ago

GuiCruz commented 3 years ago

Hello, I am trying to implement my results obtained by delft3D model on opendrift. I had to convert the sigma coordinates to z-coordinate and the opendrift seems to have read it correctly. First I used this example for the conversion (https://rdrr.io/github/edwardlavender/fvcom.tbx/src/R/depth_from_known.R ) and then interpolated the data to a regular z-coordinate. However I am getting the following error:

reader = reader_netCDF_CF_generic.Reader(fn_nc)
reader
-------------------------------------------------------
06:46:29 INFO: Opening dataset: ./wind10_225deg_+2500.nc
06:46:29 INFO: Opening file with Dataset
06:46:29 DEBUG: Finding coordinate variables.
06:46:29 DEBUG: Parsing variable: latitude_longitude
06:46:29 DEBUG: Parsing CF grid mapping dictionary: {'grid_mapping_name': 'latitude_longitude'}
06:46:29 INFO: Could not parse CF grid_mapping
06:46:29 DEBUG: Parsing variable: depth
06:46:29 DEBUG: Parsing variable: time
06:46:29 DEBUG: Parsing variable: u
06:46:29 DEBUG: Parsing variable: v
06:46:29 DEBUG: Parsing variable: w
06:46:29 DEBUG: Parsing variable: d
06:46:29 DEBUG: Parsing variable: s
06:46:29 DEBUG: Parsing variable: lon
06:46:29 DEBUG: Parsing variable: lat
06:46:29 INFO: No projection found, using lon/lat arrays
06:46:29 DEBUG: Reading lon lat 2D arrays, since projection is not given
06:46:29 INFO: Making Splines for lon,lat to x,y conversion...
06:46:31 DEBUG: Setting buffer size 18 for reader ./wind10_225deg_+2500.nc, assuming a maximum average speed of 5 m/s.

===========================
Reader: ./wind10_225deg_+2500.nc
Projection: 
  None
Coverage: [pixels]
  xmin: 0.000000   xmax: 198.000000   step: 1   numx: 199
  ymin: 0.000000   ymax: 501.000000   step: 1   numy: 502
  Corners (lon, lat):
    (-52.18, -31.83)  (-51.97, -31.83)
    (-52.19, -32.28)  (-51.98, -32.28)
Vertical levels [m]: 
  [ -0.   -0.1  -0.2  -0.3  -0.4  -0.5  -0.6  -0.7  -0.8  -0.9  -1.   -1.1
  -1.2  -1.3  -1.4  -1.5  -1.6  -1.7  -1.8  -1.9  -2.   -2.1  -2.2  -2.3
  -2.4  -2.5  -2.6  -2.7  -2.8  -2.9  -3.   -3.1  -3.2  -3.3  -3.4  -3.5
  -3.6  -3.7  -3.8  -3.9  -4.   -4.1  -4.2  -4.3  -4.4  -4.5  -4.6  -4.7
  -4.8  -4.9  -5.   -5.1  -5.2  -5.3  -5.4  -5.5  -5.6  -5.7  -5.8  -5.9
  -6.   -6.1  -6.2  -6.3  -6.4  -6.5  -6.6  -6.7  -6.8  -6.9  -7.   -7.1
  -7.2  -7.3  -7.4  -7.5  -7.6  -7.7  -7.8  -7.9  -8.   -8.1  -8.2  -8.3
  -8.4  -8.5  -8.6  -8.7  -8.8  -8.9  -9.   -9.1  -9.2  -9.3  -9.4  -9.5
  -9.6  -9.7  -9.8  -9.9 -10.  -10.1 -10.2 -10.3 -10.4 -10.5 -10.6 -10.7
 -10.8 -10.9 -11.  -11.1 -11.2 -11.3 -11.4 -11.5 -11.6 -11.7 -11.8 -11.9
 -12.  -12.1 -12.2 -12.3 -12.4 -12.5 -12.6 -12.7 -12.8 -12.9 -13.  -13.1
 -13.2 -13.3 -13.4 -13.5 -13.6 -13.7 -13.8 -13.9 -14.  -14.1 -14.2 -14.3
 -14.4 -14.5 -14.6 -14.7 -14.8 -14.9 -15.  -15.1 -15.2 -15.3 -15.4 -15.5
 -15.6 -15.7 -15.8 -15.9 -16.  -16.1 -16.2 -16.3 -16.4 -16.5 -16.6 -16.7
 -16.8 -16.9 -17.  -17.1 -17.2 -17.3 -17.4 -17.5 -17.6 -17.7 -17.8 -17.9
 -18.  -18.1 -18.2 -18.3 -18.4 -18.5 -18.6 -18.7 -18.8 -18.9 -19.  -19.1
 -19.2 -19.3 -19.4 -19.5 -19.6 -19.7 -19.8 -19.9 -20. ]
Available time range:
  start: 2020-07-09 00:05:00   end: 2020-07-09 00:20:00   step: 0:05:00
    4 times (0 missing)
Variables:
  time
  y_sea_water_velocity
  x_sea_water_velocity
  upward_sea_water_velocity
  ocean_vertical_diffusivity
  sea_water_salinity
===========================
06:42:39 DEBUG: Adding 17 config items from basemodel
06:42:40 DEBUG: Adding 4 config items from basemodel
06:42:40 DEBUG: Adding 34 config items from basemodel
06:42:40 INFO: OpenDriftSimulation initialised (version 1.4.2)
06:42:40 DEBUG: Adding 13 config items from oceandrift
06:42:40 DEBUG:   Overwriting config item seed:z
06:42:40 INFO: Opening dataset: ./wind10_225deg_+2500.nc
06:42:40 INFO: Opening file with Dataset
06:42:40 DEBUG: Finding coordinate variables.
06:42:40 DEBUG: Parsing variable: latitude_longitude
06:42:40 DEBUG: Parsing CF grid mapping dictionary: {'grid_mapping_name': 'latitude_longitude'}
06:42:40 INFO: Could not parse CF grid_mapping
06:42:40 DEBUG: Parsing variable: depth
06:42:40 DEBUG: Parsing variable: time
06:42:40 DEBUG: Parsing variable: u
06:42:40 DEBUG: Parsing variable: v
06:42:40 DEBUG: Parsing variable: w
06:42:40 DEBUG: Parsing variable: d
06:42:40 DEBUG: Parsing variable: s
06:42:40 DEBUG: Parsing variable: lon
06:42:40 DEBUG: Parsing variable: lat
06:42:40 INFO: No projection found, using lon/lat arrays
06:42:40 DEBUG: Reading lon lat 2D arrays, since projection is not given
06:42:40 INFO: Making Splines for lon,lat to x,y conversion...
06:42:42 DEBUG: Setting buffer size 18 for reader ./wind10_225deg_+2500.nc, assuming a maximum average speed of 5 m/s.
06:42:42 DEBUG: Reading shapefile: D:/Dropbox/Doutorado/LC08_RS/Ideal_model/Coastline_25x100km_GCS.shp
06:42:42 INFO: Pre-processing 2 geometries
06:42:42 DEBUG: Setting buffer size 6 for reader ./wind10_225deg_+2500.nc, assuming a maximum average speed of 1 m/s.
06:42:42 DEBUG: ./wind10_225deg_+2500.nc is unprojected, cannot use for common grid
06:42:42 DEBUG: Added reader ./wind10_225deg_+2500.nc
06:42:42 INFO: Setting SRS to latlong, since not defined before.
06:42:42 DEBUG: Calculation SRS set to: +proj=longlat +ellps=WGS84 +no_defs
06:42:42 DEBUG: Added reader shape
06:42:42 INFO: Seeding for layer: seeder (95 features)
06:42:42 INFO: Total area of all polygons: 903608.4602117687 m2
06:42:42 INFO:  Seeding 129 elements within polygon number 1
06:42:42 INFO:  Seeding 100 elements within polygon number 2
06:42:42 INFO:  Seeding 100 elements within polygon number 3
06:42:42 INFO:  Seeding 100 elements within polygon number 4
06:42:42 INFO:  Seeding 100 elements within polygon number 5
06:42:42 INFO:  Seeding 100 elements within polygon number 6
06:42:42 INFO:  Seeding 100 elements within polygon number 7
06:42:42 INFO:  Seeding 100 elements within polygon number 8
06:42:42 INFO:  Seeding 100 elements within polygon number 9
06:42:42 INFO:  Seeding 100 elements within polygon number 10
06:42:42 INFO:  Seeding 100 elements within polygon number 11
06:42:42 INFO:  Seeding 100 elements within polygon number 12
06:42:42 INFO:  Seeding 100 elements within polygon number 13
06:42:42 INFO:  Seeding 100 elements within polygon number 14
06:42:42 INFO:  Seeding 100 elements within polygon number 15
06:42:42 INFO:  Seeding 100 elements within polygon number 16
06:42:42 INFO:  Seeding 100 elements within polygon number 17
06:42:42 INFO:  Seeding 100 elements within polygon number 18
06:42:42 INFO:  Seeding 100 elements within polygon number 19
06:42:42 INFO:  Seeding 100 elements within polygon number 20
06:42:42 INFO:  Seeding 100 elements within polygon number 21
06:42:42 INFO:  Seeding 100 elements within polygon number 22
06:42:42 INFO:  Seeding 100 elements within polygon number 23
06:42:42 INFO:  Seeding 71 elements within polygon number 24
06:42:42 INFO:  Seeding 100 elements within polygon number 25
06:42:42 INFO:  Seeding 100 elements within polygon number 26
06:42:42 INFO:  Seeding 100 elements within polygon number 27
06:42:43 INFO:  Seeding 100 elements within polygon number 28
06:42:43 INFO:  Seeding 100 elements within polygon number 29
06:42:43 INFO:  Seeding 100 elements within polygon number 30
06:42:43 INFO:  Seeding 100 elements within polygon number 31
06:42:43 INFO:  Seeding 100 elements within polygon number 32
06:42:43 INFO:  Seeding 100 elements within polygon number 33
06:42:43 INFO:  Seeding 100 elements within polygon number 34
06:42:43 INFO:  Seeding 100 elements within polygon number 35
06:42:43 INFO:  Seeding 100 elements within polygon number 36
06:42:43 INFO:  Seeding 100 elements within polygon number 37
06:42:43 INFO:  Seeding 100 elements within polygon number 38
06:42:43 INFO:  Seeding 100 elements within polygon number 39
06:42:43 INFO:  Seeding 100 elements within polygon number 40
06:42:43 INFO:  Seeding 100 elements within polygon number 41
06:42:43 INFO:  Seeding 100 elements within polygon number 42
06:42:43 INFO:  Seeding 100 elements within polygon number 43
06:42:43 INFO:  Seeding 100 elements within polygon number 44
06:42:43 INFO:  Seeding 100 elements within polygon number 45
06:42:43 INFO:  Seeding 100 elements within polygon number 46
06:42:43 INFO:  Seeding 100 elements within polygon number 47
06:42:43 INFO:  Seeding 100 elements within polygon number 48
06:42:43 INFO:  Seeding 100 elements within polygon number 49
06:42:43 INFO:  Seeding 100 elements within polygon number 50
06:42:43 INFO:  Seeding 100 elements within polygon number 51
06:42:43 INFO:  Seeding 100 elements within polygon number 52
06:42:43 INFO:  Seeding 100 elements within polygon number 53
06:42:43 INFO:  Seeding 100 elements within polygon number 54
06:42:43 INFO:  Seeding 100 elements within polygon number 55
06:42:43 INFO:  Seeding 100 elements within polygon number 56
06:42:43 INFO:  Seeding 100 elements within polygon number 57
06:42:43 INFO:  Seeding 100 elements within polygon number 58
06:42:43 INFO:  Seeding 100 elements within polygon number 59
06:42:43 INFO:  Seeding 100 elements within polygon number 60
06:42:43 INFO:  Seeding 100 elements within polygon number 61
06:42:43 INFO:  Seeding 100 elements within polygon number 62
06:42:43 INFO:  Seeding 100 elements within polygon number 63
06:42:43 INFO:  Seeding 100 elements within polygon number 64
06:42:43 INFO:  Seeding 100 elements within polygon number 65
06:42:43 INFO:  Seeding 100 elements within polygon number 66
06:42:43 INFO:  Seeding 100 elements within polygon number 67
06:42:43 INFO:  Seeding 100 elements within polygon number 68
06:42:43 INFO:  Seeding 100 elements within polygon number 69
06:42:43 INFO:  Seeding 100 elements within polygon number 70
06:42:43 INFO:  Seeding 100 elements within polygon number 71
06:42:43 INFO:  Seeding 100 elements within polygon number 72
06:42:43 INFO:  Seeding 100 elements within polygon number 73
06:42:43 INFO:  Seeding 100 elements within polygon number 74
06:42:43 INFO:  Seeding 100 elements within polygon number 75
06:42:43 INFO:  Seeding 100 elements within polygon number 76
06:42:43 INFO:  Seeding 100 elements within polygon number 77
06:42:43 INFO:  Seeding 100 elements within polygon number 78
06:42:43 INFO:  Seeding 100 elements within polygon number 79
06:42:43 INFO:  Seeding 100 elements within polygon number 80
06:42:43 INFO:  Seeding 100 elements within polygon number 81
06:42:43 INFO:  Seeding 100 elements within polygon number 82
06:42:43 INFO:  Seeding 100 elements within polygon number 83
06:42:43 INFO:  Seeding 100 elements within polygon number 84
06:42:43 INFO:  Seeding 100 elements within polygon number 85
06:42:43 INFO:  Seeding 100 elements within polygon number 86
06:42:43 INFO:  Seeding 100 elements within polygon number 87
06:42:43 INFO:  Seeding 100 elements within polygon number 88
06:42:43 INFO:  Seeding 100 elements within polygon number 89
06:42:43 INFO:  Seeding 100 elements within polygon number 90
06:42:43 INFO:  Seeding 100 elements within polygon number 91
06:42:43 INFO:  Seeding 100 elements within polygon number 92
06:42:44 INFO:  Seeding 100 elements within polygon number 93
06:42:44 INFO:  Seeding 100 elements within polygon number 94
06:42:44 INFO:  Seeding 100 elements within polygon number 95
06:42:44 DEBUG: 
------------------------------------------------------
Software and hardware:
  OpenDrift version 1.4.2
  Platform: Windows, 10
  15.881317138671875 GB memory
  4 processors (Intel64 Family 6 Model 142 Stepping 9, GenuineIntel)
  NumPy version 1.20.2
  SciPy version 1.6.2
  Matplotlib version 3.1.3
  NetCDF4 version 1.5.6
  Python version 3.7.8 | packaged by conda-forge | (default, Nov 27 2020, 18:48:03) [MSC v.1916 64 bit (AMD64)]
------------------------------------------------------

06:42:44 DEBUG: No output file is specified, neglecting export_buffer_length
06:42:44 INFO: Fallback values will be used for the following variables which have no readers: 
06:42:44 INFO:  x_wind: 0.000000
06:42:44 INFO:  y_wind: 0.000000
06:42:44 INFO:  sea_surface_wave_significant_height: 0.000000
06:42:44 INFO:  sea_surface_wave_stokes_drift_x_velocity: 0.000000
06:42:44 INFO:  sea_surface_wave_stokes_drift_y_velocity: 0.000000
06:42:44 INFO:  sea_surface_wave_period_at_variance_spectral_density_maximum: 0.000000
06:42:44 INFO:  sea_surface_wave_mean_period_from_variance_spectral_density_second_frequency_moment: 0.000000
06:42:44 INFO:  surface_downward_x_stress: 0.000000
06:42:44 INFO:  surface_downward_y_stress: 0.000000
06:42:44 INFO:  turbulent_kinetic_energy: 0.000000
06:42:44 INFO:  turbulent_generic_length_scale: 0.000000
06:42:44 INFO:  sea_floor_depth_below_sea_level: 10000.000000
06:42:44 INFO: Duration, steps or end time not specified, running until end of first reader: 2020-07-09 00:20:00
06:42:44 INFO: Duration, steps or end time not specified, running until end of first reader: 2020-07-09 00:20:00
06:42:44 DEBUG: Preparing readers for simulation coverage ([-52.17961120605469, -32.16594696044922, -52.05530548095703, -31.97176170349121]) and time (2020-07-09 00:00:00 to 2020-07-09 00:00:00)
06:42:44 DEBUG:     Preparing ./wind10_225deg_+2500.nc
06:42:44 DEBUG:     Preparing shape
06:42:44 INFO: Using existing reader for land_binary_mask
06:42:44 DEBUG: ----------------------------------------
06:42:44 DEBUG: Variable group ['land_binary_mask']
06:42:44 DEBUG: ----------------------------------------
06:42:44 DEBUG: Calling reader shape
06:42:44 DEBUG: ----------------------------------------
06:42:44 DEBUG: Data needed for 9500 elements
06:42:44 DEBUG: Reader time:
        None (before)
        None (after)
06:42:44 DEBUG: Fetching variables from shape
06:42:44 DEBUG: Fetched env-before
06:42:44 DEBUG: No time interpolation needed - right on time.
06:42:44 DEBUG: Reader SRS is the same as calculation SRS - rotation of vectors is not needed.
06:42:44 DEBUG: Masking 6176 elements outside coverage
06:42:44 DEBUG: Data missing for 6176 elements.
06:42:44 DEBUG: ---------------------------------------
06:42:44 DEBUG: Finished processing all variable groups
06:42:44 DEBUG: ------------ SUMMARY -------------
06:42:44 DEBUG:     land_binary_mask: 0 (min) 0 (max)
06:42:44 DEBUG: ---------------------------------
06:42:44 DEBUG:         0 active elements
06:42:44 INFO: Moving 0 out of 9500 points from land to water
06:42:44 DEBUG: ----------------------------------------
06:42:44 DEBUG: Variable group ['land_binary_mask']
06:42:44 DEBUG: ----------------------------------------
06:42:44 DEBUG: Calling reader shape
06:42:44 DEBUG: ----------------------------------------
06:42:44 DEBUG: Data needed for 560 elements
06:42:44 DEBUG: Reader time:
        None (before)
        None (after)
06:42:44 DEBUG: Fetching variables from shape
06:42:44 DEBUG: Fetched env-before
06:42:44 DEBUG: No time interpolation needed - right on time.
06:42:44 DEBUG: Reader SRS is the same as calculation SRS - rotation of vectors is not needed.
06:42:44 DEBUG: Obtained data for all elements.
06:42:44 DEBUG: ---------------------------------------
06:42:44 DEBUG: Finished processing all variable groups
06:42:44 DEBUG: ------------ SUMMARY -------------
06:42:44 DEBUG:     land_binary_mask: 0 (min) 1 (max)
06:42:44 DEBUG: ---------------------------------
06:42:44 DEBUG:         0 active elements
06:42:44 DEBUG: Cleaning up
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-3-f8629ba90ed2> in <module>()
     24 
     25 
---> 26 o.run()

D:\Dropbox\Doutorado\Jupyter\opendrift\opendrift\models\basemodel.py in run(self, time_step, steps, time_step_output, duration, end_time, outfile, export_variables, export_buffer_length, stop_on_error)
   2506 
   2507         self.interact_with_coastline(final=True)
-> 2508         self.state_to_buffer()  # Append final status to buffer
   2509 
   2510         #############################

D:\Dropbox\Doutorado\Jupyter\opendrift\opendrift\models\basemodel.py in state_to_buffer(self)
   2618                 getattr(self.elements, var)[element_ind]
   2619         # Copy environment data to history array
-> 2620         for i, var in enumerate(self.environment.dtype.names):
   2621             if self.export_variables is not None and \
   2622                     var not in self.export_variables:

**AttributeError: 'OceanDrift' object has no attribute 'environment'**

Previously I was using only one sigma layer from delft3D, so as not to have to convert to the z coordinate, and I had no problem running the opendrift.

My netCDF file:

{'latitude_longitude': <class 'netCDF4._netCDF4.Variable'>
 int32 latitude_longitude()
     grid_mapping_name: latitude_longitude
 unlimited dimensions: 
 current shape = ()
 filling on, default _FillValue of -2147483647 used,
 'depth': <class 'netCDF4._netCDF4.Variable'>
 float32 depth(depth)
     units: m
     positive: down
     axis: Z
     standard_name: depth
 unlimited dimensions: 
 current shape = (201,)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'time': <class 'netCDF4._netCDF4.Variable'>
 int32 time(time)
     long_name: time since initialization
     units: seconds since 2020-07-01 00:00:00
     calendar: proleptic_gregorian
     standard_name: time
 unlimited dimensions: 
 current shape = (4,)
 filling on, default _FillValue of -2147483647 used,
 'u': <class 'netCDF4._netCDF4.Variable'>
 float32 u(time, depth, Y, X)
     units: meter second-1
     standard_name: y_sea_water_velocity
     long_name: Sea water y velocity
 unlimited dimensions: 
 current shape = (4, 201, 502, 199)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'v': <class 'netCDF4._netCDF4.Variable'>
 float32 v(time, depth, Y, X)
     units: meter second-1
     standard_name: x_sea_water_velocity
     long_name: Sea water x velocity
 unlimited dimensions: 
 current shape = (4, 201, 502, 199)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'w': <class 'netCDF4._netCDF4.Variable'>
 float32 w(time, depth, Y, X)
     units: meter second-1
     standard_name: upward_sea_water_velocity
     long_name: W-velocity per layer in zeta point
 unlimited dimensions: 
 current shape = (4, 201, 502, 199)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'd': <class 'netCDF4._netCDF4.Variable'>
 float32 d(time, depth, Y, X)
     units: m2 s-1
     standard_name: ocean_vertical_diffusivity
     long_name: Vertical eddy diffusivity-3D in zeta point
 unlimited dimensions: 
 current shape = (4, 201, 502, 199)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 's': <class 'netCDF4._netCDF4.Variable'>
 float32 s(time, depth, Y, X)
     units: ppt
     standard_name: sea_water_salinity
     long_name: concentrations per layer in zeta point - Salinity
 unlimited dimensions: 
 current shape = (4, 201, 502, 199)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'lon': <class 'netCDF4._netCDF4.Variable'>
 float32 lon(Y, X)
     units: 
     standard_name: longitude
     long_name: longitude
 unlimited dimensions: 
 current shape = (502, 199)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'lat': <class 'netCDF4._netCDF4.Variable'>
 float32 lat(Y, X)
     units: 
     standard_name: latitude
     long_name: latitude
 unlimited dimensions: 
 current shape = (502, 199)
 filling on, default _FillValue of 9.969209968386869e+36 used}

My netCDF file+

knutfrode commented 3 years ago

Hi,

This log is a bit hard to understand. Are you seeding elements from a shapefile with 95 polygons? And you are using another shapefile for the landmask? Can you try o.plot() after your seeding, but before you do run()?

You might also need to do the following, which you are perhaps already doing:

o.set_config('general:use_auto_landmask', False)  # Disabling the automatic GSHHG landmask

Also it would help if you can post your script.

GuiCruz commented 3 years ago

Yes, I am seeding from a shapefile with around 100 polygons, and I am also using a custom shapefile for my landmask

seed_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/45degree/seeder.shp'
fn_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/Coastline_25x100km_GCS.shp'

o = OceanDrift(loglevel=0)
o.set_config('general:use_auto_landmask', False)
o.set_config('general:coastline_action', 'previous')

reader = reader_netCDF_CF_generic.Reader(fn_nc)
reader_shp = reader_shape.Reader.from_shpfiles(fn_shp)

o.add_reader([reader])
o.add_reader([reader_shp])

lats= np.linspace(-32.15, -32.15, 200)
#lons = np.linspace(-52.18633433, -51.97683977, 200)

o.seed_from_shapefile(seed_shp,number=9500, radius = 90,time=[reader.start_time, reader.end_time])
#o.seed_from_shapefile(seed_shp,number=9500, radius = 90,time=[datetime.datetime(2020, 7, 9, 0, 5), datetime.datetime(2020, 7, 9, 0, 20)])

o.set_config('environment:fallback:land_binary_mask', 0)

o.run()
knutfrode commented 3 years ago

This seems ok. Can you insert o.plot() before run, and post the graphics you get?

On Sun, Apr 18, 2021, 20:37 Guilherme Cruz @.***> wrote:

Yes, I am seeding from a shapefile with around 100 polygons, and I am also using a custom shapefile for my landmask

seed_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/45degree/seeder.shp' fn_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/Coastline_25x100km_GCS.shp'

o = OceanDrift(loglevel=0) o.set_config('general:use_auto_landmask', False) o.set_config('general:coastline_action', 'previous')

reader = reader_netCDF_CF_generic.Reader(fn_nc) reader_shp = reader_shape.Reader.from_shpfiles(fn_shp)

o.add_reader([reader]) o.add_reader([reader_shp])

lats= np.linspace(-32.15, -32.15, 200)

lons = np.linspace(-52.18633433, -51.97683977, 200)

o.seed_from_shapefile(seed_shp,number=9500, radius = 90,time=[reader.start_time, reader.end_time])

o.seed_from_shapefile(seed_shp,number=9500, radius = 90,time=[datetime.datetime(2020, 7, 9, 0, 5), datetime.datetime(2020, 7, 9, 0, 20)])

o.set_config('environment:fallback:land_binary_mask', 0)

o.run()

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/OpenDrift/opendrift/issues/593#issuecomment-822039227, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABH25IYWDWAHUEMUAH25INLTJMRFNANCNFSM43EBN6ZQ .

GuiCruz commented 3 years ago

Here is my complete code... Not sure if you will understand is a little confused:

import datetime
import utm
from pyproj import CRS
from opendrift.readers import reader_netCDF_CF_generic
from opendrift.readers import reader_netCDF_CF_unstructured
from opendrift.readers import reader_shape
from opendrift.readers import reader_global_landmask
from opendrift.models.oceandrift import OceanDrift
from scipy.interpolate import griddata
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D #
from matplotlib import cm
from metpy.interpolate import (interpolate_to_grid, remove_nan_observations,
                               remove_repeat_coordinates)
fn_nc = './wind10_225deg_+2500.nc'
fn  = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/5min_step_25x100km_Wind10_225deg_Rev1_+2500vaz.nc'#E:/5min_step_25x100km_Wind5_225deg_Rev1_0vaz.nc'
fn_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/Coastline_25x100km_GCS.shp'
ds = nc.Dataset(fn)
x = ds.variables['XCOR'][9:92,1:51]
y = ds.variables['YCOR'][9:92,1:51]
dpo = ds.variables['DP0'][9:92,1:51]
wls = ds.variables['S1'][2305:2309,9:92,1:51]

U = ds.variables['U1'][2305:2309,:,9:92,1:51]
V = ds.variables['V1'][2305:2309,:,9:92,1:51]
W = ds.variables['WPHY'][2305:2309,:,9:92,1:51]#W-velocity per layer in zeta point
diff = ds.variables['DICWW'][2305:2309,:,9:92,1:51]#Vertical eddy diffusivity-3D in zeta point
sal = ds.variables['R1'][2305:2309,0,:,9:92,1:51]#Concentrations per layer in zeta point - Salinity

dim0 = len(U)

dim2 = len(U[0][0])

dim3 = len(U[0][0][0])

x_f = np.ravel(x)
y_f = np.ravel(y)

dpo_f = np.ravel(dpo)
gx, gy, dpo_tp = interpolate_to_grid(x_f, y_f, dpo_f, interp_type='nearest', hres=100)#sea_floor_depth_below_sea_level

r_u = []
r_v = []
r_wl = []
r_W = []
r_d = []
r_s = []

t_u = []
t_v = []
t_W = []
t_d = []
t_s = []

r_m = []
i = 0
ii = 0
while i < dim0:
    wl_fp = np.ravel(wls[i,:,:])
    while ii < 10:
        u_fp = np.ravel(U[i,ii,:,:])
        v_fp = np.ravel(V[i,ii,:,:])
        w_fp = np.ravel(W[i,ii,:,:])
        s_fp = np.ravel(sal[i,ii,:,:])
        d_fp = np.ravel(diff[i,ii,:,:])

        for t in range (len(u_fp)):
            if u_fp[t]==-999:
                u_fp[t]=0

        for t in range (len(v_fp)):
            if v_fp[t]==-999:
                v_fp[t]=0

        gx, gy, u_tp = interpolate_to_grid(x_f, y_f, u_fp, interp_type='linear', hres=100)
        gx, gy, v_tp = interpolate_to_grid(x_f, y_f, v_fp, interp_type='linear', hres=100)
        gx, gy, w_tp = interpolate_to_grid(x_f, y_f, w_fp, interp_type='linear', hres=100)
        gx, gy, s_tp = interpolate_to_grid(x_f, y_f, s_fp, interp_type='linear', hres=100)
        gx, gy, d_tp = interpolate_to_grid(x_f, y_f, d_fp, interp_type='linear', hres=100)

        v_tp = v_tp*-1
        t_u.append(u_tp)
        t_v.append(v_tp)
        t_W.append(w_tp)
        t_d.append(d_tp)
        t_s.append(s_tp)
        ii+=1
    r_u.append(t_u)
    r_v.append(t_v)
    r_W.append(t_W)
    r_d.append(t_d)
    r_s.append(t_s)

    t_u = []
    t_v = []
    t_W = []
    t_d = []
    t_s = []    

    gx, gy, wl_tp = interpolate_to_grid(x_f, y_f, wl_fp, interp_type='linear', hres=100)
    r_wl.append(wl_tp)
    ii = 0
    i+=1
    print(i)

[gy_Deg,gx_Deg] = utm.to_latlon(gx, gy, 22, 'H')
gx_1 = gx_Deg[0,:]
gy_1 = gy_Deg[:,0]

gx_Deg = np.nan_to_num(gx_Deg)
gy_Deg = np.nan_to_num(gy_Deg)
r_wl = np.nan_to_num(r_wl)
r_u = np.nan_to_num(r_u)
r_v = np.nan_to_num(r_v)
r_W = np.nan_to_num(r_W)
r_d = np.nan_to_num(r_d)
r_s = np.nan_to_num(r_s)

dim0 = len(r_u)
dim2 = len(r_u[0][0])
dim3 = len(r_u[0][0][0])
z_reg = np.arange(0,20.1,0.1)

U_z = np.zeros((dim0,len(z_reg),dim2,dim3))
V_z = np.zeros((dim0,len(z_reg),dim2,dim3))
W_z = np.zeros((dim0,len(z_reg),dim2,dim3))
D_z = np.zeros((dim0,len(z_reg),dim2,dim3))
S_z = np.zeros((dim0,len(z_reg),dim2,dim3))

z_tsh = []

i = 0
ii = 0
iii = 0
while i < dim0:
    wl = r_wl[i][:][:]
    z_tsh = wl+dpo_tp
    while ii < dim2:
        while iii < dim3:
            z_sigma = [0.05*z_tsh[ii,iii], 0.15*z_tsh[ii,iii], 0.25*z_tsh[ii,iii], 0.35*z_tsh[ii,iii], 0.45*z_tsh[ii,iii], 0.55*z_tsh[ii,iii], 0.65*z_tsh[ii,iii], 0.75*z_tsh[ii,iii], 0.85*z_tsh[ii,iii], 0.95*z_tsh[ii,iii]] 
            u_sigma = [r_u[i,0,ii,iii],r_u[i,1,ii,iii],r_u[i,2,ii,iii],r_u[i,3,ii,iii],r_u[i,4,ii,iii],r_u[i,5,ii,iii],r_u[i,6,ii,iii],r_u[i,7,ii,iii],r_u[i,8,ii,iii],r_u[i,9,ii,iii]]
            v_sigma = [r_v[i,0,ii,iii],r_v[i,1,ii,iii],r_v[i,2,ii,iii],r_v[i,3,ii,iii],r_v[i,4,ii,iii],r_v[i,5,ii,iii],r_v[i,6,ii,iii],r_v[i,7,ii,iii],r_v[i,8,ii,iii],r_v[i,9,ii,iii]]
            w_sigma = [r_W[i,0,ii,iii],r_W[i,1,ii,iii],r_W[i,2,ii,iii],r_W[i,3,ii,iii],r_W[i,4,ii,iii],r_W[i,5,ii,iii],r_W[i,6,ii,iii],r_W[i,7,ii,iii],r_W[i,8,ii,iii],r_W[i,9,ii,iii]]
            d_sigma = [r_d[i,0,ii,iii],r_d[i,1,ii,iii],r_d[i,2,ii,iii],r_d[i,3,ii,iii],r_d[i,4,ii,iii],r_d[i,5,ii,iii],r_d[i,6,ii,iii],r_d[i,7,ii,iii],r_d[i,8,ii,iii],r_d[i,9,ii,iii]]
            s_sigma = [r_s[i,0,ii,iii],r_s[i,1,ii,iii],r_s[i,2,ii,iii],r_s[i,3,ii,iii],r_s[i,4,ii,iii],r_s[i,5,ii,iii],r_s[i,6,ii,iii],r_s[i,7,ii,iii],r_s[i,8,ii,iii],r_s[i,9,ii,iii]]

            U_reg = np.interp(z_reg,z_sigma,u_sigma,right=0)
            V_reg = np.interp(z_reg,z_sigma,v_sigma,right=0)
            W_reg = np.interp(z_reg,z_sigma,w_sigma,right=0)
            D_reg = np.interp(z_reg,z_sigma,d_sigma,right=0)
            S_reg = np.interp(z_reg,z_sigma,s_sigma,right=0)

            U_z[i,:,ii,iii] = U_reg#
            V_z[i,:,ii,iii] = V_reg#
            W_z[i,:,ii,iii] = V_reg#upward_sea_water_velocity - Opendrift
            D_z[i,:,ii,iii] = D_reg#ocean_vertical_diffusivity - Opendrift
            S_z[i,:,ii,iii] = S_reg#sea_water_salinity
            iii+=1

        iii = 0
        ii+=1
    ii = 0
    i+=1
    print(i)


t = ds.variables['time'][2305:2309]

DS_Open.close()
DS_Open = nc.Dataset(fn_nc, 'w', format='NETCDF4_CLASSIC', diskless=True, persist=True)
DS_Open.Conventions = 'CF-1.7'
DS_Open.title = '25x100_Wind05ms_opendrift'
DS_Open.institution = 'Unidata'
DS_Open.source = 'WRF-1.5'
DS_Open.history = ''
DS_Open.references = ''
DS_Open.comment = ''
DS_Open.createDimension('time', None)
DS_Open.createDimension('depth', z_reg.size)
DS_Open.createDimension('X', gx_1.size)
DS_Open.createDimension('Y', gy_1.size)

crs = DS_Open.createVariable('latitude_longitude', np.int32, ())
crs.grid_mapping_name = 'latitude_longitude'

Dph_var = DS_Open.createVariable('depth', np.float32, ('depth',))
Dph_var[:] = z_reg
Dph_var.units = 'm'
Dph_var.positive = 'down'
Dph_var.axis = 'Z' # Optional
Dph_var.standard_name = 'depth'

time_var = DS_Open.createVariable('t', np.int32, ('time',))
time_var[:] = t
time_var.long_name = 'time since initialization'
time_var.units = 'seconds since 2020-07-01 00:00:00'
time_var.calendar = 'proleptic_gregorian'
time_var.standard_name = 'time'  # Optional

U_var = DS_Open.createVariable('u', datatype=np.float32,
                              dimensions=('time', 'depth', 'Y', 'X'),
                              zlib=True)
U_var[:] = U_z
U_var.units =  'meter second-1'
U_var.standard_name = 'y_sea_water_velocity'
U_var.long_name = 'Sea water y velocity'

V_var = DS_Open.createVariable('v', datatype=np.float32,
                              dimensions=('time', 'depth', 'Y', 'X'),
                              zlib=True)
V_var[:] = V_z
V_var.units =  'meter second-1'
V_var.standard_name = 'x_sea_water_velocity'
V_var.long_name = 'Sea water x velocity'

W_var = DS_Open.createVariable('w', datatype=np.float32,
                              dimensions=('time', 'depth', 'Y', 'X'),
                              zlib=True)
W_var[:] = W_z
W_var.units =  'meter second-1'
W_var.standard_name = 'upward_sea_water_velocity'
W_var.long_name = 'W-velocity per layer in zeta point'

D_var = DS_Open.createVariable('d', datatype=np.float32,
                              dimensions=('time', 'depth', 'Y', 'X'),
                              zlib=True)
D_var[:] = D_z
D_var.units =  'm2 s-1'
D_var.standard_name = 'ocean_vertical_diffusivity'
D_var.long_name = 'Vertical eddy diffusivity-3D in zeta point'

S_var = DS_Open.createVariable('s', datatype=np.float32,
                              dimensions=('time', 'depth', 'Y', 'X'),
                              zlib=True)
S_var[:] = S_z
S_var.units =  'ppt'
S_var.standard_name = 'sea_water_salinity'
S_var.long_name = 'concentrations per layer in zeta point - Salinity'

lon_var = DS_Open.createVariable('lon', np.float32, ('Y', 'X'))
lon_var[:] = gx_Deg
lon_var.units = ''
lon_var.standard_name = 'longitude'  # Optional
lon_var.long_name = 'longitude'

lat_var = DS_Open.createVariable('lat', np.float32, ('Y', 'X'))
lat_var[:] = gy_Deg
lat_var.units = ''
lat_var.standard_name = 'latitude'  # Optional
lat_var.long_name = 'latitude'

#m_var = DS_Open.createVariable('mask', np.float32, ('Y', 'X'))
#m_var[:] = m_tp
#m_var.units = ''
#m_var.standard_name = 'mask'  # Optional
#m_var.long_name = 'mask'

DS_Open.variables
seed_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/45degree/seeder.shp'

o = OceanDrift(loglevel=0)
o.set_config('general:use_auto_landmask', False)
o.set_config('general:coastline_action', 'previous')

reader = reader_netCDF_CF_generic.Reader(fn_nc)
reader_shp = reader_shape.Reader.from_shpfiles(fn_shp)

o.add_reader([reader])
o.add_reader([reader_shp])

lats= np.linspace(-32.15, -32.15, 200)
#lons = np.linspace(-52.18633433, -51.97683977, 200)

o.seed_from_shapefile(seed_shp,number=9500, radius = 90,time=[reader.start_time, reader.end_time])
#o.seed_from_shapefile(seed_shp,number=9500, radius = 90,time=[datetime.datetime(2020, 7, 9, 0, 5), datetime.datetime(2020, 7, 9, 0, 20)])

o.set_config('environment:fallback:land_binary_mask', 0)

o.plot()
o.run()
GuiCruz commented 3 years ago

This seems ok. Can you insert o.plot() before run, and post the graphics you get? … On Sun, Apr 18, 2021, 20:37 Guilherme Cruz @.***> wrote: Yes, I am seeding from a shapefile with around 100 polygons, and I am also using a custom shapefile for my landmask seed_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/45degree/seeder.shp' fn_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/Coastline_25x100km_GCS.shp' o = OceanDrift(loglevel=0) o.set_config('general:use_auto_landmask', False) o.set_config('general:coastline_action', 'previous') reader = reader_netCDF_CF_generic.Reader(fn_nc) reader_shp = reader_shape.Reader.from_shpfiles(fn_shp) o.add_reader([reader]) o.add_reader([reader_shp]) lats= np.linspace(-32.15, -32.15, 200) #lons = np.linspace(-52.18633433, -51.97683977, 200) o.seed_from_shapefile(seed_shp,number=9500, radius = 90,time=[reader.start_time, reader.end_time]) #o.seed_from_shapefile(seed_shp,number=9500, radius = 90,time=[datetime.datetime(2020, 7, 9, 0, 5), datetime.datetime(2020, 7, 9, 0, 20)]) o.set_config('environment:fallback:land_binary_mask', 0) o.run() — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#593 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABH25IYWDWAHUEMUAH25INLTJMRFNANCNFSM43EBN6ZQ .

download

knutfrode commented 3 years ago

I guess the problem is that the temporal coverage of your current reader is only 15 minutes, which is even less than a calculation time step.

GuiCruz commented 3 years ago

That was the problem! with a time_step = 5 worked perfectly...

I got just another question for you. With the model sedimentdrift it is possible to set a "hindered settling" for the sediment? I need to simulate cohesive sediment and this parameter is important for me.

Thanks again for you support! Iam getting interesting results from opendrift... I hope to publish it soon

knutfrode commented 3 years ago

There is no built in feature for "hinders", but I guess you could add an artificial island to the landmsak reader, combined with coastline-interaction=previous.

On Sun, Apr 18, 2021, 21:04 Guilherme Cruz @.***> wrote:

That was the problem! with a time_step = 5 worked perfectly...

I got just another question for you. With the model sedimentdrift it is possible to set a "hindered settling" for the sediment? I need to simulate cohesive sediment and this parameter is important for me

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/OpenDrift/opendrift/issues/593#issuecomment-822043226, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABH25I6VBI5J2PT7FNMBLPLTJMULTANCNFSM43EBN6ZQ .