This Tutorial has been established from:

ICES Cooperative Research Report Handbook of Geostatistics in R for fisheries and marine ecology

Based on: ICES Training Course Application of Geostatistics to analyse spatially explicit Survey data in an Ecosystem Approach December 2013-2014, Fontainebleau, Centre de Geosciences, Mines ParisTech chapter: Indices of spatial distribution Author: Mathieu Woillez, Ifremer Date: 14/01/2015

Bez N., and Rivoirard J., 2000, Indices of collocation between populations. In: Checkley D.M., Hunter J.R., Motos L., von der Lingen C.D. (Eds.).Workshop on the Use of Continuous Underway Fish Egg Sampler (CUFES) for mapping spawning habitat of pelagic fish. GLOBEC Rep., pp. 48-52.

Loading the two files from the official depository. It consists in two files: - the Data file containing the captures on Bay of Biscay - the Polygon delineating the area of interest

The Data file contains a series of samples with three variables of interest: - the captures of age 0 (called A0) - the captures of age 1 (called A1) - the area of Influence attached to each datum.

db.data = Db_createFromNF(loadData("halieutic", "Hake_Bob_Db.NF"))
poly.data = Polygons_createFromNF(loadData("halieutic", "Hake_Bob_Polygon.NF"))

Visualizing the data set

ggplot() + 
  plot.polygon(poly.data, fill="NA") + 
  plot.point(db.data, nameSize="A0", colour="blue") +
  plot.decoration(title="Age-0 hake density")

We define a projection, based on the center of gravity of the data. Modify the coordinates of the Data Base and of the Polygon consequently.

projec = Projection(TRUE, db.data)

err = projec$operateOnDb(db.data)
err = projec$operateOnPolygons(poly.data)

Display in the projected system

ggplot() + 
  plot.polygon(poly.data, fill="NA") + 
  plot.point(db.data, nameSize="A0", color="blue") +
  plot.decoration(title="Age-0 hake density (in projected space)")

Instantiate the class for calculating the Spatial Indices.

cg = SpatialIndices(db.data)

Calculate the center of gravity, the Inertia and the Isotropy of the fish densities. The Center of Gravity and Inertia are calculated from the data location. If the argument ‘name’ is specified, the calculations take it into account. Finally a weight can be attached to each datum (if a variable is assigned to the locator “w”).

err = cg$computeCGI("A0")
cgAxesA0 = cg$getAxes()

Get the results of the calculation of the Center of gravity (the calculations are performed in the projected space). The results are back transformed into the initial system prior to displaying them.

vec = projec$operateInvert(cg$getCenter())
cat("Center of gravity = ",vec,"\n")
## Center of gravity =  -3.780361 47.09949
cat("Inertia =           ",cg$getInertia(),"\n")
## Inertia =            2773.657
cat("Isotropy =          ",cg$getIso(),"\n")
## Isotropy =           0.38866

Compute the center of gravity of the samples (without considering any specific variable).

err = cg$computeCGI("")
cgAxesP0 = cg$getAxes()

Get the results of the calculation of the Center of gravity (in the initial system)

vec = projec$operateInvert(cg$getCenter())
cat("Center of gravity = ",vec,"\n")
## Center of gravity =  -4.091006 46.74465
cat("Inertia =           ",cg$getInertia(),"\n")
## Inertia =            12474.91
cat("Isotropy =          ",cg$getIso(),"\n")
## Isotropy =           0.3059677

Display the sample within the Polygon, together with the two inertia calculated beforehand: the one using the sample locations only (in blue) and the one using the weighted variable A0 into account (in red).

ggplot() + 
  plot.polygon(poly.data, fill="NA") + 
  plot.point(db.data, color="black") +
  plot.decoration(title="Center of gravity and Inertia of densities and Samples") +
  plot.XY(cgAxesA0[[1]], cgAxesA0[[2]], linewidth=1, col="red") + 
  plot.XY(cgAxesA0[[3]], cgAxesA0[[4]], linewidth=1, col="red") +
  plot.XY(cgAxesP0[[1]], cgAxesP0[[2]], linewidth=1, col="blue") + 
  plot.XY(cgAxesP0[[3]], cgAxesP0[[4]], linewidth=1, col="blue")

The two populations (variables A0 and A1) are displayed simultaneously.

ggplot() + 
  plot.polygon(poly.data, fill="NA") + 
  plot.point(db.data, nameSize="A0", color="red") +
  plot.point(db.data, nameSize="A1", color="blue") +
  plot.decoration(title="Age-0 and age-1 hake densities") 
## Scale for size is already present.
## Adding another scale for size, which will replace the existing scale.

The Local overlapping between two populations is addressed using the local index of collocation, that is, the non-centered correlation between their densities.

lic = cg$getLIC("A0", "A1")
cat("Local index of collocation = ",lic,"\n")
## Local index of collocation =  0.6774079

Calculate the Global Index of Collocation which looks at how geographically distinct two populations are by comparing the distance between their centers of gravity and the mean distance between individual fish taken at random and independently from each population.

gic = cg$getGIC("A0", "A1")
cat("Global index of collocation = ",gic,"\n")
## Global index of collocation =  0.8977238

Calculate: - the Abundance index - the Positive area - the Spreading area - the Equivalent area

The spreading area ‘sparea’ (expressed in square nautical miles) is derived from the curve expressing (Q-Q(T))/Q as a function of T. T being the cumulated area occupied by the density values, ranked in decreasing order. Q(T) the correspondi ng cumulated abundance and Q the overall abundance. The spreading area is simply defined as twice the area below this curve.

err = cg$spatial("A0")
## Abundance Index = 69915207.365000
## Positive Area   = 22895.454300
## Equivalent Area = 4771.683336
VDD = cg$getQT("A0")
## Spreading Area  = 5664.107384
ggplot() + 
  plot.XY(VDD[[1]], VDD[[2]], flagPoint=TRUE) + 
  plot.decoration(title="Curve: (Q - Q(T)) / Q", xlab="T", ylab="(Q-Q(T))/Q")

Calculate the patches and represent them by color.

centers = cg$getPatches("A0", 100, 10)
## Regrouping the information of A0 by patches
## - Distance to Center of Gravity = 100.000000
## - Total Number of patches = 4
## - Number of patches with abundance >  10 % = 1
## - Percentage of abundance in these patches =  91.048
## - Percentage of area in these patches =  35.199
ggplot() + 
  plot.polygon(poly.data, fill="NA") + 
  plot.point(db.data, nameColor="Patch", palette=c("blue", "red", "green", "black"),
             asFactor=TRUE) +
  plot.decoration(title="Patches")
## Warning: Using size for a discrete variable is not advised.
## Using size for a discrete variable is not advised.

Calculate the Microstructure index which measures the relative importance of structural components that have a smaller scale than the sampling lag (including random noise). This calculation is not completed in the current version of the package gstlearn.

micro = cg$getMicroStructure("A0",10.,NULL,50,400)
cat("Microstructure Index = ",micro,"\n")