In this section, the different Ichthyop processes are described. These processes are implemented in the classes within the action folder.
Parameters associated with these processes must be included in action blocks. They can be activated or deactivated by setting the enabled tag (cf. Section 2.1).
There is 2 different implementations of the growth module. They can be selected by setting the class_name parameter in the action.growth action configuration block.
To use the linear growth method described in Lett et al. (2008), choose the LinearGrowthAction.java class.
Length increment is provided as follows:
\[ \Delta L = C_1 + C_2 \times \dfrac{F}{F + K_S} \times max(T, T_{thres}) \times \Delta t \]
where \(\Delta t\) is the time-step (in days), \(C_1\) and \(C_2\) are parameters (coeff1 and coeff2), \(T_{thres}\) is a temperature threshold (threshold_temp parameter), \(F\) is the food quantity and \(K_S\) is a half-saturation constant (half_saturation parameter). If the latter is not defined or equals 0, \(Q\) is assumed to be \(1\).
The name of the food and temperature variables are provided by the food_field and temperature_field parameters.
If the {samp}SoleGrowthAction.java class is selected, the growth model used in Tanner et al. (2017) is used. It relies on the growth equation from Fonds (1979) and given by:
\[ \Delta L = C_1 \times T^{C_2} \times \Delta t \]
with \(L\) the length, \(T\) the temperature and \(\Delta t\) the time step in days, and \(C_1\) and \(C_2\) are parameters (c1 and c2 respectively).
The temperature field is provided by the temperature_field parameter.
In Ichthyop, there is the possibility to define a range of temperature and salinity beyond which the particle is killed.
The functionning og the lethal temperature module depends on whether the growth module is enabled or not.
lethal temperature can either be provided in a CSV file (lethal_temp_file parameter), which provides the lower and upper temperature values that can be supported by the particle as a function of age (in hours). The file must have the following format:
Time(hour);Cold temperature (C);Warm temperature (C)
0;14;22
48;13;22
96;12;22In this case, three age classes will be considered: \([0, 48[\), \([48, 96[\) and \([96, \infty[\)
If the lethal_temp_file parameter is not defined, single values will be used independently of the age. These values are provided in the cold_lethal_salinity_egg and warm_lethal_salinity_egg parameters.
Note that the temperature_field parameter, providing the name of the temperature variable, must also be provided.
If the growth module is enabled, two cold and warm lethal temperatures must be provided. One for eggs (cold_lethal_temperature_egg and hot_lethal_temperature_egg), one for larva (cold_lethal_temperature_larva and hot_lethal_temperature_larva).
The stage (egg or larva) is determined by the growth module, and the right temperature range is applied to the particle.
Lethal salinity can either be provided in a CSV file (lethal_salt_file parameter), which provides the lower and upper salinity values that can be supported by the particle as a function of age (in hours). The file must have the following format:
Time(hour);Fresh salinity (PSU);Salty salinity (PSU)
0;35;40
48;30;40
96;30;45In this case, three age classes will be considered: \([0, 48[\), \([48, 96[\) and \([96, \infty[\)
If the lethal_salt_file parameter is not defined, single values will be used independently of the age. These values are provided in the fresh_lethal_salinity_egg and saline_lethal_salinity_egg parameters.
Note that the salinity_field_ parameter, providing the name of the salinity variable, must also be provided.
The buoyancy module allow to vertically displace a particle, depending on it’s density and on the sea water density, following Parada et al. (2003).
The buoyancy-induced vertical velocity of the particle is given by:
\[ W_{buoy} = \dfrac{1}{24} \times g \times a \times b \times \dfrac{\rho_{water} - \rho_{part}}{\rho_{water}} \mu^{-1} \log\left(2 \dfrac{a}{b} + 0.5\right); \]
with \(a\) and \(b\) the semi-major axis and semi-minor axis of an ellipse (mean_major_axis and mean_minor_axis parameters), \(\mu\) the molecular viscosity (molecular_viscosity parameter), \(\rho_{water}\) the water density, \(\rho_{part}\) the particle density and \(g\) the gravitational acceleration (\(cm.s^{-2}\)).
If the particle density varies with age, it can be provided in a CSV file (density_file parameter), formatted as follows:
Age(hour);Density (g/cm3)
0;1.0235625
24;1.02374
48;1.023335
60;1.0239
66;1.025672If the density_file parameter is not found, a constant density will be assumed (particle_density parameter).
The buoyancy process is only applied for the early life stages. If the growth process is disabled, it is controlled by the age_max parameter (provided in days). If this parameter is not provided, the buoyancy process will always be applied.
If the growth process is enabled, the application of the buoyancy module will be automatically managed.
The MigrationAction.java manages the daily migration of the particles. It is done by providing daytime and nighttime depths, and the timing of the sunset and sunrise.
If the daytime_depth_file parameter is defined, it provides the depth of the particle during daytime. It is formatted as follows:
Age (day);Depth (m)
0.0;-20
3.0;-25
5.0;-30
8.0;-35If this parameter is not found, a constant daytime depth, provided by the daytime_depth parameter, is assumed.
Same thing for the nighttime depths, which can be set by either the nighttime_depth_file or the nighttime_depth parameters.
The sunset and sunrise hours are set by the sunset and sunrise parameters, which must have a HH::mm format.
If the growth module is deactivated, the user must provide the minimum age (in days) at which the particle starts to migrate (age_min parameter). If the growth module is activated, it manages the activation or deactivation of the daily migration.
When the target depth is greater than the total depth, the particle does not move.
The ontogenetic migration module controls the migration to different habitats as a function of age. Its implementation in Ichthyop follows the CMS one. It uses the CMS configuration file (cms_ovm_config_file parameter), which is formatted as follows:
nTime
nDepth
z1 z2 z3 z4 z5 z6 ... znDepth
t1 t2 t3 t4 t5 t6 ... tnTime
P1,1 P1,2 ... P1,nTime
P2,1 P2,2 ... P2,nTime
P3,1 P3,2 ... P3,nTime
...
PnDepth,1 PnDepth,2 ... PnDepth,nTimeThe first line provides the number of time steps in the CMS file, the second line provides the number of vertical levels in the CMS file. The third line provides the depth values and the fourth lines provides the time step values.
The remaining lines provide the probability matrix \(P_{z, t}\).
The sum of the probability matrix along the depth dimension should equal 100.
At each time step, the index of the CMS time step, \(k\), is determined by comparing the simulation and CMS times as follows:
\[ t_{CMS}(k) < t \leq t_{CMS}(k + 1) \]
When the CMS time index changes, all the particles are randomly distributed on the vertical, following the probability distribution of the given CMS time. If the CMS time index remains unchanged, nothing is done.
The wind-drift is implemented in the WindDriftFileAction.java class. This class requires that 2D (time, latitude, longitude) wind fields are provided. This is done by providing a input_path and a file_filter parameter, which specify the location and names of the wind files.
The user also provides field_time, wind_u, wind_v, longitude, latitude parameters, which specify the names of the time, zonal wind, meridional wind, longitude and latitude variables in the NetCDF files.
The user also provides a depth_application parameter, which specifies the depth at which the wind will impact the trajectories (only valid for 3D simulations), a wind_factor (\(F\), multiplication factor), an angle (\(\theta\)) and a wind_convention parameter (\(\varepsilon\), +1 if convention is ocean-based, i.e. wind-to, -1 if convention is atmospheric based, i.e. wind-from).
The changes in particle longitude (\(\lambda\)) and latitude (\(\phi\)) is provided as follows:
\[ \Delta \lambda_0 = \frac{\Delta t \times U_W }{\frac{R \pi}{180} \times \cos(\frac{\pi \phi}{ 180})} \]
\[ \Delta \phi_0 = \frac{\Delta t \times V_W}{ \frac{R \pi}{180}} \]
\[ \Delta \lambda = \varepsilon \times F \times \left(\Delta \lambda_0 \times \cos\left(\frac{\theta \pi}{180}\right) - \Delta \phi_0 \times \sin\left(\frac{\theta \pi}{180}\right)\right) \]
\[ \Delta \phi = \varepsilon \times F \times \left(\Delta \lambda_0 \times \sin\left(\frac{\theta \pi}{180}\right) + \Delta \phi_0 \times \sin\left(\frac{\theta \pi}{180}\right)\right) \]
\(R\) is the Earth Radius and equals 6367.74 km. In the code, \(\dfrac{R \pi}{180}\), which is the distance in m of a 1 degree cell, is approximated to 111138 m
Ichthyop allows particles to randomly swim. This is achieved by using the SwimmingAction.java. The velocity may vary with the age of the particle.
The user provides a CSV absolute velocity file (velocity_file parameter), which must me formatted as follows:
Age (days);Speed (m/s)
0;0.1
5;0.2
15;0.3The user also provides a boolean parameter (constant_velocity) that specifies whether the velocity should remain as defined in the CSV file, or if the absolute velocity should be randomly selected as follows:
\[ \|U\| = \|U\|_{file} \times \kappa \]
with \(\kappa\) a random value in the \([0, 2]\) interval and \(\|U\|_{file}\) the velocity defined in the CSV file.
At each time-step, the zonal and meridonal velocities are defined as follows:
\[ U = \kappa' \varepsilon \|U\| \]
\[ V = \varepsilon' \sqrt{\|U\|^2 - U^2} \]
with \(\kappa'\) a random value between 0 and 1, and \(\varepsilon\) and \(\varepsilon'\) random values either equal to 1 or -1.
Ichthyop can take into account the effects of waves on the particles trajectories, following the Stokes drift equations of Stokes (2009):
\[ U_{S} = w \times k \times a^2 \times \exp\left(2 \times k \times z\right) \]
In Ichthyop, the user provides zonal and meridional wave stokes drift and the wave periods. The horizontal displacement of particles is computed as follows:
\[ \|U_{wave}\| = U_{wave}^2 + V_{wave}^2 \]
\[ \lambda_{wave} = \|U_{wave}\| \times T_{wave} \]
\[ k_{wave} = \dfrac{2 \pi}{\lambda_{wave}} \]
\[ \Delta X = F \times U_{cur} \times \Delta t \times \exp\left(2 \times k_{wave} \times z\right) \]
\[ \Delta Y = F \times V_{cur} \times \Delta t \times \exp\left(2 \times k_{wave} \times z\right) \]
with \(U_{wave}\) and \(V_{wave}\) the zonal and meridional stokes components, \(T_{wave}\) the wave period, \(U_{cur}\) and \(V_{cur}\) the zonal and meridional ocean current components, \(z\) the depth and \(F\) a multiplication factor provided by the user.
In the bouncing mode, the particle will bounce on the coast. First, whether the bouncing occurs on a meridional or a zonal coastline is determined.
In case of a meridional coastline, as shown in Figure 6.1, the calculation of the new position is performed as follows.
Let’s assume that the particle is at the position \((x, y)\) and is moved at the postion \((x + \Delta x, y + \Delta y)\). We suppose that
\[ \Delta y = \Delta y_1 + \Delta y_2, \tag{6.1}\]
where \(\Delta y_1\) is the meridional distance between the particle and the coastline, and \(\Delta y_2\) is the distance that the particle will spend on land.
In the bouncing mode, the position increment can be written as
\[ \Delta_{cor} y = \Delta y_1 - \Delta y_2 \tag{6.2}\]
By replacing \(\Delta y_2\) using Equation 6.1, we can write:
\[ \Delta_{cor} y = \Delta y_1 - (\Delta y - \Delta y_1) \]
\[ \Delta_{cor} y = 2 \Delta y_1 - \Delta y \]
Active swimming has been implemented in Ichthyop following the work of Romain Chaput. Three active swimming behaviours have been implemented: the rheotaxis orientation (i.e. against the current), the cardinal orientation (i.e toward a given direction) and the reef orientation, i.e. orientation toward points of interests.
These three implementations involve the computation of a swimming velocity and direction. The former is common to all three methods, but the directions depend on the method considered.
The orientation processes all share common features. They both depend on a swimming velocity and random directions. The methods described below differ on the way the random directions are drafted.
Swimming velocity is computed following Staaterman, Paris, and Helgers (2012):
#|echo: false \[ V = V_{hatch} + (V_{settle} - V_{hatch}) ^ {\log(age) / \log(PLD)} \]
with \(V_{hatch}\) and \(V_{settle}\) the larval velocity at hatching and settle, \(A\) the age of the larva and \(PLD\) the transport duraction.
Von Mises distribution is used in all three methods.The Von Mises distribution is given by:
\[ f(\theta, \mu, \kappa) = \dfrac {\exp(\kappa \cos(\theta - \mu))} {2 \pi I_{0}(\kappa)} \]
where \(I_{0}(\kappa)\) is the modified Bessel function of the first kind of order 0, \(\mu\) is angle where the distribution is centerred and \(\kappa\) is the concentration parameter. The distribution is as follows:
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import i0
x = np.linspace(-np.pi, np.pi, 200)
def von_misses(kappa):
mu = 0
y = np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa))
return y
plt.figure()
for i in [0.5, 1, 2, 5, 10, 20]:
plt.plot(x, von_misses(i), label=f'$\kappa = {i}$')
plt.legend()
plt.xlim(x.min(), x.max())
plt.show()
For computation purposes, all the Von Mises drafts performed in Ichthyop are done by using \(mu = 0\). The angles are thus centerred around 0. Then, the \(mu\) value is added.
\(\theta = 0\) is eastward, \(\theta = \frac{\pi}{2}\) is northward, etc.
Given a swimming velocity \(V\) and a direction \(\theta\), the larva displacement (in \(m\)) is computed as follows:
\[ \Delta X = V \times \cos(\theta) \times \Delta t \]
\[ \Delta Y = V \times \sin(\theta) \times \Delta t \]
with \(\Delta t\) the time step. Next, the corresponding change in longitude (\(\lambda\)) and latitude (\(\varphi\)) is computed as follows:
\[ \Delta \lambda = \dfrac{\Delta X}{111138 \times \cos{\varphi}} \]
\[ \Delta \varphi = \dfrac{\Delta Y}{111138 } \]
In cardinal orientation, the user provides a fixed heading \(\theta_{card}\) and a fixed \(\kappa\) parameter. Then, at each time step, a new angle is randomly drafted following a Von Misses distribution \(f(\theta, \theta_{card}, \kappa)\).
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import os
from glob import glob
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
config = 'card'
if os.path.isdir('ichthyop_output'):
pattern = os.path.join('ichthyop_output', config, '*nc')
else :
pattern = os.path.join('..', '..', 'ichthyop_output', config, '*nc')
pattern
filelist = glob(pattern)
filelist
data = xr.open_mfdataset(filelist, decode_times=False)
data
# +
lon = data['lon'].values
lat = data['lat'].values
mort = data['mortality'].values
lon = np.ma.masked_where(mort > 0, lon)
lat = np.ma.masked_where(mort > 0, lat)
ntime, ndrifters = lon.shape
# +
time = np.arange(ntime)
drifters = np.arange(ndrifters)
d2d, t2d = np.meshgrid(drifters, time)
# +
plt.figure()
projin = ccrs.PlateCarree()
projout = ccrs.PlateCarree()
ax = plt.axes(projection = projout)
ax.scatter(lon[:, :], lat[:, :], c=t2d, marker='.', transform=projin, s=0.5, cmap=mpl.colormaps['jet'])
feat = ax.add_feature(cfeature.LAND)
feat = ax.add_feature(cfeature.COASTLINE)
ax.set_extent([49, 55.25, -13.13, -9.45], crs=projin)
In the rheotaxis orientation method, the particles swim against the current. The user only provides a kappa parameter.
First, the angle of the current is computed as follows:
\[ \theta_{current} = \arctan2(V_{current}, U_{current}) \]
Then, the angle that the particle must follow is given by adding \(\pi\):
\[ \theta_{direction} = \theta_{current} + \pi \]
Finally, a random angle is picked up following a Von Mises distribution \(f(\theta, \theta_{direction}, \kappa_{reef})\)
(ref-rheo)=
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import os
from glob import glob
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
config = 'rheo'
if os.path.isdir('ichthyop_output'):
pattern = os.path.join('ichthyop_output', config, '*nc')
else :
pattern = os.path.join('..', '..', 'ichthyop_output', config, '*nc')
pattern
filelist = glob(pattern)
filelist
data = xr.open_mfdataset(filelist, decode_times=False)
data
# +
lon = data['lon'].values
lat = data['lat'].values
mort = data['mortality'].values
lon = np.ma.masked_where(mort > 0, lon)
lat = np.ma.masked_where(mort > 0, lat)
ntime, ndrifters = lon.shape
# +
time = np.arange(ntime)
drifters = np.arange(ndrifters)
d2d, t2d = np.meshgrid(drifters, time)
# +
plt.figure()
projin = ccrs.PlateCarree()
projout = ccrs.PlateCarree()
ax = plt.axes(projection = projout)
ax.scatter(lon[:, :], lat[:, :], c=t2d, marker='.', transform=projin, s=0.5, cmap=mpl.colormaps['jet'])
feat = ax.add_feature(cfeature.LAND)
feat = ax.add_feature(cfeature.COASTLINE)
ax.set_extent([49, 55.25, -13.13, -9.45], crs=projin)
In the reef orientation method, the larva will target the closest target area (for instance coral reef). These areas are defined in an XML zone file by a polygon and a zone-specific \(\kappa\) parameter. The user also provides the sensory detection threshold of the larva (maximum detection distance \(\beta\)).
If the distance between the particle and the barycenter of the closest reef (\(D\)) is below the detection thereshold \(\beta\), the larva will swim in the direction of the reef.
(ref-orientation)=
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 15
plt.figure(figsize = (10, 10))
xnew = 0
ynew = 0
#draw point at orgin
plt.plot(xnew, ynew, color = 'red', marker = 'o')
plt.gca().annotate('$P_{t}$', xy=(0 + 0.05, 0 - 0.07), xycoords='data', color='red')
#draw circle
r = 1.5
angles = np.linspace(0 * np.pi, 2 * np.pi, 100 )
xs = r * np.cos(angles)
ys = r * np.sin(angles)
plt.plot(xs, ys, color = 'k', ls='--', lw=0.5)
angle_old = np.pi / 6.
rold = 1.
xold = rold * np.cos(angle_old)
yold = rold * np.sin(angle_old)
plt.plot(xold, yold, marker='o', color='blue')
plt.gca().annotate('$P_{t - 1}$', xy=(xold + 0.05, yold), xycoords='data', color='blue')
angle_reef = np.pi + np.pi / 3.
print(np.rad2deg(angle_reef))
rreef = 1
xreef = rreef * np.cos(angle_reef)
yreef = rreef * np.sin(angle_reef)
print(yreef)
plt.plot(xreef, yreef, marker='o', color='orange')
plt.gca().annotate('$R$', xy=(xreef + 0.1, yreef), xycoords='data', color='orange')
plt.gca().annotate(r'$D$', xy=(0.5*(xreef) + 0.02, 0.5*yreef), xycoords='data', color='orange', ha='left')
plt.axvline(xnew, color='k', ls='--', lw=0.5)
plt.axhline(ynew, color='k', ls='--', lw=0.5)
p = np.polyfit([xold, xnew], [yold, ynew], deg=1)
xtemp = np.linspace(xold, xold -2, 100)
plt.plot(xtemp, np.polyval(p, xtemp), color='black', ls='--')
angle_current = angle_old + np.pi
tmp_angle = np.linspace(0, angle_current, 100)
rtmp = 0.2
plt.plot(0 + rtmp * np.cos(tmp_angle), 0 + rtmp * np.sin(tmp_angle), color='k', ls='--')
plt.gca().annotate(r'$\theta_{actual}$', xy=(xnew + 0.2, ynew+0.05), xycoords='data', color='k')
plt.plot([xnew, xreef], [ynew, yreef], color='orange', ls='--')
angle_current = angle_reef
tmp_angle = np.linspace(0, angle_current, 100)
rtmp = 0.3
plt.plot(0 + rtmp * np.cos(tmp_angle), 0 + rtmp * np.sin(tmp_angle), color='orange', ls='--')
plt.gca().annotate(r'$\theta_{reef}$', xy=(xnew - 0.3, ynew+0.1), xycoords='data', color='orange',ha='right')
angle_current = angle_old + np.pi
tmp_angle = np.linspace(angle_current, angle_reef, 100)
rtmp = 0.4
plt.plot(0 + rtmp * np.cos(tmp_angle), 0 + rtmp * np.sin(tmp_angle), color='plum', ls='--')
plt.gca().annotate(r'$\theta_{turning}$', xy=(-0.43, -0.34), xycoords='data', color='plum')
off = 0.05
plt.xlim(xreef - off, xold + off)
plt.ylim(yreef - off, yold + off)
plt.gca().set_aspect('equal')
plt.axis('off')
plt.show()239.99999999999997
-0.8660254037844384
First, the angle of the current trajectory, \(\theta_{actual}\), is computed by using the particle position at the previous time step (blue point) and the actual position (red point).
\[ \Delta_X = (X_{t - 1} - X_{t}) \]
\[ \Delta_Y = (Y_{t - 1} - Y_{t}) \]
\[ \theta_{actual} = \arctan2(\Delta Y, \Delta X) + \pi \]
The direction toward the reef, \(\theta_{reef}\) is also computed.
\[ \Delta_X = (X_{reef} - X_{t}) \]
\[ \Delta_Y = (Y_{reef} - Y_{t}) \]
\[ \theta_{reef} = \arctan2(\Delta Y, \Delta X) \]
The angles are computed in the \((X, Y)\) space. Therefore, longitude and latitude coordinates are converted in \((X, Y)\) using the latlon2xy Dataset methods.
The turning angle \(\theta_{turning}\) is given by:
\[ \theta_{turning} = \theta_{reef} - \theta_{actual} \]
The turning angle is then ponderated by the ratio of the distance from the reef to the detection threshold as follows:
\[ \theta_{ponderated} = \left(1 - \dfrac{D}{\beta}\right) \theta_{turning} \]
\[ \theta_{ponderated} = \left(1 - \dfrac{D}{\beta}\right) \left(\theta_{reef} - \theta_{actual}\right) \]
Therefore, the closest to the reef, the strongest the turning angle.
Then, a random angle is picked up following a Von Mises distribution \(f(\theta, \theta_{ponderated}, \kappa_{reef})\)
An example of a trajectory is provided below. In this case, two target destinations are provided (black boxes). The same \(\kappa\) value was used for both ares (1.2) and the \(\beta\) parameter has been set equal to 3 km.
(ref-orientation-2)=
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import os
from glob import glob
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
config = 'reef'
if os.path.isdir('ichthyop_output'):
pattern = os.path.join('ichthyop_output', config, '*nc')
else :
pattern = os.path.join('..', '..', 'ichthyop_output', config, '*nc')
pattern
filelist = glob(pattern)
filelist
data = xr.open_mfdataset(filelist, decode_times=False)
data
# +
lon = data['lon'].values
lat = data['lat'].values
mort = data['mortality'].values
lon = np.ma.masked_where(mort > 0, lon)
lat = np.ma.masked_where(mort > 0, lat)
ntime, ndrifters = lon.shape
# +
time = np.arange(ntime)
drifters = np.arange(ndrifters)
d2d, t2d = np.meshgrid(drifters, time)
# +
plt.figure()
projin = ccrs.PlateCarree()
projout = ccrs.PlateCarree()
ax = plt.axes(projection = projout)
ax.scatter(lon[:, :], lat[:, :], c=t2d, marker='.', transform=projin, s=0.5, cmap=mpl.colormaps['jet'])
feat = ax.add_feature(cfeature.LAND)
feat = ax.add_feature(cfeature.COASTLINE)
xp = np.array([51.625, 52.625, 52.625,51.625, 51.625])
yp = np.array([-1.179878E1, -1.179878E1, -1.079878E1, -1.079878E1, -1.179878E1])
ax.plot(xp, yp, transform=projin, color='k', ls='--')
xp = np.array([51.625, 52.625 - 0.5, 52.625 - 0.5,51.625, 51.625]) - 1
yp = np.array([-1.179878E1 + 0.5, -1.179878E1 + 0.5, -1.079878E1, -1.079878E1, -1.179878E1 + 0.5]) - 0.25
ax.plot(xp, yp, transform=projin, color='k', ls='--')
ax.set_extent([49, 55.25, -13.13, -9.45], crs=projin)