Commit ddc2f393 authored by nmenard's avatar nmenard
Browse files

Ajout convertisseur eco intégration vers SonarNetCDF

parent 40e5748b
......@@ -10,6 +10,7 @@ from pymovies_3d.core.echointegration.sampleechointegration import (
sample_echo_integration,
)
from pymovies_3d.core.echointegration.ei_bind import ei_bind
from pymovies_3d.core.export.hacei2xsf import hacei2xsf
......@@ -258,6 +259,34 @@ def ei_hacfile(
shutil.rmtree(chemin_save)
def ei_hacfile_netcdf(
path_hac_survey,
hacfilename,
path_config,
path_save) :
"""
Process EI over specified hac file and specified echosounders (ME70, EK80 and horizontal EK80)
inputs:
- path_hac_survey: hac directory
- hacfilename: name of hac file
- path_config: path of M3D configuration (with ESU definition, layers and EI threshold)
- path_save: path for resulting EI file
output:
- one SonarNetcdf file is created in path_save directory for each specified echosounde
"""
file_ini=path_hac_survey + "/" + hacfilename
chemin_save=path_save + "/" + hacfilename[:-4]
if not os.path.exists(chemin_save):
os.makedirs(chemin_save)
ei = hacei2xsf()
ei.savePath = chemin_save
ei.run(file_ini, path_config)
def ei_survey_multi_threshold(
path_hac_survey,
path_config,
......
#!/usr/bin/python
# -*- coding: utf-8 -*-
import netCDF4 as nc
import math
import os
import calendar
import datetime
import pymovies.pyMovies as mv
import pymovies_3d.core.export.xsfutils as xsfutils
from sonar_netcdf import sonar_groups as xsf
class hacei2xsf:
def __init__(self):
self.savePath = ""
self.datasets = {}
self.gridGroups = {}
# 0 : Sv, 1 : Sa
self.backscatter_type = 0
def __del__(self):
for ds in self.datasets.values():
ds.close()
def process_sounder(self, ds, sounder, nb_beams, nb_layersByType, nb_freqs):
rootGroup = xsf.RootGrp()
rootGroup.create_group(parent_group=ds)
sonarGroup = xsf.SonarGrp()
sonarGroup_ds = sonarGroup.create_group(parent_group=ds)
sonarGroup.sonar_serial_number = "M3D Sonar"
self.gridGroups[sounder.m_SounderId] = {}
gridGroupIdx = 1
for layerType, count in nb_layersByType.items() :
if count > 0:
print("Create grip group for type " + str(layerType) + ", with " + str(count) + " layer")
gridGroup = xsfutils.GridGroup(
parentGroup=sonarGroup_ds,
index = gridGroupIdx,
beamDimensions=nb_beams,
frequencyDimensions=nb_freqs,
rangeDimensions=count
)
gridGroup.backscatter_type = self.backscatter_type
#Range = 0, Depth = 1
gridGroup.range_axis_interval_type = 1 if layerType == mv.LayerType.SurfaceLayer else 0
self.gridGroups[sounder.m_SounderId][layerType] = gridGroup
gridGroupIdx += 1
def run(
self,
chemin_ini,
chemin_config,
nameTransect=None,
dateStart=None,
timeStart=None,
dateEnd=None,
timeEnd=None,
):
"""
Effectue l'écho-intégration du fichier ou du répertoire indiqué, et l'enregistre en .pickle
entrées:
- chemin_ini: chemin du répertoire ou du fichier à traiter
- chemin_config: chemin de la configuration M3D à utiliser (ESU, couches etc)
- nameTransect: nom de sauvegarde (optionnel). Si None, les fichiers pickles seront
nommés results_000X_EI.pickle
- dateStart et dateEnd: dates de début et fin des séquences à écho-intégrer, au format 'dd/mm/yyyy' ou ['dd/mm/yyyy','dd/mm/yyyy'...]
Si dateStart est None, tout le répertoire est écho-intégré
- timeStart et timeEnd: heures de début et fin des séquences à écho-intégrer, au format 'hh:mm:ss' ou ['hh:mm:ss','hh:mm:ss'...]
sortie:
- Un fichier pickle par chunk lu est créé
"""
# load configuration
if not os.path.exists(chemin_config):
print("****** stop: chemin de configuration non valide")
return
mv.moLoadConfig(chemin_config)
num_bloc = 1
if nameTransect is not None:
str_Name = "_" + nameTransect + "_"
else:
str_Name = ""
if dateStart is None:
range_date = [-1]
process_all = True
else:
if type(dateStart) == str:
dateStart = [dateStart]
if type(dateEnd) == str:
dateEnd = [dateEnd]
if type(timeStart) == str:
timeStart = [timeStart]
if type(timeEnd) == str:
timeEnd = [timeEnd]
range_date = range(len(dateStart))
if (len(dateEnd) + len(timeStart) + len(timeEnd)) / 3 != len(dateStart):
print("****** stop: les listes des dates/heures de début/fin n ont pas la même taille")
return
process_all = False
for x in range_date:
slice_end = False
if not process_all:
goto = calendar.timegm(
time.strptime("%s %s" % ((dateStart[x]), (timeStart[x])), "%d/%m/%Y %H:%M:%S")
)
mv.moOpenHac(chemin_ini)
util.hac_goto(goto)
timeEndEI = calendar.timegm(
time.strptime("%s %s" % (dateEnd[x], timeEnd[x]), "%d/%m/%Y %H:%M:%S")
)
else:
# chargement de tous les fichiers HAC contenus dans le repertoire
# et lecture du premier chunk
print("chemin_ini : " + str(chemin_ini))
mv.moOpenHac(chemin_ini)
timeEndEI = calendar.timegm(datetime.datetime.now().timetuple())
# Activate 'EchoIntegrModule'
EchoIntegrModule = mv.moEchoIntegration()
EchoIntegrModule.setEnable(True)
# récupération de la définition des couches
LayersDef = EchoIntegrModule.GetEchoIntegrationParameter().m_tabLayerDef
nbLayers = len(LayersDef)
nbLayersByType = {}
layerIndexes = {}
for indexLayer in range(nbLayers):
layerType = LayersDef[indexLayer].m_layerType
indexes = layerIndexes.get(layerType, [])
indexes.append(indexLayer)
layerIndexes[layerType] = indexes
nbLayersByType[layerType] = nbLayersByType.get(layerType, 0) + 1
# recupération du esu manager
# Time_seconds = 0, Distance_nautical_miles = 1, Distance_meters = 2, Number_of_ping = 3
esuMgr = mv.moGetEsuManager()
esuParamater = esuMgr.GetEsuManagerParameter()
if esuParamater.m_esuCutType == mv.ESUCutType.EsuCutByTime :
ping_axis_interval_type = 0
ping_axis_interval_value = esuParamater.m_time
elif esuParamater.m_esuCutType == mv.ESUCutType.EsuCutByDistance :
ping_axis_interval_type = 1
ping_axis_interval_value = esuParamater.m_distance
elif esuParamater.m_esuCutType == mv.ESUCutType.EsuCutByPingNb :
ping_axis_interval_type = 3
ping_axis_interval_value = esuParamater.m_pingNumber
FileStatus = mv.moGetFileStatus()
while not FileStatus.m_StreamClosed:
mv.moReadChunk()
AllData = EchoIntegrModule.GetOutputList()
# récupération des volumes de chaque couche (pour 1 ping)
LayersVolumes = EchoIntegrModule.GetLayersVolumes()
nbSounder = len(LayersVolumes[0])
# nb d'ESUs traitees
nbResults = len(AllData)
if nbResults > 0:
print("nb results : " + str(nbResults))
for indexResults in range(nbResults):
ESUResult = AllData[indexResults]
nbBeams = len(ESUResult.m_tabChannelResult)
print("EchoIntegrationOutput for sounder : " + str(ESUResult.m_SounderId))
timeESU = (ESUResult.m_timeEnd.m_TimeCpu + ESUResult.m_timeEnd.m_TimeFraction / 10000)
# On récupère le sondeur
sounder = self.get_sounder(ESUResult.m_SounderId)
# de la même façon créer un tableau de SoftChannels correspondant à l'index du beam (nbeams = Somme des softChan des transducteurs)
allChans = []
allTrans = []
for indexTrans in range(sounder.m_numberOfTransducer):
transducer = sounder.GetTransducer(indexTrans)
# indexBeam = beam in SonarNetCDF
for indexChan in range(transducer.m_numberOfSoftChannel):
allChans.append(transducer.getSoftChannelPolarX(indexChan))
allTrans.append(transducer)
nbChans = len(allChans)
# On récupère la liste de toutes les fréquences et les noms des transducteurs associés, pour tous les transducteurs
allFreqs = []
beamRefs = []
for esuResultIdx in range(len(ESUResult.m_tabChannelResult)):
esuResult = ESUResult.m_tabChannelResult[esuResultIdx]
for freqIdx in range(len(esuResult.m_tabFrequencies)):
allFreqs.append(esuResult.m_tabFrequencies[freqIdx])
beamRefs.append(allTrans[esuResultIdx].m_transName)
nbFreqs = len(allFreqs)
if ESUResult.m_SounderId not in self.gridGroups:
ds = nc.Dataset(
self.savePath
+ "_"
+ str(ESUResult.m_SounderId)
+ ".xsf.nc",
"w",
)
self.datasets[ESUResult.m_SounderId] = ds
self.process_sounder(ds, sounder, nbBeams, nbLayersByType, nbFreqs)
gridGroups = self.gridGroups[ESUResult.m_SounderId]
for gridType, gridGroup in gridGroups.items():
gridGroup.ping_axis_interval_type = ping_axis_interval_type
gridGroup.ping_axis_interval_value = ping_axis_interval_value
# On remplit les fréquences
for indexFreq in range(nbFreqs):
gridGroup.frequency[indexFreq] = allFreqs[indexFreq]
gridGroup.beam_reference[indexFreq] = beamRefs[indexFreq]
for indexChan in range(nbChans):
gridGroup.beam_type[indexChan] = allChans[indexChan].m_beamType
gridGroup.beam[indexChan] = allTrans[indexChan].m_transName
# On récupère le Grid Group pour le remplir
gridGroups = self.gridGroups[ESUResult.m_SounderId]
for gridType, gridGroup in gridGroups.items():
# index Ping, méthode plus simple
indexPing = len(gridGroup.sound_speed_at_transducer)
# Vitesse du son
gridGroup.sound_speed_at_transducer[indexPing] = sounder.m_soundVelocity
# stabilisation ? dans le contexte M3D => 1 si ME70, 0 sinon
gridGroup.beam_stabilisation[indexPing] = (1 if sounder.m_isMultiBeam else 0)
# 0 dans le contexte M3D
gridGroup.non_quantitative_processing[indexPing] = 0
# gridGroup.tx_transducer_depth[
# indexPing
# ] = mainBeamData.m_compensateHeave
indexAllFreq = 0
for indexBeam in range(nbBeams):
tabChannelResult = ESUResult.m_tabChannelResult[indexBeam]
freqs = tabChannelResult.m_tabFrequencies
gridGroup.cell_ping_time[indexPing, indexBeam] = timeESU
# On remplit le Bottom Range
gridGroup.detected_bottom_range[indexPing, indexBeam] = tabChannelResult.m_bottomDepth
transducer = allTrans[indexBeam]
softChan = allChans[indexBeam]
gridGroup.sample_interval[indexPing, indexBeam] = transducer.m_timeSampleInterval * 1e-64
gridGroup.blanking_interval[indexPing, indexBeam] = 0
gridGroup.transmit_power[indexPing, indexBeam] = softChan.m_transmissionPower
gridGroup.beamwidth_transmit_major[indexPing, indexBeam] = math.degrees(softChan.m_beam3dBWidthAthwartRad)
gridGroup.beamwidth_transmit_minor[indexPing, indexBeam] = math.degrees(softChan.m_beam3dBWidthAlongRad)
gridGroup.beamwidth_receive_major[indexPing, indexBeam] = math.degrees(softChan.m_beam3dBWidthAthwartRad)
gridGroup.beamwidth_receive_minor[indexPing, indexBeam] = math.degrees(softChan.m_beam3dBWidthAlongRad)
gridGroup.gain_correction[indexPing, indexBeam] = softChan.m_beamSACorrection
gridGroup.equivalent_beam_angle[indexPing, indexBeam] = math.pow(10, softChan.m_beamEquTwoWayAngle * 0.1)
# Ecriture des données dépendante du txBeam
# CW => CW(0), FM => LFM(1), pas HFM(2) dans M3D
gridGroup.transmit_type[indexPing, indexBeam] = 1 if transducer.m_pulseShape == 2 else 0
# => internalDelayTranslation => modify transducteur, on laisse à 0
gridGroup.sample_time_offset[indexPing, indexBeam] = 0
gridGroup.transducer_gain[indexPing, indexBeam] = softChan.m_beamGain
gridGroup.transmit_bandwidth[indexPing, indexBeam] = softChan.m_bandWidth
gridGroup.tx_beam_rotation_phi[indexPing, indexBeam] = math.degrees(softChan.m_mainBeamAthwartSteeringAngleRad)
gridGroup.rx_beam_rotation_phi[indexPing, indexBeam] = math.degrees(softChan.m_mainBeamAthwartSteeringAngleRad)
gridGroup.tx_beam_rotation_theta[indexPing, indexBeam] = math.degrees(softChan.m_mainBeamAlongSteeringAngleRad)
gridGroup.rx_beam_rotation_theta[indexPing, indexBeam] = math.degrees(softChan.m_mainBeamAlongSteeringAngleRad)
gridGroup.tx_beam_rotation_psi[indexPing, indexBeam] = 0
gridGroup.rx_beam_rotation_psi[indexPing, indexBeam] = 0
# M - Nominal duration of the transmit pulse. This is not the effective pulse duration.
gridGroup.transmit_duration_nominal[indexPing, indexBeam] = transducer.m_pulseDuration
# MA - Effective duration of the received pulse. This is the duration of the square pulse containing the same energy as the actual receive pulse. This parameter is either theoretical or comes from a calibration exercise and adjusts the nominal duration of the transmitted pulse to the measured one. During calibration it is obtained by integrating the energy of the received signal on the calibration target normalised by its maximum energy. Necessary for type 1, 2,3 and 4 conversion equations.
if transducer.m_pulseShape == 2:
T = softChan.m_softChannelComputeData.GetTransmitSignalObject(transducer, softChan).m_effectivePulseLength
gridGroup.transmit_frequency_start[indexPing, indexBeam] = softChan.m_startFrequency
gridGroup.transmit_frequency_stop[indexPing, indexBeam] = softChan.m_endFrequency
else:
T = transducer.m_pulseDuration / 1000000.0
gridGroup.transmit_frequency_start[indexPing, indexBeam] = softChan.m_acousticFrequency
gridGroup.transmit_frequency_stop[indexPing, indexBeam] = softChan.m_acousticFrequency
effectivePulseLength = T * math.pow(10, (2.0 * softChan.m_beamSACorrection) / 10.0)
gridGroup.receive_duration_effective[indexPing, indexBeam] = effectivePulseLength
# Cell Latitude/Longitude/Depth
for idx, indexLayer in enumerate(layerIndexes.get(gridType, [])) :
layerResult = tabChannelResult.m_tabLayerResult[0][indexLayer]
gridGroup.cell_latitude[indexPing, idx, indexBeam] = layerResult.m_latitude
gridGroup.cell_longitude[indexPing, idx, indexBeam] = layerResult.m_longitude
gridGroup.cell_depth[idx, indexBeam] = layerResult.m_Depth
for indexFreq in range(len(freqs)):
freqResult = tabChannelResult.m_tabLayerResult[indexFreq]
for idx, indexLayer in enumerate(layerIndexes.get(gridType, [])) :
# Sa o Sv
gridGroup.integrated_backscatter[indexPing, idx, indexAllFreq] = (
freqResult[indexLayer].m_Sv
if self.backscatter_type == 0 else
freqResult[indexLayer].m_Sa
)
indexAllFreq += 1
if (timeESU > timeEndEI) & (not process_all):
slice_end = True
break
# Suppression des résultats d'EI consommés
EchoIntegrModule.ReleaseOutputList()
num_bloc = num_bloc + 1
if slice_end:
break
EchoIntegrModule.setEnable(False)
return
def get_sounder(self, sounder_id):
sounders = mv.moGetSounderDefinition()
nb_sounders = sounders.GetNbSounder()
for i in range(nb_sounders):
sounder = sounders.GetSounder(i)
if sounder.m_SounderId == sounder_id:
return sounder
\ No newline at end of file
from sonar_netcdf import sonar_groups as xsf
from sonar_netcdf.nc_tools import find_type
import netCDF4 as nc
import numpy as np
import enum
class BeamType(enum.IntEnum):
BeamTypeSplit = 1
BeamTypeSplit3 = 17
BeamTypeSplit3CN = 65
class NcGroup:
def create_attr(self, attr):
fname = "create_" + attr
f = getattr(self.xsfGroup, fname, None)
if f is not None:
setattr(self, attr, f(self.g))
else:
print(fname + " not found !")
def create_attributes(self, ncType, exclude = []):
attributeList = [func for func in dir(ncType) if callable(getattr(ncType, func)) and func.startswith("create_")]
attributeList.remove("create_group")
attributeList.remove("create_dimension")
for a in exclude:
fname = "create_" + a
if fname in attributeList:
attributeList.remove(fname)
for fname in attributeList:
attr = fname[7:]
self.create_attr(attr)
class SoundSpeedProfileGroup(NcGroup):
def __init__(self, parentGroup):
self.xsfGroup = xsf.SoundSpeedProfileGrp()
self.g = self.xsfGroup.create_group(parent_group=parentGroup)
soundSpeedProfileDimensions = {}
soundSpeedProfileDimensions[xsf.SoundSpeedProfileGrp.PROFILE_DIM_NAME] = None
soundSpeedProfileDimensions[xsf.SoundSpeedProfileGrp.SAMPLECOUNT_DIM_NAME] = None
self.xsfGroup.create_dimension(group=self.g, values=soundSpeedProfileDimensions)
self.create_attributes(ncType=xsf.SoundSpeedProfileGrp)
class AttitudeGroup(NcGroup):
def __init__(self, parentGroup, ident):
self.xsfGroup = xsf.AttitudeSubGroup()
self.g = self.xsfGroup.create_group(parent_group=parentGroup, ident=ident)
dimensions = {}
dimensions[xsf.AttitudeSubGroup.TIME_DIM_NAME] = None
self.xsfGroup.create_dimension(group=self.g, values=dimensions)
self.create_attributes(ncType=xsf.AttitudeSubGroup)
class PositionGroup(NcGroup):
def __init__(self, parentGroup):
self.xsfGroup = xsf.PositionSubGroup()
self.g = self.xsfGroup.create_group(parent_group=parentGroup, ident="1")
dimensions = {}
dimensions[xsf.PositionSubGroup.TIME_DIM_NAME] = None
self.xsfGroup.create_dimension(group=self.g, values=dimensions)
self.create_attributes(ncType=xsf.PositionSubGroup)
class PlatformGroup(NcGroup):
def __init__(self, parentGroup, sounder):
self.xsfGroup = xsf.PlatformGrp()
self.g = self.xsfGroup.create_group(parent_group=parentGroup)
dimensions = {}
dimensions[xsf.PlatformGrp.TRANSDUCER_DIM_NAME] = sounder.m_numberOfTransducer
dimensions[xsf.PlatformGrp.POSITION_DIM_NAME] = 1
dimensions[xsf.PlatformGrp.MRU_DIM_NAME] = 1
self.xsfGroup.create_dimension(group=self.g, values=dimensions)
self.create_attributes(ncType=xsf.PlatformGrp)
# create attitude and position sub groups
#attitudeGroup = xsf.AttitudeGrp().create_group(parent_group=parentGroup)
self.attitudes = {}
positionGroup = xsf.PositionGrp().create_group(parent_group=parentGroup)
self.position = PositionGroup(parentGroup=positionGroup)
#self.mru_ids[0] = 0
self.mru_offset_x[0] = 0
self.mru_offset_y[0] = 0
self.mru_offset_z[0] = 0
self.mru_rotation_x[0] = 0
self.mru_rotation_y[0] = 0
self.mru_rotation_z[0] = 0
#self.position_ids[0] = 0
self.position_offset_x[0] = 0
self.position_offset_y[0] = 0
self.position_offset_z[0] = 0
# position des transducteurs
for transducerIdx in range(0, sounder.m_numberOfTransducer):
transducer = sounder.GetTransducer(transducerIdx)
self.transducer_function[transducerIdx] = 3
self.transducer_ids[transducerIdx] = transducer.m_transName
platform = transducer.GetPlatform()
self.transducer_offset_x[transducerIdx] = platform.along_ship_offset
self.transducer_offset_y[transducerIdx] = platform.athwart_ship_offset
self.transducer_offset_z[transducerIdx] = platform.depth_offset
self.transducer_rotation_x[transducerIdx] = math.degrees(transducer.m_transFaceAlongAngleOffsetRad)
self.transducer_rotation_y[transducerIdx] = math.degrees(transducer.m_transFaceAthwarAngleOffsetRad)
self.transducer_rotation_z[transducerIdx] = math.degrees(transducer.m_transRotationAngleRad)
class BeamGroup(NcGroup):
def __init__(self, parentGroup, index, numberOfSoftChannel, numberOfQuadrants, isSv):
sample_t_type = np.int16 if numberOfQuadrants == 1 else np.float32
angle_t_type = np.int16
vlen_type_dict = {
"sample_t" : sample_t_type,
"angle_t" : angle_t_type
}
self.xsfGroup = xsf.BeamGroup1Grp()
self.g = self.xsfGroup.create_group(parent_group=parentGroup, vlen_type_dict=vlen_type_dict, ident="Beam_group" + str(index))
beamDimensions = {}
beamDimensions[xsf.BeamGroup1Grp.BEAM_DIM_NAME] = numberOfSoftChannel
beamDimensions[xsf.BeamGroup1Grp.SUBBEAM_DIM_NAME] = numberOfQuadrants
beamDimensions[xsf.BeamGroup1Grp.PING_TIME_DIM_NAME] = None
beamDimensions[xsf.BeamGroup1Grp.TX_BEAM_DIM_NAME] = numberOfSoftChannel
self.xsfGroup.create_dimension(group=self.g, values=beamDimensions)
self.backscatter_r = self.xsfGroup.create_backscatter_r(group=self.g)
self.echoangle_major = self.xsfGroup.create_echoangle_major(group=self.g)
self.echoangle_minor = self.xsfGroup.create_echoangle_minor(group=self.g)
if numberOfQuadrants == 1: # Sv of power
if sample_t_type == np.int16 :
self.backscatter_r.missing_value = -32768
if isSv :
self.backscatter_r.scale_factor = np.float32(0.01)
if angle_t_type == np.int16 :
self.echoangle_major.missing_value = -32768
self.echoangle_major.scale_factor = np.float32(0.01)
self.echoangle_major.valid_range = np.int16([-18000, 18000])
self.echoangle_minor.missing_value = -32768
self.echoangle_minor.scale_factor = np.float32(0.01)
self.echoangle_minor.valid_range = np.int16([-18000, 18000])
self.create_attributes(ncType=xsf.BeamGroup1Grp, exclude=["backscatter_r", "echoangle_major", "echoangle_minor"])
# Disable automatic conversion
self.backscatter_r.set_auto_maskandscale(False)
self.backscatter_i.set_auto_maskandscale(False)
self.echoangle_major.set_auto_maskandscale(False)
self.echoangle_minor.set_auto_maskandscale(False)
class GridGroup(NcGroup):
def __init__(
self,
parentGroup,
index,
beamDimensions=None,
frequencyDimensions=None,
rangeDimensions=None,
):
# sample_t_type = np.int16 if numberOfQuadrants == 1 else np.float32
sample_t_type = np.int16
angle_t_type = np