Commit e82385a1e9198f4fe9d52e401f6d0540292c1e82

Authored by Paul Nicot
1 parent 0cd511e2
Exists in master and in 2 other branches mingan, simon

Ajout des fichiers matlab tr_sort.m (compilation des positions par bouées et par mission) et utc.m

Ajout du script tr_postition.py (script python pour récupérer les positions spots en continue)
script à faire tourner sur Brandypot
Showing 3 changed files with 573 additions and 0 deletions   Show diff stats
tr_position.py 0 → 100644
... ... @@ -0,0 +1,408 @@
  1 +#!/usr/bin/env python
  2 +#! python
  3 +
  4 +import sys
  5 +import xml.dom.minidom as minidom
  6 +import urllib
  7 +import math
  8 +import csv
  9 +import fiona
  10 +import shapely.geometry
  11 +from time import localtime, strftime, ctime
  12 +
  13 +
  14 +
  15 +#
  16 +#Liste globale contenant tous les objets
  17 +#instancies a partir de la classe Spot
  18 +#
  19 +spots = []
  20 +
  21 +
  22 +class Spot:
  23 + '''Classe pour les GPS SPOT'''
  24 + spotCount=0
  25 + rayonTerre = 6371.0 # Rayon moyen de la terre
  26 +
  27 + distancePlusPres = 100000000.0 # Init a 1 million de km !
  28 + capPlusPres = 0.0
  29 + identPlusPres = 0
  30 + parentPlusPres = 0
  31 +
  32 + spotPlusPres = None
  33 + dejaTraite = False
  34 +
  35 + def __init__(self,esnName):
  36 + self.esnName = esnName
  37 + Spot.spotCount += 1
  38 +
  39 + def setEsn(self,esn):
  40 + self.esn = esn
  41 +
  42 + def setTimestamp(self, timestamp):
  43 + self.timestamp = timestamp
  44 +
  45 + def setTimeGMT(self, timeGMT):
  46 + self.timeGMT = int(timeGMT)
  47 +
  48 + def setTimeStampEST(self, timeGMT):
  49 + self.timeStampEST = ctime(int(timeGMT))
  50 +
  51 + def setLatitude(self,latitude):
  52 + self.latitude = float(latitude)
  53 +
  54 + def setLongitude(self, longitude):
  55 + self.longitude = float(longitude)
  56 +
  57 + def setMessageType(self, messagetype):
  58 + self.messagetype = messagetype
  59 +
  60 + def setFlag(self, flag):
  61 + self.flag = flag
  62 +
  63 +
  64 +
  65 +############################################################
  66 +#
  67 +#Calculer la distance en km entre deux points
  68 +#le long d'une ligne de vent (rhumb line)
  69 +#Le calcul se fait entre le spot courant (self) et le spot cible (ref)
  70 +#
  71 +#Ref : http://williams.best.vwh.net/avform.htm#Rhumb
  72 +#
  73 +############################################################
  74 +
  75 + def calcDistance(self, ref, plusPres):
  76 + if ref == self:
  77 + return 0
  78 +
  79 + if ref == None:
  80 + return 0
  81 +
  82 + if plusPres == True and ref.dejaTraite == True:
  83 + return 0
  84 +
  85 +
  86 + lat1 = math.radians(self.latitude)
  87 + lon1 = math.radians(self.longitude)
  88 + lat2 = math.radians(ref.latitude)
  89 + lon2 = math.radians(ref.longitude)
  90 + dlat = math.radians(ref.latitude - self.latitude)
  91 + dlon = math.fabs(math.radians(ref.longitude - self.longitude))
  92 +
  93 + dPhi = math.log(math.tan(lat2/2+math.pi/4)/math.tan(lat1/2+math.pi/4))
  94 +
  95 + if dPhi <> 0 :
  96 + q = dlat/dPhi
  97 + else :
  98 + q = math.cos(lat1)
  99 +
  100 +
  101 + #if dLon > 180 degres take shorter rhumb across anti-meridian:
  102 + if math.fabs(dlon) > math.pi :
  103 + if dlon > 0 :
  104 + dlon= -(2*math.pi - dlon)
  105 + else:
  106 + dlon = (2*math.pi + dlon)
  107 +
  108 + dist = math.sqrt(dlat*dlat + q*q*dlon*dlon) * Spot.rayonTerre
  109 +
  110 + return dist
  111 +
  112 +############################################################
  113 +#
  114 +#Calculer le cap (bearing) pour atteindre le spot cible (ref)
  115 +#a partir du spot courant (self)
  116 +#
  117 +############################################################
  118 +
  119 + def calcCap(self, ref):
  120 + lat1 = math.radians(self.latitude)
  121 + lat2 = math.radians(ref.latitude)
  122 + dlon = math.radians(ref.longitude - self.longitude)
  123 +
  124 + dPhi = math.log(math.tan(lat2/2+math.pi/4)/math.tan(lat1/2+math.pi/4))
  125 +
  126 +
  127 + if math.fabs(dlon) > math.pi :
  128 + if dlon > 0 :
  129 + dlon= -(2*math.pi - dlon)
  130 + else:
  131 + dlon = (2*math.pi + dlon)
  132 +
  133 + return (math.degrees(math.atan2(dlon,dPhi)) + 360.0) % 360
  134 +
  135 +
  136 + def printInfoPlusPres(self) :
  137 + if self.spotPlusPres == None:
  138 + return
  139 +
  140 + print 'Voisin de Spot {} est:'.format(self.esnName)
  141 + print 'Spot {} : Lat={} , Lon = {} , Dist = {} , Cap = {}{}'.format(self.identPlusPres, self.spotPlusPres.latitude,self.spotPlusPres.longitude,self.distancePlusPres, self.capPlusPres, '\n')
  142 +
  143 +
  144 +#############################################################
  145 +#
  146 +#Lecture du fichier XML contenant les information des spots et
  147 +#Instanciation de tous les objets Spot a partir de ces info.
  148 +#
  149 +#A la sortie de cette fonction, la liste globale spots[] contient tous
  150 +#tous les spots necessaires au traitement
  151 +#
  152 +#############################################################
  153 +
  154 +def getTags(url):
  155 + """
  156 + Initialiser tous les spots a partir du fichier XML
  157 + """
  158 +
  159 + doc = minidom.parse(urllib.urlopen(url))
  160 + node = doc.documentElement
  161 + messages = doc.getElementsByTagName("message")
  162 +
  163 +# Initialisation des listes de donnes
  164 +
  165 + esns = []
  166 + esnNames = []
  167 + timestamps = []
  168 + timeInGMTSeconds = []
  169 + latitudes = []
  170 + longitudes = []
  171 + messagetype = []
  172 + spottemp = []
  173 +
  174 +# On boucle sur tous les tags messages du fichier XML
  175 +# et on extrait les informations
  176 +
  177 + if len(messages) == 0:
  178 + sys.exit(1)
  179 +
  180 + for message in messages:
  181 + esnObj = message.getElementsByTagName("messengerId")[0]
  182 + esns.append(esnObj)
  183 + esnNameObj = message.getElementsByTagName("messengerName")[0]
  184 + esnNames.append(esnNameObj)
  185 + timestampObj = message.getElementsByTagName("dateTime")[0]
  186 + timestamps.append(timestampObj)
  187 + timeGMTObj = message.getElementsByTagName("unixTime")[0]
  188 + timeInGMTSeconds.append(timeGMTObj)
  189 + latObj = message.getElementsByTagName("latitude")[0]
  190 + latitudes.append(latObj)
  191 + lonObj = message.getElementsByTagName("longitude")[0]
  192 + longitudes.append(lonObj)
  193 + msgtypeObj = message.getElementsByTagName("messageType")[0]
  194 + messagetype.append(msgtypeObj)
  195 +
  196 +
  197 +#On instancie une liste temporaire d'objets de classe Spot
  198 +#qui contient l'information de tous les spots
  199 + for esnName in esnNames:
  200 + nodes = esnName.childNodes
  201 + for node in nodes:
  202 + spottemp.append(Spot(node.data))
  203 +
  204 +#On conserve l'info dans chaque objet
  205 +
  206 + for i,s in enumerate(spottemp):
  207 + nodes = esns[i].childNodes
  208 + for node in nodes:
  209 + s.setEsn(node.data)
  210 + nodes = timestamps[i].childNodes
  211 + for node in nodes:
  212 + s.setTimestamp(node.data)
  213 + nodes = timeInGMTSeconds[i].childNodes
  214 + for node in nodes:
  215 + s.setTimeGMT(node.data)
  216 + s.setTimeStampEST(node.data)
  217 + nodes = latitudes[i].childNodes
  218 + for node in nodes:
  219 + s.setLatitude(node.data)
  220 + nodes = longitudes[i].childNodes
  221 + for node in nodes:
  222 + s.setLongitude(node.data)
  223 + nodes = messagetype[i].childNodes
  224 + for node in nodes:
  225 + s.setMessageType(node.data)
  226 +
  227 +#
  228 +#On produit une liste unique des identifiants des spots
  229 +#
  230 + spotlist = sorted(set([s.esnName for s in spottemp ]))
  231 +
  232 +#
  233 +#On boucle sur la liste des identifiants uniques
  234 +#et pour chaque spot on conserve la donnee la plus recente
  235 +#(max timestamp)
  236 +#
  237 + for i in spotlist :
  238 + spotmax = None
  239 + for s in spottemp :
  240 + if s.esnName == i:
  241 + if spotmax == None:
  242 + spotmax = s
  243 + if int(s.timeGMT) > int(spotmax.timeGMT) :
  244 + spotmax = s
  245 +
  246 + spots.append(spotmax)
  247 +
  248 +#
  249 +#Rendu ici, la liste spots contient une entree par spot
  250 +#et cette entree est la plus recente (plus grand timestamp)
  251 +# Fin de la fonction
  252 +#
  253 +
  254 +############################################################
  255 +#
  256 +#Produire un fichier html decrivant la sequence de repechage
  257 +#La page se rafraichie automatiquement au 5 minutes sur le fureteur
  258 +#(voir balise META)
  259 +#
  260 +############################################################
  261 +
  262 +def ecrireHtml(fichier) :
  263 + f = open(fichier,"w")
  264 +
  265 + #
  266 + #Creation du tableau pour la sequence de repechage
  267 + #
  268 + f.write("<html><head><META HTTP-EQUIV=\"refresh\" CONTENT=\"60; URL=http://demeter.uqar.ca/spots/index.html\"></head><body>\n")
  269 + f.write("Derniere mise-a-jour : " + strftime("%a, %d %b %Y %X", localtime())+ "\n")
  270 + f.write("<table border=\"1\">\n")
  271 + f.write("<tr><th>D&eacute;part</th><th>Cible</th><th>Latitude</th><th>Longitude</th><th>Dist.(km)</th><th>Cap</th><th>Date-heure</th></tr>\n")
  272 +
  273 + s = spots[0]
  274 + while not s == None:
  275 + if not s.spotPlusPres == None :
  276 + f.write("<tr><td>" + str(s.esnName) + "</td><td>" + str(s.identPlusPres) + "</td><td>" + str(s.spotPlusPres.latitude) + "</td>" + \
  277 + "<td>" + str(s.spotPlusPres.longitude) + "</td><td>" + str("{0:.3f}".format(s.distancePlusPres)) + "</td>" + \
  278 + "<td>" + str("{0:.1f}".format(s.capPlusPres)) + "</td>" + \
  279 + "<td>" + s.spotPlusPres.timeStampEST + "</td></tr>" )
  280 +
  281 + s = s.spotPlusPres
  282 +
  283 + f.write("</table>\n")
  284 +
  285 + #
  286 + #Creation du tableau des distances par rapport au Spot 0
  287 + #
  288 + f.write("<p />Distances du Spot " + str(spots[0].esnName))
  289 + f.write("<table border=\"1\">\n")
  290 + f.write("<tr><th>Cible</th><th>Latitude</th><th>Longitude</th><th>Dist. (km)</th><th>Cap</th><th>Date-heure</th></tr>\n")
  291 +
  292 + s = spots[0]
  293 + for b in spots[1:] :
  294 + f.write("<tr><td>" + str(b.esnName) + "</td><td>" + str(b.latitude) + "</td>" + \
  295 + "<td>" + str(b.longitude) + "</td><td>" + str("{0:.2f}".format(s.calcDistance(b,False))) + "</td>" + \
  296 + "<td>" + str("{0:.1f}".format(s.calcCap(b))) + "</td>" + \
  297 + "<td>" + b.timeStampEST + "</td></tr>" )
  298 +
  299 + f.write("</table><p /><p/>\n")
  300 + f.write("<a href=\"./spots.kml\">spots.kml</a>\n")
  301 +
  302 +
  303 + f.write("</body></html>")
  304 + f.close()
  305 +
  306 +
  307 +############################################################
  308 +#
  309 +#Produire un fichier kml pour tous les spots
  310 +#
  311 +############################################################
  312 +
  313 +def ecrireKml(fichier) :
  314 + f = open(fichier,"w")
  315 +
  316 + f.write("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
  317 + f.write("<kml xmlns=\"http://www.opengis.net/kml/2.2\">\n")
  318 + f.write("<Document>\n")
  319 + for s in spots:
  320 + f.write("<Placemark><name>Spot " + str(s.esnName) + "</name><description>(" + str(s.latitude) + \
  321 + "," +str(s.longitude) \
  322 + + ")" + s.timeStampEST + "</description><Point><coordinates>" + str(s.longitude) + "," + \
  323 + str(s.latitude) + ",0</coordinates></Point></Placemark>")
  324 +
  325 + f.write("</Document></kml>")
  326 + f.close()
  327 +
  328 +############################################################
  329 +#
  330 +# Produire un fichier csv pour tous les spots
  331 +#
  332 +############################################################
  333 +
  334 +def ecrirecsv(fichier) :
  335 +
  336 + f = open(fichier, "ab")
  337 + writer = csv.writer(f,delimiter=';')
  338 +
  339 + # charge le shapefile du trait de cote du GSL
  340 + s = fiona.open("/home/nicopa01/data/spot/shp/GSL_GSHHS_f.shp")
  341 + shapefile_record = s.next()
  342 + shape = shapely.geometry.asShape(shapefile_record['geometry'])
  343 + x_min, y_min, x_max, y_max = shape.bounds
  344 +
  345 + for k in spots:
  346 + point = shapely.geometry.Point(k.longitude,k.latitude)
  347 +
  348 + if shape.contains(point):
  349 + k.flag = "TER"
  350 + else:
  351 + k.flag = "EAU"
  352 +
  353 + if (k.longitude < x_min) or (k.longitude > x_max) \
  354 + or (k.latitude < y_min) or (k.latitude > y_max):
  355 + k.flag = "OUT"
  356 +
  357 + data = [ str(k.esnName), str(k.latitude), str(k.longitude), int(k.timeGMT), str(k.messagetype), str(k.flag)]
  358 +
  359 + writer.writerow(data)
  360 +
  361 + f.close
  362 +
  363 +
  364 +
  365 +#############################################################
  366 +#
  367 +# Debut du programme principal
  368 +#
  369 +#############################################################
  370 +
  371 +if __name__ == "__main__":
  372 +# url = ' http://share.findmespot.com/messageService/guestlinkservlet?glId=03kp21oFnOj8bRHaPQPPAjf2TuO19VRa3'
  373 + url = ' https://api.findmespot.com/spot-main-web/consumer/rest-api/2.0/public/feed/0zgtS2Yrxgzl8in5LvVoAg0ae76LfQcdO/message.xml'
  374 +# positions = 'positions_mai_2013.xml'
  375 + getTags(url)
  376 +
  377 +#
  378 +#On boucle sur tous les spots pour calculer la distance
  379 +#entre eux.
  380 +#On forme une chaine de pres en pres
  381 +#
  382 +if len(spots) == 0:
  383 + sys.exit(1)
  384 +
  385 +a = spots[0]
  386 +while not a == None:
  387 + dist = 1000000.0
  388 + s = None
  389 + for b in spots:
  390 + distt = a.calcDistance(b,True)
  391 + if distt > 0 and distt < dist :
  392 + dist = distt
  393 + s = b
  394 + if not s == None:
  395 + a.distancePlusPres = dist
  396 + a.identPlusPres = s.esnName
  397 + a.capPlusPres = a.calcCap(s)
  398 + a.spotPlusPres = s
  399 + a.dejaTraite = True
  400 + a = s
  401 +
  402 +#ecrireHtml("./index.html")
  403 +#ecrireKml("./spots.kml")
  404 +ecrirecsv("/home/nicopa01/data/spot/spots_pos.csv")
  405 +#ecrirecsv("./spots.csv")
  406 +
  407 +sys.exit(0)
  408 +
... ...
tr_sort.m 0 → 100644
... ... @@ -0,0 +1,98 @@
  1 +function tr_sort(num_spot,mission,debut,h1,fin,h2)
  2 +% Compilation de données des bouées dérivantes par mission et date
  3 +%
  4 +% Attention, heure UTC!! Enregistrements des positions à partir du 20
  5 +% février 2015, certaines positions de spots sont manquantes... (début
  6 +% OHA)
  7 +%
  8 +% -» num_spot = numéro du/des spots à traiter (ex: [1] ou [45 46])
  9 +% -» mission = diminutif de la mission (3 caractères en capitale)
  10 +% -» debut = année, mois, jour de début de dérive (ex: 20150315)
  11 +% -» h1 = heure,min de début de dérive (hhmm)
  12 +% -» fin = année, mois, jour de fin de dérive
  13 +% -» h2 = heure,min de fin de dérive (hhmm)
  14 +% (facultatif, si non précisé -» données compilées jusqu'à maintenant)
  15 +%
  16 +% exemple : sort_spot([23 24],'BDC',20150612,1255,20150623,1835)
  17 +%
  18 +% Le fichier brute de données est situé sur Brandypot:
  19 +% /sas/usagers/share_lasso/data/MEOPAR/spots/spots_pos_clean.csv
  20 +%
  21 +% Enregistrements sous forme de structure sur Brandypot:
  22 +% sas/usagers/share_lasso/data/MEOPAR/drifter
  23 +
  24 +rootdir = '/sas/usagers/share_lasso/data/MEOPAR';
  25 +
  26 +fid1 = fopen([ rootdir 'spots_pos_clean.csv']);
  27 +spot_clean = textscan(fid1,'%f %f %f %f %s', 'delimiter',';');
  28 +fclose(fid1);
  29 +
  30 +spot = spot_clean{1,1};
  31 +lat = spot_clean{1,2}; lat(find(lat == -99999)) = NaN;
  32 +lon = spot_clean{1,3}; lon(find(lon == -99999)) = NaN;
  33 +nU = spot_clean{1,4}; % temps UNIX
  34 +nM = nU/86400 + datenum(1970,1,1); % temps matlab
  35 +
  36 +% date de début de dérive (en heure UTC)
  37 +dat1 = num2str(debut); H1 = num2str(h1,'%04i');
  38 +y1 = str2double(dat1(1:4)); m1 = str2double(dat1(5:6)); d1 = str2double(dat1(7:8));
  39 +h1 = str2double(H1(1:2)); mm1 = str2double(H1(3:4));
  40 +yr = dat1(3:4);
  41 +
  42 +dat1 = datenum(y1,m1,d1,h1,mm1,00);
  43 +
  44 +% date de fin de dérive (en heure UCT)
  45 +% si non spécifié, prise de données jusqu'à maintenant
  46 +if(~exist('fin','var'))
  47 + % reference à l'heure UTC
  48 + t = utc; dat2 = datenum(t(1),t(2),t(3),t(4),t(5),t(6));
  49 +else
  50 + dat2 = num2str(fin); H2 = num2str(h2,'%04i');
  51 + y2 = str2double(dat2(1:4)); m2 = str2double(dat2(5:6)); d2 = str2double(dat2(7:8));
  52 + h2 = str2double(H2(1:2)); mm2 = str2double(H2(3:4));
  53 +
  54 + dat2 = datenum(y2,m2,d2,h2,mm2,00);
  55 +end
  56 +
  57 +for i = 1:length(num_spot);
  58 +
  59 + filename = ['s' num2str(num_spot(i),'%03i') '_' mission yr ];
  60 +
  61 + I = find(spot(:) == num_spot(i));
  62 + tmp = [nM(I) lon(I) lat(I) nU(I)];
  63 +
  64 + K = find(tmp(:,1) > dat1 & tmp(:,1) < dat2);
  65 +
  66 + data = struct('timeM',{0},'lon',{''},'lat',{''},'timeU',{0},...
  67 + 'start',{''},'end',{''},'mission',{''},...
  68 + 'name',{''},'spot',{''});
  69 +
  70 + data.timeM = tmp(K,1); data.lon = tmp(K,2); data.lat = tmp(K,3);
  71 + data.timeU = tmp(K,4);
  72 +
  73 + data.start = [datestr(data.timeM(1)) ' UTC'];
  74 + if(~exist('fin','var'))
  75 + data.end = [''];
  76 + else
  77 + data.end = [datestr(data.timeM(end)) ' UTC'];
  78 + end
  79 +
  80 + data.spot = [ num_spot(i) ];
  81 + data.mission = [ mission ];
  82 + data.name = [ filename ];
  83 +
  84 + if(~exist([ rootdir '/drifter/' mission '/csv' ]));
  85 + mkdir([ rootdir '/drifter/' mission '/csv' ]);
  86 + mkdir([ rootdir '/drifter/' mission '/mat' ]);
  87 + end
  88 +
  89 + pathname = [ rootdir '/drifter/' mission '/mat/' filename '.mat'];
  90 + save(pathname,'data');
  91 +
  92 + csvname = [ rootdir '/drifter/' mission '/csv/' filename '.csv'];
  93 + CSV = struct2cell(data);
  94 + CSV = [ CSV(1,1) CSV(2,1) CSV(3,1) CSV(4,1)];
  95 + CSV = cell2mat(CSV);
  96 + dlmwrite(csvname,CSV,'precision', 11);
  97 +
  98 +end
... ...
utc.m 0 → 100644
... ... @@ -0,0 +1,67 @@
  1 +function [gmt,clk,dt] = utc;
  2 +
  3 +% UTC returns actual UTC-Time on a UNIX-System
  4 +%
  5 +% UTC_Time = UTC
  6 +%
  7 +% UTC_Time = [ YYYY MM DD hh mm ss ]
  8 +%
  9 +% use UNIX: date --utc +"%Y %m %d %H %M %S"
  10 +%
  11 +% A second Output returns the Vector of Matlab's CLOCK
  12 +%
  13 +% A third Output returns the Deviation: Local - UTC
  14 +%
  15 +% see also: CLOCK
  16 +%
  17 +
  18 +w = [];
  19 +
  20 +if isunix
  21 + [s,w] = unix('date --utc +"%Y %m %d %H %M %S"');
  22 + clk = clock;
  23 +else
  24 + warning('UNIX required.')
  25 + clk = clock;
  26 +end
  27 +
  28 +
  29 +ok = ~isempty(w);
  30 +if ok
  31 + w = eval(['[' w ']'],'NaN');
  32 + ok = ( isnumeric(w) & isequal(size(w),[1 6]) );
  33 +end
  34 +
  35 +if ok
  36 + gmt = w;
  37 +else
  38 + gmt = NaN * ones(1,6);
  39 +end
  40 +
  41 +
  42 +clk = floor(clk);
  43 +
  44 +
  45 +if ( nargout == 0 ) & ~isnan(gmt(1))
  46 +
  47 + fprintf(1,'%s\n',datestr(gmt,0))
  48 +
  49 + clear gmt
  50 +
  51 +elseif nargout == 3
  52 +
  53 + dt = cat( 1 , clk([3 4 5 6]) , gmt([3 4 5 6]) );
  54 +
  55 + dt = dt(1,:) - dt(2,:); % CLK - GMT
  56 +
  57 + if abs(dt(1)) > 1
  58 + dt(1) = -sign(dt(1));
  59 + end
  60 +
  61 + dt(4) = dt(4) .* ( abs(dt(4)) > 1 ); % Seconds
  62 +
  63 + dt = sum( dt .* [ 1 1/24 1/24/60 1/24/3600 ] );
  64 +
  65 +end
  66 +
  67 +
... ...