Commit 509db6045c4274912fafa946c4f6d61f4d267c80

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

maj tr_position.py + ajout de tr_snake.py (creation de cartes des suivi dans l'estuaire)

Showing 2 changed files with 190 additions and 191 deletions   Show diff stats
tr_position.py 100644 → 100755
... ... @@ -10,14 +10,12 @@ import fiona
10 10 import shapely.geometry
11 11 from time import localtime, strftime, ctime
12 12  
13   -
14   -
15 13 #
16 14 #Liste globale contenant tous les objets
17 15 #instancies a partir de la classe Spot
18 16 #
19   -spots = []
20 17  
  18 +spots = []
21 19  
22 20 class Spot:
23 21 '''Classe pour les GPS SPOT'''
... ... @@ -28,6 +26,7 @@ class Spot:
28 26 capPlusPres = 0.0
29 27 identPlusPres = 0
30 28 parentPlusPres = 0
  29 + ecartDistance = 0
31 30  
32 31 spotPlusPres = None
33 32 dejaTraite = False
... ... @@ -62,85 +61,6 @@ class Spot:
62 61  
63 62  
64 63  
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 64 #############################################################
145 65 #
146 66 #Lecture du fichier XML contenant les information des spots et
... ... @@ -253,81 +173,7 @@ def getTags(url):
253 173  
254 174 ############################################################
255 175 #
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
  176 +# implemente les positions des spots dans le csv
331 177 #
332 178 ############################################################
333 179  
... ... @@ -346,21 +192,30 @@ def ecrirecsv(fichier) :
346 192 point = shapely.geometry.Point(k.longitude,k.latitude)
347 193  
348 194 if shape.contains(point):
349   - k.flag = "TER"
  195 + k.flag = "LND"
350 196 else:
351   - k.flag = "EAU"
352   -
  197 + k.flag = "OCN"
  198 +
353 199 if (k.longitude < x_min) or (k.longitude > x_max) \
354 200 or (k.latitude < y_min) or (k.latitude > y_max):
355 201 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 202  
  203 + data = [ str(k.esnName), float(k.latitude), float(k.longitude), int(k.timeGMT), str(k.messagetype), str(k.flag)]
  204 +
359 205 writer.writerow(data)
360 206  
361 207 f.close
362 208  
363   -
  209 +def cleancsv(fichier_clean) :
  210 + # nettoyage des doublons
  211 + with open('/home/nicopa01/data/spot/spots_pos.csv','r') as in_file,\
  212 + open(fichier_clean,'w') as out_file:
  213 + seen = set()
  214 + for line in in_file:
  215 + if line in seen: continue # skip duplicate
  216 + seen.add(line)
  217 + out_file.write(line)
  218 +
364 219  
365 220 #############################################################
366 221 #
... ... @@ -374,35 +229,9 @@ if __name__ == &quot;__main__&quot;:
374 229 # positions = 'positions_mai_2013.xml'
375 230 getTags(url)
376 231  
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 232 ecrirecsv("/home/nicopa01/data/spot/spots_pos.csv")
405   -#ecrirecsv("./spots.csv")
  233 +cleancsv("/home/nicopa01/data/spot/spots_pos_clean.csv")
  234 +
406 235  
407 236 sys.exit(0)
408 237  
... ...
tr_snake.py 0 → 100755
... ... @@ -0,0 +1,170 @@
  1 +#!/usr/bin/env python
  2 +#! python
  3 +
  4 +from mpl_toolkits.basemap import Basemap
  5 +from time import gmtime, strftime
  6 +import matplotlib.pyplot as plt
  7 +import time
  8 +import numpy as np
  9 +import pandas as pd
  10 +from mpl_toolkits.axes_grid1 import make_axes_locatable
  11 +
  12 +strtime = strftime("%Y-%m-%d %H:%M", gmtime())
  13 +nametime = strftime("%Y%m%d%H%M", gmtime())
  14 +# temps unix
  15 +timestamp = int(time.time())
  16 +# sec dans une journee
  17 +sday = 60 * 60 * 24
  18 +# intervalle de temps
  19 +dd = 5 # nombre de jours
  20 +dt = timestamp - 5*sday # intervalle de temps depuis dd jours
  21 +diff = timestamp - dt #
  22 +
  23 +# colormap
  24 +cmap=plt.get_cmap('OrRd')
  25 +
  26 +
  27 +# Base de donnees ----------------------------
  28 +df = pd.read_csv('/home/nicopa01/data/spot/spots_pos.csv',delimiter=';')
  29 +df = df.drop_duplicates(take_last=True)
  30 +
  31 +
  32 +# Baie des Chaleurs ---------------------------
  33 +# traces des x derniers jours dans la baie
  34 +BDC = df[(df.SPOT > 699) & (df.SPOT < 800) & (df.FLAG == "OCN")
  35 + & (df.TIME > dt)]
  36 +
  37 +plt.clf()
  38 +fig, ax = plt.subplots()
  39 +
  40 +# caracteristiques de la carte
  41 +m = Basemap(projection='stere',lon_0=-68.5,lat_0=48,\
  42 + llcrnrlat=47.5,urcrnrlat=48.7,\
  43 + llcrnrlon=-66.9,urcrnrlon=-64.2,\
  44 + rsphere=6371.2,resolution='h')
  45 +
  46 +m.drawmapboundary(fill_color='lightblue') # fill to edge
  47 +m.drawcountries()
  48 +m.fillcontinents(color='grey',lake_color='lightblue',zorder=0)
  49 +m.drawrivers()
  50 +m.drawcoastlines(linewidth=0.8)
  51 +m.drawparallels(np.arange(40,60,0.5),labels=[1,0,0,0])
  52 +m.drawmeridians(np.arange(-80,-40,1),labels=[0,0,0,1])
  53 +
  54 +# trace les chemins de bouees
  55 +for spot, group in BDC.groupby(['SPOT']):
  56 + latitude = group.LAT.values
  57 + longitude = group.LON.values
  58 + t = (group.TIME.values - dt) / diff
  59 + x,y = m(longitude, latitude)
  60 + cax = plt.scatter(x,y,c=t, s=45, cmap=cmap,linewidth=0,vmin=0, vmax=1)
  61 + plt.plot(x,y,'k-',linewidth=0.5)
  62 +
  63 +# colorbar
  64 +cbar = fig.colorbar(cax,fraction=0.03, pad=0.03)
  65 +cbar.set_ticks([0,0.2,0.4,0.6,0.8,1])
  66 +cbar.update_ticks()
  67 +cbar.ax.set_yticklabels(['5','4','3','2','1','0'])
  68 +cbar.set_label('jours', rotation=270)
  69 +
  70 +# legende
  71 +plt.annotate(strtime + ' UTC', xy=(1, 1), xycoords='data', size=20,
  72 + xytext=(0.55, 0.1), textcoords='axes fraction',
  73 + horizontalalignment='right', verticalalignment='top')
  74 +plt.title(str(dd) + " derniers jours de derive de surface (BDC)")
  75 +# enregistrement
  76 +plt.savefig('/home/nicopa01/data/spot/suivi/BDC/BDC' + str(dd) + '_' + nametime + '.png',dpi=300)
  77 +
  78 +
  79 +# Estuaire du Saint-Laurent ---------------------------
  80 +# traces des x derniers jours dans l'estuaire
  81 +SLE = df[(df.SPOT > 1) & (df.SPOT < 300) & (df.FLAG == "OCN")
  82 + & (df.TIME > dt) | (df.SPOT > 899) & (df.FLAG == "OCN")
  83 + & (df.TIME > dt)]
  84 +
  85 +plt.clf()
  86 +fig, ax = plt.subplots()
  87 +
  88 +# caracteristiques de la carte
  89 +m = Basemap(projection='stere',lon_0=-68.5,lat_0=48,\
  90 + llcrnrlat=47.8,urcrnrlat=49.5,\
  91 + llcrnrlon=-70,urcrnrlon=-67,\
  92 + rsphere=6371.2,resolution='h')
  93 +
  94 +m.drawmapboundary(fill_color='lightblue') # fill to edge
  95 +m.drawcountries()
  96 +m.fillcontinents(color='grey',lake_color='lightblue',zorder=0)
  97 +m.drawrivers()
  98 +m.drawcoastlines(linewidth=0.8)
  99 +m.drawparallels(np.arange(40,60,0.5),labels=[1,0,0,0])
  100 +m.drawmeridians(np.arange(-80,-40,1),labels=[0,0,0,1])
  101 +
  102 +# trace les chemins de bouees
  103 +for spot, group in SLE.groupby(['SPOT']):
  104 + latitude = group.LAT.values
  105 + longitude = group.LON.values
  106 + t = (group.TIME.values - dt) / diff
  107 + x,y = m(longitude, latitude)
  108 + cax = plt.scatter(x,y,c=t, s=45, cmap=cmap,linewidth=0,vmin=0, vmax=1)
  109 + plt.plot(x,y,'k-',linewidth=0.5)
  110 +
  111 +# colorbar
  112 +cbar = fig.colorbar(cax,fraction=0.03, pad=0.03)
  113 +cbar.set_ticks([0,0.2,0.4,0.6,0.8,1])
  114 +cbar.update_ticks()
  115 +cbar.ax.set_yticklabels(['5','4','3','2','1','0'])
  116 +cbar.set_label('jours', rotation=270)
  117 +
  118 +# legende
  119 +plt.annotate(strtime + ' UTC', xy=(1, 1), xycoords='data', size=20,
  120 + xytext=(0.9, 0.1), textcoords='axes fraction',
  121 + horizontalalignment='right', verticalalignment='top')
  122 +plt.title(str(dd) + " derniers jours de derive de surface (SLE)")
  123 +#enregistrement
  124 +plt.savefig('/home/nicopa01/data/spot/suivi/SLE/SLE' + str(dd) + '_' + nametime + '.png',dpi=300)
  125 +
  126 +
  127 +# Old Harry ---------------------------
  128 +# traces des x derniers jours autour de Old Harry
  129 +OHA = df[(df.FLAG == "OCN") & (df.TIME > dt)]
  130 +
  131 +plt.clf()
  132 +fig, ax = plt.subplots()
  133 +
  134 +# caracteristiques de la carte
  135 +m = Basemap(projection='stere',lon_0=-60.5,lat_0=49,\
  136 + llcrnrlat=47.1,urcrnrlat=51.2,\
  137 + llcrnrlon=-62.9,urcrnrlon=-56.2,\
  138 + rsphere=6371.2,resolution='h')
  139 +
  140 +m.drawmapboundary(fill_color='lightblue') # fill to edge
  141 +m.drawcountries()
  142 +m.fillcontinents(color='grey',lake_color='lightblue',zorder=0)
  143 +m.drawrivers()
  144 +m.drawcoastlines(linewidth=0.8)
  145 +m.drawparallels(np.arange(40,60,1),labels=[1,0,0,0])
  146 +m.drawmeridians(np.arange(-80,-40,2),labels=[0,0,0,1])
  147 +
  148 +# trace les chemins de bouees
  149 +for spot, group in BDC.groupby(['SPOT']):
  150 + latitude = group.LAT.values
  151 + longitude = group.LON.values
  152 + t = (group.TIME.values - dt) / diff
  153 + x,y = m(longitude, latitude)
  154 + cax = plt.scatter(x,y,c=t, s=45, cmap=cmap,linewidth=0,vmin=0, vmax=1)
  155 + plt.plot(x,y,'k-',linewidth=0.5)
  156 +
  157 +# colorbar
  158 +cbar = fig.colorbar(cax,fraction=0.03, pad=0.03)
  159 +cbar.set_ticks([0,0.2,0.4,0.6,0.8,1])
  160 +cbar.update_ticks()
  161 +cbar.ax.set_yticklabels(['5','4','3','2','1','0'])
  162 +cbar.set_label('jours', rotation=270)
  163 +
  164 +# legende
  165 +plt.annotate(strtime + ' UTC', xy=(1, 1), xycoords='data', size=20,
  166 + xytext=(0.92, 0.1), textcoords='axes fraction',
  167 + horizontalalignment='right', verticalalignment='top')
  168 +plt.title(str(dd) + " derniers jours de derive de surface (BDC)")
  169 +# enregistrement
  170 +plt.savefig('/home/nicopa01/data/spot/suivi/OHA/OHA' + str(dd) + '_' + nametime + '.png',dpi=300)
... ...