Skip to content

Commit

Permalink
clean up a bunch of warnings
Browse files Browse the repository at this point in the history
  • Loading branch information
gauteh committed Nov 4, 2019
1 parent 48b85e5 commit 0c852d5
Show file tree
Hide file tree
Showing 12 changed files with 122 additions and 131 deletions.
1 change: 1 addition & 0 deletions docs/.gitignore
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
build
_autosummary
source/autoapi
1 change: 1 addition & 0 deletions docs/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@

autoapi_type = 'python'
autoapi_dirs = [ '../../opendrift' ]
autoapi_keep_files = False # set to True when debugging autoapi generated files

extensions = [
"sphinx_rtd_theme",
Expand Down
15 changes: 8 additions & 7 deletions opendrift/elements/elements.py
Original file line number Diff line number Diff line change
@@ -1,17 +1,17 @@
# This file is part of OpenDrift.
#
#
# OpenDrift is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2
#
#
# OpenDrift is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
#
# You should have received a copy of the GNU General Public License
# along with OpenDrift. If not, see <http://www.gnu.org/licenses/>.
#
#
# Copyright 2015, Knut-Frode Dagestad, MET Norway

from collections import OrderedDict
Expand All @@ -33,19 +33,20 @@ class LagrangianArray(object):
types, e.g. oil, fish eggs...) to the core variables described below.
Attributes:
variables: An OrderedDict where keys are names of the
variables/properties of the current object. The values
of the OrderedDict are dictionaries with names such as
'dtype', 'unit', 'standard_name' (CF), 'default' etc.
All variable names will be added dynamically as attributes of
All variable names will be added dynamically as attributes of
the object after initialisation. These attributes will be
numpy ndarrays of same length, or scalars. The core variables
are:
ID: an integer identifying each particle.
status: 0 for active particles and a positive integer when deactivated
lon: longitude (np.float32)
lat: latitude (np.float32)
z: vertical position of the particle in m,
z: vertical position of the particle in m,
positive upwards (above sea surface)
"""

Expand Down Expand Up @@ -154,7 +155,7 @@ def extend(self, other):
new_data = getattr(other, var)

# If both arrays have an identical scalar, it remains a scalar
if (not isinstance(new_data, np.ndarray) and
if (not isinstance(new_data, np.ndarray) and
not isinstance(present_data, np.ndarray) and
present_data == new_data):
continue
Expand Down
3 changes: 2 additions & 1 deletion opendrift/models/basemodel.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@

import numpy as np
import scipy
import pyproj
import configobj, validate
try:
import matplotlib
Expand All @@ -52,7 +53,7 @@
print('matplotlib and/or cartopy is not available, can not make plots')

import opendrift
from opendrift.readers.basereader import pyproj, BaseReader, vector_pairs_xy
from opendrift.readers.basereader import BaseReader, vector_pairs_xy
from opendrift.readers import reader_from_url
from opendrift.models.physics_methods import PhysicsMethods

Expand Down
120 changes: 56 additions & 64 deletions opendrift/models/openberg.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,27 +21,19 @@
import numpy as np
import sys
from scipy.interpolate import interp1d

try:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib import animation
except:
logging.info('Basemap is not available, can not make plots')

import pyproj

from opendrift.models.basemodel import OpenDriftSimulation
from opendrift.elements.elements import LagrangianArray
from opendrift.readers.basereader import pyproj, BaseReader

from opendrift.readers.basereader import BaseReader


# Defining the iceberg element properties
class IcebergObj(LagrangianArray):
"""Extending LagrangianArray with variables relevant for iceberg objects.
"""

# We add the properties to the element class
variables = LagrangianArray.add_variables([
('wind_drift_factor', {'dtype': np.float32, # The fraction of the wind speed at which an iceberg is moved
Expand All @@ -63,10 +55,10 @@ class OpenBerg(OpenDriftSimulation):
"""

ElementType = IcebergObj

required_variables = ['x_wind', 'y_wind',
'x_sea_water_velocity', 'y_sea_water_velocity',
'land_binary_mask']
'land_binary_mask']

fallback_values = {'x_wind': 0.0,
'y_wind': 0.0,
Expand All @@ -81,119 +73,119 @@ class OpenBerg(OpenDriftSimulation):
def __init__(self, d=None, label=None, *args, **kwargs):

self.name = 'OpenBerg'
self.label=label
self.label=label

self.required_profiles = ['x_sea_water_velocity',
'y_sea_water_velocity'] # Get vertical current profiles

self.required_profiles_z_range = [-190, 0] # [min_depth, max_depth]

# Calling general constructor of parent class
super(OpenBerg, self).__init__(*args, **kwargs)


def update(self):
"""Update positions and properties of icebergs."""
"""Update positions and properties of icebergs."""
# Move icebergs at given % of wind speed
self.advect_wind()

# Move icebergs as per Anne Barker et al, 2004, with weighted average of current
# Move icebergs as per Anne Barker et al, 2004, with weighted average of current
# down to draft of berg.
# For current 10n meters down, where n is the depth index from 0 to 19

# RETRIEVE CURRENT SPEED:
y_sea_water_vel = self.environment_profiles['y_sea_water_velocity']
x_sea_water_vel = self.environment_profiles['x_sea_water_velocity']
x_sea_water_vel = self.environment_profiles['x_sea_water_velocity']

# FIND THE WEIGHTED AVERAGE CURRENT SPEED ACROSS THE ICEBERG KEEL
# FIND THE WEIGHTED AVERAGE CURRENT SPEED ACROSS THE ICEBERG KEEL
net_x_swv = np.zeros(x_sea_water_vel.shape[1])
net_y_swv = np.zeros(y_sea_water_vel.shape[1])
for indx in range(len(self.uw_weighting)):
net_y_swv = net_y_swv +y_sea_water_vel[indx,:]*self.uw_weighting[indx]
net_x_swv = net_x_swv +x_sea_water_vel[indx,:]*self.uw_weighting[indx]

self.update_positions(net_x_swv,net_y_swv)


def prepare_run(self):
""" Model spesific preparations.
Set the weighting for modelled current depths as per Table 5 of Barker 2004,
'Determination of Iceberg Draft, Mass and Cross-Sectional Areas',
Proceedings of The Fourteenth (2004) International Offshore and
Proceedings of The Fourteenth (2004) International Offshore and
Polar Engineering Conference.
NB! This version of OpenBerg does not allow for seeding of iceberg elements
of different sizes.
Also controles that the model handles readers without block data correctly.
"""
# Retrieve profile provided in z dimension by reader
# Retrieve profile provided in z dimension by reader
variable_groups, reader_groups, missing_variables = \
self.get_reader_groups(['x_sea_water_velocity','y_sea_water_velocity'])
self.get_reader_groups(['x_sea_water_velocity','y_sea_water_velocity'])

if len(reader_groups) == 0:
# No current data -> fallback values used
self.uw_weighting = np.array([1])
return
return

# Obtain depth levels from reader:
reader_name = reader_groups[0][0]
profile = np.abs(np.ma.filled(self.readers[reader_name].z))
# Make sure that interpolation is left to reader if no blocks are used .
# NB! This is a workaround, two additional points should be removed in basereader!

# Make sure that interpolation is left to reader if no blocks are used .
# NB! This is a workaround, two additional points should be removed in basereader!
if self.readers[reader_name].return_block == False:
self.use_block = False

# If current data is missing in at least one dimension, no weighting is performed:
if len(missing_variables) > 0:
logging.warning('Missing current data, weigthing array set to [1]')
self.uw_weighting = np.array([1])
return
return

# No need to create weighting array if only one z-level is provided:
if len(profile) == 1:
self.uw_weighting = np.array([1])
return

# Make copies to prevent outside value to be modified
water_line_length = self.elements_scheduled.water_line_length.copy()
depth = self.elements_scheduled.keel_depth.copy()

# Calculate the weighting array corresponding to the iceberg profile
uw_weighting = self.composite_iceberg(water_line_length=water_line_length, depth=depth)





#### Interpolate weighting array to match z-dimension of current reader ###

# Array of "Barker-depths" (10n)m, , where n is the depth index from 0 to 19
x = np.linspace(0,(len(uw_weighting)-1)*10,len(uw_weighting))
x = np.linspace(0,(len(uw_weighting)-1)*10,len(uw_weighting))

# Create a linear interpolator
interpol = interp1d(x,uw_weighting)
# Obtain corresponding depth levels from reader:
interpol = interp1d(x,uw_weighting)

# Obtain corresponding depth levels from reader:
self.reader_z_profile = profile[(profile >= 0) & (profile <= x.max())]

if len(profile[profile < 0]) > 0:
logging.warning('Current reader containing currents above water!'
'Weighting of current profile may not work!')

# Interpolate and normalize weights:
###### NB! Interpolator only includes values from reader within the range of the

###### NB! Interpolator only includes values from reader within the range of the
###### Barker-depth array. Meaning values just outside is not used for interpolation.
###### Ex.: If x.max=195 and the current reader includes data for
###### the z-profile: [0,3,10,15,25,50,75,100,150,200] only data from the
###### Ex.: If x.max=195 and the current reader includes data for
###### the z-profile: [0,3,10,15,25,50,75,100,150,200] only data from the
###### z-levels [0,3,10,15,25,50,75,100,150] are used.
interpol_weight = interpol(self.reader_z_profile)

interpol_weight = interpol(self.reader_z_profile)
normalized_weight = interpol_weight/interpol_weight.sum()

if normalized_weight.all():
self.uw_weighting = normalized_weight

super(OpenBerg, self).prepare_run()

def composite_iceberg(self, water_line_length=90.5, depth=60):
Expand All @@ -203,30 +195,30 @@ def composite_iceberg(self, water_line_length=90.5, depth=60):
in table 5 from Barker et. al.(2004).
"""
# Parameters from Baker et. al.:

a_param = [9.5173,11.1717,12.4798,13.6010,14.3249,13.7432,13.4527,15.7579,
14.7259,11.8195,11.3610,10.9202,10.4966,10.0893,9.6979,9.3216,8.9600,
8.6124,8.2783,7.9571]

b_param = [-25.94,-107.50,-232.01,-344.60,-456.57,-433.33,-519.56,-1111.57,-1125.00,
-852.90,-931.48,-1007.02,-1079.62,-1149.41,-1216.49,-1280.97,
-1342.95,-1402.52,-1459.78,-1514.82]

d = int(round(depth/10))

if d < 0: #Incase depth is given as a negative value
d = -d
d = -d


if d > len(a_param):
d = len(a_param)
print('##### OpenBerg does not support icebergs with keel depths greater than 200m!\n' +
'Using a composite iceberg with given waterline length and keel depth 200m')
'Using a composite iceberg with given waterline length and keel depth 200m')

area_list=[]

for i in range(0,d):

A = self.lin_func(a_param[i],b_param[i],water_line_length)
area_list.append(A)

Expand All @@ -237,10 +229,10 @@ def composite_iceberg(self, water_line_length=90.5, depth=60):
weigthing_array = np.array(area_list)/sum(np.array(area_list))

return weigthing_array


def lin_func(self,a,b,L):
"""Returns value of linear function A=aL+b."""
A = a*L + b

return A
12 changes: 6 additions & 6 deletions opendrift/models/openoil.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,12 @@
import numpy as np
from datetime import datetime
import logging
import pyproj
import matplotlib.pyplot as plt

from opendrift.models.basemodel import OpenDriftSimulation
from opendrift.elements import LagrangianArray
import opendrift.models.noaa_oil_weathering as noaa
from opendrift.readers.basereader import pyproj

try:
from itertools import izip as zip
Expand Down Expand Up @@ -191,7 +191,7 @@ def __init__(self, weathering_model='default', *args, **kwargs):

self.oil_weathering_model = weathering_model

# Update config with oiltypes
# Update config with oiltypes
oiltypes = [str(a) for a in self.oiltypes]
self._add_config('seed:oil_type', oiltypes,
'Oil type', overwrite=True)
Expand Down Expand Up @@ -219,7 +219,7 @@ def update_surface_oilfilm_thickness(self):
# Using stereographic coordinates to get regular X and Y
psproj = pyproj.Proj(
'+proj=stere +lat_0=%s +lat_ts=%s +lon_0=%s' %
(meanlat, meanlat, meanlon))
(meanlat, meanlat, meanlon))
X,Y = psproj(self.elements.lon[surface], self.elements.lat[surface])
mass_bin, x_edge, y_edge, binnumber = binned_statistic_2d(
X, Y, self.elements.mass_oil[surface],
Expand Down Expand Up @@ -446,7 +446,7 @@ def prepare_run(self):

self.oil_water_interfacial_tension = \
self.oiltype.oil_water_surface_tension()[0]
logging.info('Oil-water surface tension is %f Nm' %
logging.info('Oil-water surface tension is %f Nm' %
self.oil_water_interfacial_tension)
else:
logging.info('Using default oil-water tension of 0.03Nm')
Expand Down Expand Up @@ -1007,7 +1007,7 @@ def seed_from_geotiff_thickness(self, filename, number=50000,
'using present time: %s' % time)

ds = gdal.Open(filename)

srcband = ds.GetRasterBand(1)
data = srcband.ReadAsArray()

Expand Down Expand Up @@ -1059,7 +1059,7 @@ def seed_from_geotiff_thickness(self, filename, number=50000,
total_area[cat-1] = np.sum(areas)
layers[cat-1].DeleteFeature(outer)
layers[cat-1].ResetReading()

# Calculate how many elements to be seeded for each category
areas_weighted = total_area*thickness_microns
numbers = number*areas_weighted/np.sum(areas_weighted)
Expand Down
Loading

0 comments on commit 0c852d5

Please sign in to comment.