OceanParcels—Optimization
Lagrangian tracking of virtual ocean particles
Abstract
The OceanParcels Lagrangian Simulation package is used to optimize the lateral and longitudinal velocities of particles at given starting points such that they will be set on a fixed position after floating on the Agulhas Current of South Africa for three days.
The parameters are adjustable and easily changed. The optimization method that is used is Python built-in Powell's Method.
Installing Parcels
Installation was mainly done by following the guidelines in Parcels Documentation. However, my Python was outdated on my updated Ubuntu system, and apt-get update didn't update but to the old version of scipy et al.
$ cat requirements.txt
numpy>=1.9.1
scipy>=0.16.0
netcdf4>=1.1.9
flake8>=2.1.0
pytest>=2.7.0
pytest-ipynb
py>=1.4.27
pymbolic
cgen
cachetools>=1.0.0
xarray>=0.5.1
six>=1.10.0
enum34
progressbar2
The old versions of scipy needed to be removed using apt-get remove and after installed using pip command. I installed pip install --user numpy scipy matplotlib ipython jupyter pandas sympy nose and sudo python -m pip install --upgrade pip
OceanParcels needs recent versions to work but the version numbers are described in requirements.txt, which is shown on the image.
Using the commands
In [1]: import scipy
In [2]: scipy.__version__
Running
Running Parcels is straightforward, and the examples given in the package or Parcels Jupyter Notenbook.
The compilded Runge--Kutta kernel is fast, and reliable.
Letters for Simulations
First we need to define the letters in a separate file called molLetter.py. This script is used often in our Simulations.
lonB=[30,29, 29, 29, 30, 31, 31 ]
latB=[-37.5,-38,-37,-36, -39 , -38.5, -37.8 ]
lonC=[29.5, 29, 29, 29.7, 31.0, 30.5, 31 ]
latC=[-38.9,-38,-37.5,-37, -38.4 , -38.9, -37.2 ]
lonD1=[ 18, 18, 18, 18, 19, 20, 20.2, 20, 19 ] #Created from L
lonD = [x+2 for x in lonD1] #Translate
latD=[-39,-38,-37,-36, -36.5,-37,-37.6,-38,-38.8 ]
lonJ=[ 21.1, 20.4, 20.3, 20., 22, 22, 22, 22, 22, ]
latJ=[-39 ,-38.7,-38.3, -38 ,-37.5, -36,-37,-38,-38.5 ]
lonL=[29,29, 29, 29, 30, 31, 32 ]
latL=[-39,-38,-37,-36, -39 , -39, -39 ]
lonO=[ 25, 24, 24 , 25, 26, 27, 27, 26]
latO=[-39, -38.5, -37.5, -37, -37, -37.5, -38.5, -39]
lonM=[ 18, 18, 18, 18, 20, 22, 22, 22, 22, ]
latM=[-39,-38,-37,-36,-37.5,-36,-37,-38,-39 ]
from parcels import Grid, ParticleSet, Variable, JITParticle, AdvectionRK4
from scripts import plotTrajectoriesFile
import numpy as np
import math
from datetime import timedelta
from operator import attrgetter
import matplotlib
import sys
import molLetters
filenames = {'U': "examples/GlobCurrent_example_data/20*.nc",
'V': "examples/GlobCurrent_example_data/20*.nc"}
variables = {'U': 'eastward_eulerian_current_velocity',
'V': 'northward_eulerian_current_velocity'}
dimensions = {'lat': 'lat',
'lon': 'lon',
'time': 'time'}
grid = Grid.from_netcdf(filenames, variables, dimensions)
lon = molLetters.lonD+molLetters.lonO+molLetters.lonC
lat = molLetters.latD+molLetters.latO+molLetters.latC
pset = ParticleSet.from_list(grid=grid,
pclass=JITParticle,
lon=lon,
lat=lat)
pset.execute(AdvectionRK4,
runtime=timedelta(days=.1),
dt=timedelta(minutes=2),
interval=timedelta(hours=1),
output_file=pset.ParticleFile(name="LetterGlobCurrentParticles"))
plotTrajectoriesFile('LetterGlobCurrentParticles.nc',
tracerfile='examples/GlobCurrent_example_data/
20020101000000-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc',
tracerlon='lon',
tracerlat='lat',
tracerfield='eastward_eulerian_current_velocity',
mode='2d'
)
To set up the letters we run the Parcels code for a very short period of time. The letters were set by hand, and no automation or optimization was used. The number of points is such limited that it is not that time consuming.
Here, we use only the traditional Runge—Kutta Kernel because the runtime is such a short that the particles do not move significantly.
However, the plotting script needs to be adjusted to see the points clearly. As in import part of the Python code is shown, the plotting script is included in the scripts folder. We changed
plt.plot(np.transpose(lon), np.transpose(lat), '.', color='.-')
to
plt.plot(np.transpose(lon), np.transpose(lat), '.', color='k')
on line 70 of scripts/plotParticles.py to show black dots instead of dotted line..
Kernels
Parcels allows to use and add user defined Kernels, e.g. to take the wind into account. Following the guides in Ocean Parcels Jupyter Notebook/ periodic boundaries we implement periodic boundary conditions in longitudinal direction. This was needed only on a preliminary simulations, and the final optimization algorithm does not need this.
def periodicBC(particle, grid, time, dt):
particle.lon = particle.lon
if particle.lon < 15:
particle.lon = 32.5
if particle.lon > 32.5:
particle.lon = 15
Here, the longitudinal position of particles is restricted to between 15 and 32.5 degrees East. Following the Jupyter notebook/ Adding a custom behaviour kernel, we define the wind velocity Kernel for each particles independently:
def WestVel(particle, grid, time, dt):
particle.lon += particle.vlon * dt / 1.1e5 # 1.1e5 m in 1 degree latitude
particle.lat += particle.vlat * dt / 1.1e5 # 1.1e5 m in 1 degree latitude, not correct but ok for the test.
We use the vlon and vlat properties of particle objects to set a wind velocity (or motor direction) to a particle. However, the particles are now our own defined VelParticles instead of JITParticles or SciPyParticles. VelParticles are defined as adding new variables to JITParticle:
class VelParticle(JITParticle):
vlon = Variable('vlon', initial=0., dtype=np.float32)
vlat = Variable('vlat', initial=0., dtype=np.float32)
Now, while setting up the particle's position, we need to use VelParticles in ParticleSet -command:
pset = ParticleSet.from_list(grid=grid
pclass=VelParticle, # (JITParticle or ScipyParticle)
lon=lon,
lat=lat)
We need to define the Kernels to be real Kernels
k_WestVel = pset.Kernel( WestVel )
k_dist = pset.Kernel( initialVel )
k_BC = pset.Kernel( periodicBC )
to be used in Simulation:
pset.execute(AdvectionRK4 + k_WestVel + k_BC + k_dist,
runtime=timedelta(days=3),
dt=timedelta(minutes=2),
interval=timedelta(hours=1),
output_file=pset.ParticleFile(name="GlobCurrentParticles"))
The Kernel called initialVel is defined, again, in a separate file. This is due to the fact that the Kernels are to be compiled and thus changing the values in Kernels seems to be difficult. While it is in a separate file it is possible to edit the file and recompile on a every optimization step.
The file initialVelocity.py reads
def initialVel(particle, grid, time, dt):
if particle.id == 0:
particle.vlon = 1.71528674607
particle.vlat = 0.0969265914145
if particle.id == 1:
particle.vlon = 1.2606618507
particle.vlat = -0.273957677955
.
.
.
if particle.id == 22:
particle.vlon = -0.694711688616
particle.vlat = -0.144419453392
if particle.id == 23:
particle.vlon = -0.624652167656
particle.vlat = 0.130573448637
The longitudinal and latitudinal velocities of each particles are defined in the File. We have 24 particles, and thus the last index is 23.
Distance / Cost Function
We use Manhattan, taxi or 1-norm to define the distance between current position of particles and the indented position of particles. This is to be minimized. The distance function is
def distance(final, target):
# Use 1-norm (Manhattan, or taxi norm) because it is very fast.
elon = np.zeros( len(final.lon) )
elat = np.zeros( len(lon) )
for i in range( len( lon )):
elon[i] = abs( final.lon[i] - target.lon[i] )
elat[i] = abs( final.lat[i] - target.lat[i] )
#print( lon[i], lat[i] )
return elon, elat
The function uses an object called final and target. These are defined as
class target:
name = 0
def __init__(self, name, lon, lat):
self.name = name
self.lon = lon
self.lat = lat
and
job = target("Job", molLetters.lonJ+molLetters.lonO+molLetters.lonB, molLetters.latJ+molLetters.latO+molLetters.latB )
and, finally, after the Simulation we extract the new positions as
lons = np.array([particle.lon for particle in pset.particles])
lats = np.array([particle.lat for particle in pset.particles])
fin = target("Job", lons, lats )
## Calculate and print the distance (error function)
elon, elat = distance( fin, job )
Thus, this gives the longitudinal and latitudinal distance from current point to the point we are aiming at.
Optimization
Next, we need to optimize the Velocities in the Kernel. Each particle needs different velocity, and this are optimized. First we note that using bad values (e.g. if the latitudinal coordinate of a particle is too far away from Agulhas current, the Lagrangian does not work. This is due to the fact the current is not defined too far away.
Thus, the Simulation fails at some initial velocities. We can handle this by Python's try command:
try:
pset.execute(AdvectionRK4 + k_WestVel + k_BC + k_dist,
runtime=timedelta(days=3),
dt=timedelta(minutes=2),
interval=timedelta(hours=1),
output_file=pset.ParticleFile(name="GlobCurrentParticles"))
except:
print( "Virhe" )
If the Simulation raises an exception, this will print an error message "Virhe".
The output, also the distance to the target is then printed to the screen, or to the stdout stream.
elon, elat = distance( fin, job )
printArray( elon )
printArray( elat )
This is done in a function execution code, in a separate file. The optimization code reads the output, uses Powell's optimization method to set new variables, writes those to the initialVelocity.py script and runs the function execution code, as is shown in the Image.
The optimization is set up rather easily in Python. Just write
res = minimize( fun_min, x0, method='powell',
options={'xtol':1e-3, 'disp': True})
This need the function fun_min which is defined below:
def fun_min(x):
#Input x: lon + lat = 24 + 24
#Convert x to particles:
particles = []
for i in range( 24 ):
particles.append( particle(i, x[2*i], x[2*i+1] ))
#print( particles[i].vlon )
#print( particles[i].vlat )
#Write the input file
createInitialVeloFile( particles )
#Do the Parcels simulation
out = check_output(["python", "molDistance.py"])
#print( out )
elon, elat, ok = getElonElat( out )
summa = 1e6
if ( ok ):
summa = ( sum( elon) + sum( elat ) )
return summa
This converts the x variable to particles object which is used to print the initialVelocity.py file. The output of simulation is then parsed and if the simulation was not successlul, the fun_min returns a very big number. Because this is actually 24 × 2 dimensional problem it could be simplified optimizing only two variables a time.
To write the initialVelocity.py file we use
def createInitialVeloFile( p ):
f = open( 'initialVelocity.py','w')
f.write('def initialVel(particle, grid, time, dt):' + '\n')
for i in range( len(p)):
f.write( ' if particle.id == ' + str( p[ i ].id ) + ':\n' )
f.write( ' particle.vlon = ' + str( p[ i ].vlon ) + '\n')
f.write( ' particle.vlat = ' + str( p[ i ].vlat ) + '\n')
The output is parsed using scripts
def getFloat( inp, si ):
ind1 = string.find( inp, '[', si )
ind2 = string.find( inp, ']', si+1 )
#print( si, " ", ind1, " ", ind2 )
Str = inp[ind1+1: ind2]
Arr = filter( bool, Str.split(' ') ) #Remove empty string elements
e = [float( li ) for li in Arr] #Convert to float
return e, ind2
def getElonElat( inp ):
#Check if inp contains error: "Virhe"
ok = True
elon = []
elat = []
if ( string.find( inp, 'Virhe' ) > 0 ):
ok = False
else:
elon, si = getFloat( inp, 1)
elat, si = getFloat( inp, si)
return elon, elat, ok
The optimization is very rude, and is written above, it can further optimize rather easily. For example, the Letter O is similar in the optimization and thus needs to be optimized only once.
Plotting the Animation
The animation is plotted using scripts/plotParticles.py script. However, the script needed some minor modifications.
The points were originally multicolored and plotted using dashes. This was changed to black dots, also
elif mode == '2d':
if type == 'array':
plt.plot(np.transpose(lon), np.transpose(lat), '.', color='k') #Used to be .- EDITED By MOL
plt.contourf(np.squeeze(X), np.squeeze(Y), np.squeeze(P))
The animation is saved using Imagemagick and
anim.save('testiAnimaatio.gif',writer='imagemagick', fps=60)
The animation is plotted as an animated gif. The gif frames were extracted using GIMP and saved as individual png images. These images are opened in Blender Video Sequence Editor (VSE) to have some modifications. We add the letters and slow the animation when the text should be readable.
The final animation is then extracted as png images and as a video to be uploaded to Youtube. The png images are imported again into GIMP and optimized as a gif animation. This can be automated using Sceme, as we have shown in a separate html document. Unfortunately, this is only in Finnish.
Finally, we get the small gif animation and high resolution Youtube video.
Visualizating the Wind
To visualize the driving force (or very local wind) of moving particles, we use Python's Quiver plot option. It allows to plot a vector field on a Figure, or as in this case, the vector field represents only a fixed number of points.
The driving force is constant, and thus we need to apply its set_offsets property in Animation script.
The vector field of the drifting Agulhas Current is not shown. It could be visualized in a similar way using Quiver plot.
from netCDF4 import Dataset
import numpy as np
import JobDocVelocity as inp
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import struct
class particle(object):
id=0
vlon=0
vlat=0
def animate(t):
scat.set_offsets(np.matrix((lon[:, t], lat[:, t])).transpose())
arrs.set_offsets(np.matrix((lon[:, t], lat[:, t])).transpose())
scale = 1.
# Background image (water)
#
tracerfile = 'examples/GlobCurrent_example_data/
20020101000000-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc'
tfile = Dataset(tracerfile, 'r')
tracerfield='eastward_eulerian_current_velocity'
X = tfile.variables['lon']
Y = tfile.variables['lat']
P = tfile.variables[tracerfield]
plt.contourf(np.squeeze(X), np.squeeze(Y), np.squeeze(P))
#
# Read the data and set up variables
#
pfile = Dataset('JobDoc.nc', 'r')
lon = pfile.variables['lon']
lat = pfile.variables['lat']
z = pfile.variables['z']
u,v = [],[]
for i in range(len(lon)):
p = particle()
p.id = i
inp.initialVel( p, 0, 0, 0)
u.append( p.vlon )
v.append( p.vlat )
fig = plt.figure(1)
ax = plt.axes()
scat = ax.scatter(lon[:, 0], lat[:, 0], s=60, cmap=plt.get_cmap('autumn'), color='k')
arrs = ax.quiver( lon[:,0], lat[:,0], u, v )
frames = np.arange(1, lon.shape[1])
ax.set_xlim([15,32.5])
ax.set_ylim([-40,-33])
#
# Animation
#
anim = animation.FuncAnimation(fig, animate, frames=frames,
interval=100, blit=False)
anim.save('testiAnimaatioVekt.gif',writer='imagemagick', fps=60)
plt.show()
First, we need to import the libraries, and this time we import data file also.
The function
def animate(t):
scat.set_offsets(np.matrix((lon[:, t], lat[:, t])).transpose())
arrs.set_offsets(np.matrix((lon[:, t], lat[:, t])).transpose())
The background image is loaded in shown as a contour plot.
The data is used to fill u and v variables. Those are read using inp.initialVel function, and appended to the list, or array.
The plots are easily made, and the scatter plot and quiver plot understands arrays, so that they are efficient and easy to use.
The animation function needs figure to update, animate function to update the frames, and some other parameters.
Also, we used Imagemagick to print the animation as an animated gif to be used in image processing program.
Results
The data and files are shown in a matrix to clarify the transitions from a position to another.
To further empahasize the initial velocity (or constant power applied) of the fluid parcel representing the bit of a font we show that on a similar grid as befire.
| MoL | Doc | Job | |
|---|---|---|---|
| Mol |
|
| |
| Doc |
|
| |
| Job |
|
|