OpenDrift / opendrift

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

Postprocessing feature of radionuclides model #1128

Open ledgint opened 1 year ago

ledgint commented 1 year ago

I am trying to use postprocessing feature of radionuclides model, but running this code i had error.


from opendrift.readers import reader_ROMS_native, reader_global_landmask
from opendrift.models.radionuclides import RadionuclideDrift
from datetime import timedelta, datetime
import numpy as np

o = RadionuclideDrift(loglevel=0, seed=0)

reader_barents = reader_ROMS_native.Reader('https://thredds.met.no/thredds/dodsC/fou-hi/barents_eps_sdepth_be')

reader_landmask = reader_global_landmask.Reader()

o.add_reader([reader_barents,reader_landmask])

o.set_config('drift:vertical_mixing', True)

o.set_config('vertical_mixing:diffusivitymodel','environment')  # include settling without vertical turbulent mixing

o.set_config('vertical_mixing:timestep', 600.) # seconds
o.set_config('drift:horizontal_diffusivity', 10)

o.set_config('radionuclide:species:LMM', True)

o.set_config('radionuclide:particle_diameter',5.e-6)  # m
o.set_config('radionuclide:transformations:Kd',2.e0) # (m3/kg)

o.set_config('radionuclide:sediment:resuspension_depth',2.)
o.set_config('radionuclide:sediment:resuspension_depth_uncert',0.1)
o.set_config('radionuclide:sediment:resuspension_critvel',0.15)

o.set_config('radionuclide:transfer_setup','Bokna_137Cs')

o.set_config('general:coastline_action', 'previous')
o.set_config('general:seafloor_action','lift_to_seafloor')
o.set_config('seed:LMM_fraction',.45)
o.set_config('seed:particle_fraction',.55)

td=datetime.today()-timedelta(days=10)
time = datetime(td.year, td.month, td.day, 0)
latseed= 73.02;   lonseed=51.35

ntraj=5000
iniz=np.random.rand(ntraj) * -20. 

o.seed_elements(lonseed, latseed, z=iniz, radius=1000,number=ntraj,
                time=time, 
                )

o.run(steps=24*1, time_step=3600, time_step_output=3600)

print(o)
print('Final speciation:')
for isp,sp in enumerate(o.name_species):
    print ('{:32}: {:>6}'.format(sp,sum(o.elements.specie==isp)))

print('Number of transformations:')
for isp in range(o.nspecies):
    print('{}'.format(['{:>9}'.format(np.int32(item)) for item in o.ntransformations[isp,:]]) )

o.conc_lat  = reader_barents.get_variables('latitude',
                                       x=[reader_barents.xmin,reader_barents.xmax],
                                       y=[reader_barents.ymin,reader_barents.ymax])['latitude'][:]
o.conc_lon  = reader_barents.get_variables('longitude',
                                       x=[reader_barents.xmin,reader_barents.xmax],
                                       y=[reader_barents.ymin,reader_barents.ymax])['longitude'][:]
o.conc_topo  = reader_barents.get_variables('sea_floor_depth_below_sea_level',
                                       x=[reader_barents.xmin,reader_barents.xmax],
                                       y=[reader_barents.ymin,reader_barents.ymax])['sea_floor_depth_below_sea_level'][:]
o.conc_mask  = reader_barents.land_binary_mask

o.write_netcdf_radionuclide_density_map('radio_conc.nc', pixelsize_m=500.,
                                      zlevels=[-2.],
                                      activity_unit='Bq',
                                      horizontal_smoothing=True,
                                      smoothing_cells=1,
                                      time_avg_conc=True,
                                      deltat=2., # hours
                                      llcrnrlon=51.20, llcrnrlat=72.45,
                                      urcrnrlon=50.55, urcrnrlat=73.50,
                                     )

Error is:


11:29:57 INFO    opendrift.models.radionuclides:997: Postprocessing: Write density and concentration to netcdf file
11:29:57 INFO    opendrift.models.radionuclides:1034: z_array: ['-10000.0', '-2.0', '0.0']
Traceback (most recent call last):
  File "c:\GIT\opendrift\Radionuclides barents.py", line 86, in <module>
    o.write_netcdf_radionuclide_density_map('radio_conc.nc', pixelsize_m=500.,
  File "c:\GIT\opendrift\opendrift\models\radionuclides.py", line 1043, in write_netcdf_radionuclide_density_map
    self.get_radionuclide_density_array(pixelsize_m, z_array,
  File "c:\GIT\opendrift\opendrift\models\radionuclides.py", line 1304, in get_radionuclide_density_array
    outsidex, outsidey = max(x_array)*1.5, max(y_array)*1.5
                         ^^^^^^^^^^^^
ValueError: max() arg is an empty sequence
knutfrode commented 3 months ago

Perhaps @magnesim has a suggestion here?

magnesim commented 3 months ago

I think lower left corner actually must be the lower left corner. Try to set both the llcrnrlon<urcrnrlon and llcrcrlat<urcrnrlat