dimanche 1 décembre 2019

BeSt - XSLT - mySQL

Il est possible de télécharger des listes de rues et de communes belges sur le site Open Data de BoSa (SPF Stratégie et Appui) : best-full-latest.zip (BeSt = Belgian Street).

Ce .zip contient des fichiers XML avec diverses informations concernant les rues et les communes de Belgique. C'est la référence officielle. Ces données sont aussi accessibles via un webservice (REST).

En ne prenant que la Wallonie, on a un peu plus que 61.000 rues et 263 communes. Ce qui va permettre de jouer un peu avec SQL.
J'ignore quelle est la meilleure manière d'injecter ces données XML dans une base de données mySQL. Alors, je le fais à ma manière... J'utilise un fichier de transformation XSL pour extraire les données qui m'intéressent et les mettre sous une forme facilement injectable dans une DB mySQL. Je crée une nouvelle base de données, best_street, Je crée deux tables (et comme je n'y connais rien, je leur donne de mauvais noms... : street_table et municipality_table.
$ mysql -u root -p
CREATE DATABASE best_street;
GRANT select, show view, lock tables, reload, replication client on *.* to 'bibi'@'localhost';
USE best_street;
CREATE TABLE street_table(
          street_id INT,
          street_name VARCHAR(100),
          municipality_id INT,
          PRIMARY KEY (street_id));
CREATE TABLE municipality_table(
          municipality_id INT,
          municipality_name VARCHAR(100),
          PRIMARY KEY (municipality_id));


Je récupère deux fichiers dans le .zip : WalloniaStreetname20191018.xml et WalloniaMunicipality20191018.xml. Ils contiennent des trucs du genre (attention, les fichiers des rues sont très grands, Firefox s'engorge) :

et


Que l'on peut traiter (sous Linux) avec
$ xsltproc x.xls y.xml > z.sql
Les fichiers .xsl étant, pour les rues,
<?xml version="1.0"?>
<xsl:stylesheet version="1.0"
 xmlns:tns="http://fsb.belgium.be/mappingservices/FullDownload/v1_00"
 xmlns:com="http://vocab.belgif.be/ns/inspire/"
 xmlns:xsl="http://www.w3.org/1999/XSL/Transform">
<xsl:output method="text"/>
<xsl:template match="/">
INSERT INTO street_table VALUES
      <xsl:for-each select="tns:StreetResponseBySource/tns:Streetname">
   (<xsl:value-of select="com:streetnameCode/com:objectIdentifier"/>,"<xsl:value-of select="com:streetname/com:spelling"/>",<xsl:value-of select="com:isAssignedBy/com:Municipality/com:objectIdentifier"/>),</xsl:for-each>
</xsl:template>

</xsl:stylesheet> 

Et pour les communes,
<?xml version="1.0"?>
<xsl:stylesheet version="1.0"
        xmlns:tns="http://fsb.belgium.be/mappingservices/FullDownload/v1_00"
        xmlns:com="http://vocab.belgif.be/ns/inspire/"
        xmlns:xsl="http://www.w3.org/1999/XSL/Transform">
<xsl:output method="text"/>
<xsl:template match="/">
INSERT INTO municipality_table VALUES
      <xsl:for-each select="tns:MunResponseBySource/tns:Municipality">
          (<xsl:value-of select="com:municipalityCode/com:objectIdentifier"/>,"<xsl:value-of select="com:municipalityName/com:spelling"/>"),</xsl:for-each>
</xsl:template>

</xsl:stylesheet>

Le résultat n'est pas parfait et quelques adaptations cosmétiques sont nécessaires pour en faire de bons fichiers SQL. Il y a aussi des 'streetname' qui contiennent des '"' et le .xsl est trop primitif pour en tenir compte. Et, en fait, ce sont des données incorrectes. Cela donne quelque chose comme :
INSERT INTO street_table VALUES
   (7700015,"Chemin des Prés",25005),
   (7700016,"Chemin des Ramiers",25005),
   (7700017,"Chemin des Soeurs",25005),
   (7700018,"Chemin du Grand Champ",25005),
   (7700019,"Chemin du Petit Brou",25005),
...
   (7741188,"Rue Jean-Louis Dumont",64074),
   (7741189,"Rue Joseph Nicolas",64074),
   (7741190,"Rue Joseph Noville",64074),
   (7741254,"Rue des Bâtis",64075),
   (7736185,"Tribomont",63035),
   (7736549,"Route de Manhay",63045);

et
INSERT INTO municipality_table VALUES
   (52074,"Aiseau-Presles"),
   (61003,"Amay"),
   (63001,"Amel"),
   (92003,"Andenne"),
   (56001,"Anderlues"),
   (91005,"Anhée"),
...
   (25110,"Waterloo"),
   (25112,"Wavre"),
   (63084,"Welkenraedt"),
   (84075,"Wellin"),
   (91141,"Yvoir");

Il faut aussi noter qu'ils n'ont pas été fichus de donner des street_id distincts couvrant les trois régions. Il y a des duplicates du côté des 36000 et 37000 (!).
Il suffit maintenant d'injecter ces données dans notre DB :
SOURCE insert_streets.sql;
SOURCE insert_municipality.sql;

Notre DB est prête, on peut commencer à jouer avec...
mysql> SELECT * FROM street_table LIMIT 5;
+-----------+------------------------+-----------------+
| street_id | street_name            | municipality_id |
+-----------+------------------------+-----------------+
|   7700001 | Aux Tiennes            |           25005 |
|   7700002 | Avenue des Bouleaux    |           25005 |
|   7700003 | Avenue des Cerisiers   |           25005 |
|   7700004 | Avenue des Combattants |           25005 |
|   7700005 | Avenue des Pruniers    |           25005 |
+-----------+------------------------+-----------------+
5 rows in set (0.00 sec)

mysql> SELECT * FROM municipality_table ORDER BY municipality_name DESC LIMIT 5;
+-----------------+-------------------+
| municipality_id | municipality_name |
+-----------------+-------------------+
|           91141 | Yvoir             |
|           84075 | Wellin            |
|           63084 | Welkenraedt       |
|           25112 | Wavre             |
|           25110 | Waterloo          |
+-----------------+-------------------+
5 rows in set (0.00 sec)

mysql> 
Par exemple, je peux avoir la liste de commune ayant le moins ou le plus de rues :
mysql> SELECT count(*), m.municipality_name from street_table s LEFT JOIN municipality_table m USING (municipality_id) GROUP BY municipality_id ORDER BY count(*) LIMIT 10;
+----------+-------------------+
| count(*) | municipality_name |
+----------+-------------------+
|       54 | Martelange        |
|       58 | Hélécine          |
|       61 | Stoumont          |
|       62 | Gouvy             |
|       63 | Lincent           |
|       64 | Flobecq           |
|       65 | Verlaine          |
|       65 | Rumes             |
|       66 | Donceel           |
|       66 | Oreye             |
+----------+-------------------+
10 rows in set (0.38 sec)

mysql> SELECT count(*), m.municipality_name from street_table s LEFT JOIN municipality_table m USING (municipality_id) GROUP BY municipality_id ORDER BY count(*) DESC LIMIT 10;
+----------+----------------------------+
| count(*) | municipality_name          |
+----------+----------------------------+
|     2117 | Charleroi                  |
|     1998 | Liège                      |
|     1730 | Namur                      |
|     1280 | Mons                       |
|     1086 | Tournai                    |
|      923 | La Louvière                |
|      812 | Ottignies-Louvain-la-Neuve |
|      738 | Wavre                      |
|      677 | Mouscron                   |
|      649 | Ath                        |
+----------+----------------------------+
10 rows in set (0.10 sec)

mysql> 

Cette dernière requête étant :
SELECT count(*), m.municipality_name from street_table s
            LEFT JOIN municipality_table m USING (municipality_id)
            GROUP BY municipality_id
            ORDER BY count(*)
            DESC LIMIT 10;



lundi 18 novembre 2019

pps

Je me demandais à quel point les signaux PPS de deux modules GPS à quelques euros étaient synchronisé...
Trigger sur le channel-1. Il y a du jitter sur l'autre trace.
La réception n'est pas idéale à l'endroit où j'ai effectué la mesure, sur une terrasse entre les maisons. Trigger sur le channel-2. Idem.
Un autre point à noter est que les modules semblent interférer entre eux. J'ai utilisé deux alimentations différentes et j'ai dû en démarrer un avant l'autre... Une question intéressante serait de savoir si le jitter diminue quand on a une meilleure réception avec plus de satellites visibles en même temps.

dimanche 15 septembre 2019

NanoVNA

Un 'Vector Network Analyzer' à 40.00 euros sur eBay... Ce petit bidule 'stand alone' permet de mesurer la réflexion (S11) et la transmission de signaux entre 50 kHz et 900 MHz.

Mesure en transmission d'un filtre 'SAW' passe-bande 433 MHz.

Le nanoVNA est équipé d'un petit 'touchscreen' mais on peut également communiquer via le port USB (USBc). Sur Linux, on peut utiliser 'minicom -D /dev/ttyACM0', par exemple. Et sur Windows, 'putty'.

L'interface est bien conçue. CH0 mesure la réflexion (S11) et CH1, la transmission (S21) en valeur imaginaires. On définit la plage de fréquences (sweep) qui sera divisée en 101 points. On peut afficher 4 traces. Pour chaque trace, on peut choisir le canal, le type d'affichage (logmag, smith, phase), etc...

Via le port USB, on a une console série dans laquelle on peut taper des commandes et afficher les résultats. Ainsi, par exemple, on peut facilement obtenir les mesures en transmission entre 400 et 500 MHz.
ch> sweep 400000000 500000000 101
ch> data 1
-0.004106074 0.001582601
-0.004401572 0.001613250
-0.004341298 0.001819048
-0.004164417 0.001013602
-0.006044433 0.001735574
...
Un petit script PERL vite écrit, convertit les valeurs en dB :
#!/usr/bin/perl
$f0 = 400.0; $f1 = 500.0;
$inc = ($f1 - $f0) / 101;
$df = 0;
while ($line = <STDIN>)
 {
 chomp $line;
 ($a, $b) = split /\s/,$line;
 $mag = sqrt ($a*$a + $b*$b);
 if ($mag == 0) {$mag = 0.001;}
 printf("%6.2f\t%f\n", $f0 + $df, 20*log($mag)/log(10.0));
 $df += $inc;
 }
et on obtient le même graphe que l'écran du nanoVNA avec gnuplot :
gnuplot <<END
set term png size 1200,800
set output 'S21-logmag.png'
set title 'S21-logmag of a "NMRF FBp - 433s" filter'
set xlabel 'Frequency (MHz)'
set ylabel 'S21 dB'
plot 'S21-logmag' with lines
END
Graphe calculé à partir des données obtenues via l'interface série.

C'est magique!

Quelques exemples pour la route :
ch>
ch> help
Commands: help exit info echo systime threads reset freq offset time dac saveconfig
 clearconfig data dump frequencies port stat sweep test touchcal touchtest pause
 resume cal save recall trace marker edelay 

ch> info
Kernel:       4.0.0
Compiler:     GCC 5.4.1 20160919
Architecture: ARMv6-M
Core Variant: Cortex-M0
Port Info:    Preemption through NMI
Platform:     STM32F072xB Entry Level Medium Density devices
Board:        NanoVNA
Build time:   Aug  2 2019 - 16:40:01

ch> trace
0 LOGMAG CH0 1.000000000 7.000000000
1 LOGMAG CH1 1.000000000 7.000000000
2 SMITH CH0 1.000000000 0.000000000
3 PHASE CH1 1.000000000 4.000000000
ch> trace 99
trace {0|1|2|3|all} [logmag|phase|smith|linear|delay|swr|off] [src]
ch> trace 1 phase 0                                                                
ch> trace                                                                          
0 LOGMAG CH0 1.000000000 7.000000000                                               
1 PHASE CH0 1.000000000 4.000000000                                                
2 SMITH CH0 1.000000000 0.000000000                                                
3 PHASE CH1 1.000000000 4.000000000
                                               
ch> sweep                                                                          
400000000 500000000 101                                                            
ch> sweep 420000000 440000000 101                                                  
ch> sweep                                                                          
420000000 440000000 101                                            
ch>

ch> data 1                                                                         
-0.004733881 0.003933506                                                           
-0.004605576 0.004513134                                                           
-0.004246967 0.003887655                                                           
-0.004040227 0.003405604                                                           
-0.003938771 0.004323835                                                           
-0.004104804 0.004129481
...
data0 pour les mesures de réflexion, data 1 pour la transmission. Après, ce sont les données de calibration.

Quelques liens

mardi 27 août 2019

XSLT - Garmin GPX - gnuplot

Un Garmin donne la trace GPS suivante : AiguillesRouges.gpx
<?xml version="1.0" encoding="UTF-8" standalone="no" ?>
<gpx xmlns="http://www.topografix.com/GPX/1/1" creator="MapSource 6.16.3" version="1.1"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd">
  
  <metadata>
    <link href="http://www.garmin.com">   
      <text>Garmin International</text>
    </link>
    <time>2019-07-28T08:08:05Z</time>
    <bounds maxlat="46.024869801476598" maxlon="6.926982309669256" minlat="45.935957757756114" minlon="6.771008670330048"/>
  </metadata>
  
  <trk>
    <name>Tracé actuel: 23 JUIL 2019 07:39</name>
    <extensions>
      <gpxx:TrackExtension xmlns:gpxx="http://www.garmin.com/xmlschemas/GpxExtensions/v3">
        <gpxx:DisplayColor>Black</gpxx:DisplayColor>
      </gpxx:TrackExtension>
    </extensions>
    <trkseg>
      <trkpt lat="45.936855375766754" lon="6.771008670330048">
        <ele>872.52999999999997</ele>
        <time>2019-07-23T05:39:31Z</time>
      </trkpt>
      <trkpt lat="45.936975907534361" lon="6.77103029564023">
        <ele>873.98000000000002</ele>
        <time>2019-07-23T05:39:42Z</time>
      </trkpt>
      ...etc...
gpx.xsl, un petit bout de code XSLT
<?xml version="1.0"?>
<xsl:stylesheet version="1.0"
        xmlns:xsl="http://www.w3.org/1999/XSL/Transform"
        xmlns:gpx="http://www.topografix.com/GPX/1/1">
  <xsl:output method="text"/>
  <xsl:template match="/">
    <xsl:for-each select="gpx:gpx/gpx:trk/gpx:trkseg/gpx:trkpt">
        <xsl:value-of select="@lon"/>
        <xsl:text> </xsl:text>
        <xsl:value-of select="@lat"/>
        <xsl:text> </xsl:text>
        <xsl:value-of select="gpx:ele"/>
        <xsl:text>
</xsl:text>
    </xsl:for-each>
  </xsl:template>
</xsl:stylesheet>
Et sur Linux,
$ xsltproc gpx.xsl AiguillesRouges.gpx > AigRou.xyz
donne les coordonnées X, Y et Z en 3 colonnes :
6.771008670330048 45.936855375766754 872.52999999999997
6.77103029564023 45.936975907534361 873.98000000000002
6.771165663376451 45.937128961086273 876.38
6.771242022514343 45.937253851443529 877.34000000000003
6.771263061091304 45.93731545843184 877.34000000000003
6.771297678351402 45.937489299103618 879.25999999999999
etc...
Qu'il suffit de traiter avec le script gnuplot suivant :
gnuplot << END
set term png size 1200,900
set grid xtics ytics ztics
set output 'AigRouXYZ.png'
splot "AigRou.xyz" with lines
set output 'AigRouXY.png'
plot "AigRou.xyz" using 1:2 with lines
set output 'AigRouZ.png'
plot "AigRou.xyz" using 3 with lines
END
pour obtenir les graphes suivants :

samedi 27 avril 2019

ADS-B avec cantenna

Ayant trouvé un bidon d'huile d'olive de 10 litres, j'envisage d'en faire une 'cantenna' pour l'ADS-B... À l'image des 'antennes Ricoré' pour le WiFi.
Le bidon fait environ 15x24 cm² avec 30 cm de profondeur. L'ADS-B, à 1090 MHz, a une longueur d'onde d'environ 27.5 cm. Le bidon n'est pas très loin d'avoir une hauteur de lambda/2 et en mettant un bout de fil de cuivre de lambda/4 (~7cm) à lambda/4 du fond, ça devrait être bon. Aussitôt dit, aussitôt fait.
Je pars faire des tests dans un endroit relativement élevé, un terril au dessus de Liège. Et j'en registre les avions avec dump1090 sur un OrangePiZero alimenté par un powerbank. Trois fois pendant une demi-heure dans des directions différentes en comparant avec une antenne 'PCB' trouvée sur eBay pendant une autre demi-heure.
Le résultat est assez bon. La cantenna capte les avions beaucoup plus loin que l'antenne PCB (mais dans une direction privilégiée). J'ignore si c'est l'amplification du signal par l'antenne ou si c'est parce que la réception est moins dérangée par d'autres émissions proches venant d'autres directions. Il faudrait mesurer le nombre de messages par seconde pour voir s'il y a beaucoup de collisions. On se rapproche de l'horizon radar théorique (D_miles = 2.2 * sqrt(altitude_m)).
De Liège, je capte des avions au delà de Paris et presque jusque Londres.

dimanche 17 février 2019

AIS again

Petite analyse d'une captation d'une minute avec un Orange Pi Zero et une clé RTL-SDR V3 autour de 162 MHz de 16000000 échantillons à 262144 échantillons par seconde (environ une minute donne un fichier de 32 MB).
/usr/bin/rtl_sdr -f 162000000 -s 262144 -n 16000000 core.iq
Je prends juste le quart du fichier pour accélérer le traitement dans GNU Octave (dd). Pour lire les données IQ 8 bits fournies par rtl_sdr, j'utilise une fonction dans loadFile.m :
function y = loadFile(filename)
%  y = loadFile(filename)
%
% reads  complex samples from the rtlsdr file 
%

fid = fopen(filename,'rb');
y = fread(fid,'uint8=>double');

y = y-127;
y = y(1:2:end) + i*y(2:2:end);
Avec plotspec.m,
% plotspec(x,Ts) plots the spectrum of the signal x
% Ts = time (in seconds) between adjacent samples in x
function plotspec(x,Ts)
N=length(x);                               % length of the signal x
t=Ts*(1:N);                                % define a time vector
ssf=(ceil(-N/2):ceil(N/2)-1)/(Ts*N);       % frequency vector
fx=fft(x(1:N));                            % do DFT/FFT
fxs=fftshift(fx);                          % shift it for plotting
subplot(2,1,1), plot(t,x)                  % plot the waveform
xlabel('seconds'); ylabel('amplitude')     % label the axes
subplot(2,1,2), plot(ssf,abs(fxs))         % plot magnitude spectrum
xlabel('frequency'); ylabel('magnitude')   % label the axes
je peux observer que j'ai bien enregistré des messages et qu'ils se situent bien à -25 KHz et + 25 KHz autour de 162 MHz.
Je vais maintenant tacher d'isoler un message à +25 KHz (up.iq8) et un autre à -25 KHz (down.iq8). Pour ce faire, je vais visualiser la magnitude du signal.
Zoomons...
Zoomons encore...
On peut maintenant extraire juste les messages du fichier core.iq8. Sachant qu'il y a deux bytes par échantillon et que le graphe nous donne les offsets absolus.
$ # 134000..141000  & 274000..281000
$ dd if=core.iq8 of=up.iq8 bs=2000 skip=134 count=7
$ dd if=core.iq8 of=down.iq8 bs=2000 skip=274 count=7
En cherchant un peu, on trouve un message à +25 Khz dont on s'intéresse à l'évolution de la phase :
En regardant la phase et en zommant :
Et la même chose pour un message à -25 KHz :
qui donne une évolution de la phase descendante :
On voit que, dans les deux cas, on a des sauts quand la phase passe au dessus de 360° ou en dessous de 0°. On va rendre tout cela continu pour, ensuite, soustraire/ajouter le changement de phase dû au +/- 25 KHz.
%  y = contzup(phi)
%
% rend la phase continue
%
function y = contzup(phi)
N=length(phi);                               % length of the signal x
phi=phi+(pi/2);
y=phi;
t=0;
for i = 1:N-1
        if ( phi(i) > phi(i+1) )
                        t = t + 2*pi;
        endif
        y(i+1) = (phi(i+1) + t)-(((2*pi)*25000*i)/262144);
%       y(i+1) = (phi(i+1) + t);
endfor


Détail, le début avec le préambule et les premiers octets.


avec le grid à 9600 samples/sec :



Comme je ne parviens pas vraiment à écrire un programme me permettant d'obtenir les trames NMEA à partir des IQ sortant de la clé RTL-SDR en bouquinant sur la théorie. J'envisage de partir d'un programme qui fonctionne et de l'analyser... Par exemple, rtl_ais donne d'assez bons résultats. D'un terril au dessus de Liège avec une antenne dipôle, je capte des bateaux sur la Meuse jusque Tihange (Huy) dans un sens et au delà de Maastricht dans l'autre. Et même quelques message provenant du côté de Genk, sur le canal Albert. Le code n'est pas très long et est assez lisible mais, c'est vraiment de la sorcellerie. Il va falloir du temps pour s'en imprégner et comprendre son fonctionnement.

À suivre...