Manipulate projections (Halieutic dataset)¶
In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.proj as prj
import gstlearn.document as gdoc
import matplotlib.pyplot as plt
import numpy as np
The ICES Working Group WGACEGG dataset: Acoustic Noise Scale (NASC) densities of anchovies and sardines per 1 nautical mile sampling unit, observed in spring 2018 and 2021.
Column descriptions:
- survey: campaign
- year: year
- time: date/time of data recording x: longitude in decimal degrees (WGS84)
- y: latitude in decimal degrees (WGS84)
- NASC: acoustic density in m².MN⁻²
- sp: species (anchovy = "ENGR-ENC"
- sardine = "SARD-PIL")
For more information on the origin of these data: Doray, M., Van Der Kooij, J., Boyra, G. (Eds.), 2021. ICES Survey Protocols - Manual for acoustic surveys coordinated under the ICES Working Group on Acoustic and Egg Surveys for Small Pelagic Fish (WGACEGG). https://doi.org/10.17895/ICES.PUB.7462.
Load the data in long/lat (alias WGS84 or EPSG:4326)
In [2]:
# Read the data
csv = csv = gl.CSVformat.create(True, charSep=";")
file_csv = gdoc.loadData("halieutic", "AC-SPRING-IBBB-NASC_ANE-PIL_2018-2021.csv")
data = gl.Db.createFromCSV(file_csv, csv=csv)
data.setLocators(["x", "y"], gl.ELoc.X)
data.setLocator("NASC", gl.ELoc.Z)
# Read the polygon
file_csv = gdoc.loadData("halieutic", "WGACEGGspringPolygon.csv")
poly = gl.Polygons.createFromCSV(file_csv, csv)
print("Polygon extension:", poly.getExtensionAsVD(), " degrees")
# Read the boundaries
name = gdoc.loadData("boundaries", "world.poly")
world = gl.Polygons.createFromNF(name)
print("Boundaries extension:", np.round(world.getExtensionAsVD(), 0), " degrees")
# Plot
fig, ax = gp.init(flagEqual=True)
ax.set_xlim([-15, 10])
ax.set_ylim([35, 52])
ax.symbol(data, s=1)
ax.polygon(poly)
ax.polygon(
world, edgecolor="black", facecolor="lightblue", flagFace=True, flagLabels=True
)
plt.show()
Polygon extension: [-9.94955409 -1.23270531 35.99409516 48.02438771] degrees Boundaries extension: [-180. 180. -90. 84.] degrees
Transform all coordinates to EPSG:2154 and display only some countries
In [3]:
crsFrom = "EPSG:4326" # WGS84 (lat/lon)
crsTo = "EPSG:2154"
# Transform the data coordinates
datat = data.clone()
x, y = prj.proj(datat["x"], datat["y"], crsFrom, crsTo)
datat["x"] = x
datat["y"] = y
# Transform polygon coordinates
polyt = poly.clone()
x, y = prj.proj(polyt.getX(0), polyt.getY(0), crsFrom, crsTo)
polyt.setX(0, x)
polyt.setY(0, y)
ext = polyt.getExtensionAsVD()
print("Polygon extension:", np.round(ext, 0), " meters")
# Transform boundaries coordinates (only for France, Portugal and Spain)
ids = [119, 214, 215]
worldt = gl.Polygons()
nid = 0
for id in ids:
worldt.addPolyElem(world.getPolyElem(id))
x, y = prj.proj(worldt.getX(nid), worldt.getY(nid), crsFrom, crsTo)
worldt.setX(nid, x)
worldt.setY(nid, y)
nid += 1
print("Boundaries extension:", np.round(worldt.getExtensionAsVD(), 0), " meters")
# Plot
fig, ax = gp.init(flagEqual=True)
ax.symbol(datat, s=1)
ax.polygon(polyt)
ax.polygon(
worldt,
edgecolor="black",
facecolor="lightblue",
flagFace=True,
flagLabels=True,
labels=["France", "Portugal", "Spain"],
)
plt.show()
Polygon extension: [-427232. 368841. 5479556. 6797888.] meters Boundaries extension: [-393135. 1072736. 5463375. 7117128.] meters
In [4]:
# Plot the two basemaps together from the initial data
fig, axs = gp.init(ny=2)
axs[0].baseMap(data, size=0.5)
axs[1].baseMap(data, size=0.5, flagProj=True, crsFrom=crsFrom, crsTo=crsTo)
plt.show()
In [ ]: