Automatizovaný příjem snímků NOAA

Nejprve byla anténa QFH, pak softwarové rádio s čipsetem Realtek RTL2832U, pak byl software na ovládání rádia zvaný gqrx, který je založený na gnuradio a pak tu byl uživatel, který musel čekat, až přiletí satelit, čekat čtvrt hodiny a sehnat několik příkazů na vytvoření obrázku pomocí programu wxtoimg. Uživatele to přestalo bavit a tak si napsal skript, který sám sestaví obrázek a nahraje jej na web server. Pak ale musel napsat skript, který ty obrázky seřadí podle data pořízení. Jenže ani to mu moc nepomohlo a nezbylo nic jiného, než přemýšlet jak automatizovat úplně všechno. První problém bylo sehnat čas příletu družice na viditelnou část oblohy. Z GPredict se nějak viditelně vydolovat nedal. Díky dobré duši, která implementovala GPredict do jazyka PHP a zveřejnila kód se dala první část odškrtnout. Zdrojové kódy stáhnete zde. Predict.php všechno zařizuje, pro vaše potřeby využijte souboru visible_passes.php ve složce examples. Tento soubor si jen poupravte dle svých potřeb. Já potřeboval čas příletu a dobu, po kterou satelit zůstává na viditelné části oblohy. Tento skript neběží na webserveru, je volán cronem vždy po 4. hodině ranní, kdy začíná první část přeletů viditelných od nás a zkontroluje přelety na dalších 24 hodin. Zatím běží na klasickém PC, běh na Raspberry Pi ještě vyzkouším. Pokud skript neběží, možná jsem udělal chybu zde při úpravách, ale je téměř samovysvětlující, takže není problém jej opravit/předělat podle svých potřeb.

noaa.php

<?php
/**
 * This is an example of how to use Predict for determining upcoming visible
 * passes of the International Space Station.  Please read the inline comments
 * for details.
 *
 * Run this from the root checkout, not from examples, or the include paths
 * will not work.
 */

require_once 'Predict.php';
require_once 'Predict/Sat.php';
require_once 'Predict/QTH.php';
require_once 'Predict/Time.php';
require_once 'Predict/TLE.php';

function sortByAos($a, $b) {
   if ($a[2] == $b[2]) {
        return 0;
    }
    return ($a[2] < $b[2]) ? -1 : 1;
}
// Track execution time of this script
$start = microtime(true);

// The observer or groundstation is called QTH in ham radio terms
$predict  = new Predict();
$qth      = new Predict_QTH();
$qth->alt = 300; // Altitude in meters

// Munich, example east of the meridian
$qth->lat = 48.1505; // Lat North
$qth->lon = 11.5809; // Lon East

// The iss.tle file is the first 3 lines of
// http://celestrak.com/NORAD/elements/stations.txt
// Make sure you update this content, it goes out of date within a day or two
exec("wget http://celestrak.com/NORAD/elements/weather.txt -O ~USERNAME/.wxtoimg/weather.txt");
$tleFile = file('http://celestrak.com/NORAD/elements/noaa.txt'); // Load up the NOAA from celestrak
$indexes = array(19, 18, 15);
$list = array();

foreach ($indexes as $index){
	$index *= 3;  //3 rows for satellite
	$tle     = new Predict_TLE(substr($tleFile[$index], 5, 2), $tleFile[$index+1], $tleFile[$index+2]); // Instantiate it
	$sat     = new Predict_Sat($tle); // Load up the satellite data
	$now     = Predict_Time::get_current_daynum(); // get the current time as Julian Date (daynum)

	// You can modify some preferences in Predict(), the defaults are below
	//
	 $predict->minEle     = 0; // Minimum elevation for a pass
	 $predict->timeRes    = 5; // Pass details: time resolution in seconds
	 $predict->numEntries = 10; // Pass details: number of entries per pass
	 $predict->threshold  = 90; // Twilight threshold (sun must be at this lat or lower)

	// Get the passes and filter visible only, takes about 4 seconds for 10 days
	$results  = $predict->get_passes($sat, $qth, $now, 1);
	$filtered = $predict->filterVisiblePasses($results);

	$zone   = 'Europe/Berlin'; // Pacific time zone
	$format = 'H:i mdy';         // Time format from PHP's date() function
	$epoch  = 'U';

	// Format the output similar to the heavens-above.com website
	//echo "<tr><td>Sat\t<td>AOS\t<td>ELE\t<td>LOS\t</tr>";
	foreach ($filtered as $pass) {
	    if (round($pass->max_el) < 17){
	       continue;
	    }
	    $list[] = array($pass->satname, Predict_Time::daynum2readable($pass->aos, $zone, $format), Predict_Time::daynum2readable($pass->aos, $zone, $epoch), Predict_Time::daynum2readable($pass->los, $zone, $epoch), round($pass->max_el));
	}
}

//delete the pass with worse max elevation from overlaping ones
usort($list, 'sortByAos');
$i=0;
while ($i < (sizeOf($list)-1)){
	if ($list[$i][3] >= $list[$i+1][2]) {
		if ($list[$i][4] > $list[$i+1][4]){
			array_splice($list, $i+1, 1);
		} else {
			array_splice($list, $i, 1);
		}		
		continue;
	}
	$i++;
}

//echo $wake_date > /sys/class/rtc/rtc0/wakealarm

foreach ($list as $pass) {
    $cmd = "echo \"cd /home/USERNAME && ./apt ".$pass[0]. " ";
    $cmd .= $pass[3] - $pass[2] - 60;
    $cmd .= "\" | at -m ";
    $cmd .= "\"". $pass[1] . " +1 minute\"\n";
    
    exec($cmd);
    //echo $cmd;
}
//echo "Execution time:  " . number_format((microtime(true) - $start) * 1000, 2) . "ms\n"; 
exit;

Zde vidíte, že php skript stáhnete TLE pro sebe a wxtoimg, Predict vypočítá potřebné časy, které php zapíše do příkazu at. At spustí mnou připravený skript, který zaznamená vysílání ze satelitu a přechroupe záznam do podoby obrázku. Ten nahraje na webovou prezentaci meteostanice pomocí nástroje curl. Nezapomeňte upravit username.
Bash skript původně zpracovával nahrávky z gqrx, z této doby pochází ta skvělá konstrukce na čtení datového razítka na začátku skriptu. Využívá programu rtl_fm z balíčku od osmocom. Dále tam najdete sox (třeba nainstalovat) pro práci se zvukovým záznamem a samozřejmě WXtotImg (také doinstalovat).

apt.sh

#!/bin/bash

bandwidth=150k
login='LOGIN:PASSWORD'

#Freq +- 0.006 ?
case $1 in
15)freq=137.54M;;
17)freq=137.5M;;
18)freq=137.9065M;;
19)freq=137.092M;;
*)exit 1;;
esac

filename=`date -u +%Y%m%d-%H%M%S`.s16
timeout $2 /usr/local/bin/rtl_fm -f $freq -s $bandwidth > $filename
#timeout $2 /usr/local/bin/rtl_sdr $filename.bin -f $freq -s 1.8e6
if [ $? -ne 124 ]; then
	rm $filename
	exit 1
fi

if [ ! -f $filename ]; then
	echo "Soubor $filename neexistuje."
	exit 2
fi

#set variables 
date=`echo $filename | cut -d- -f1`
time=`echo $filename | cut -d- -f2 | cut -d. -f1`
year=`echo $date | cut -c1-4`
month=`echo $date | cut -c5-6`
day=`echo $date | cut -c7-8`
hour=`echo $time | cut -c1-2`
minute=`echo $time | cut -c3-4`
second=`echo $time | cut -c5-6`

output=$year$month$day$hour$minute$second.wav
hvc=noaa$1-$year$month$day-$hour$minute$second-hvc.png
no=noaa$1-$year$month$day-$hour$minute$second-no.png
mcir=noaa$1-$year$month$day-$hour$minute$second-mcir.png
tmp=noaatmp.wav
map=map.png
thumbnail=noaa.png

#QTH
lat=48.1505
lon=11.5809
ele=300
#dawn=7
#dusk=17
#day/night is determined by python script

sox -r $bandwidth -c 1 -L $filename $tmp
sox $tmp $output rate 11025

touch -r $filename $output

/usr/local/bin/wxmap -T "NOAA $1" -L"$lat/$lon/$ele" -G /home/USERNAME/.wxtoimg -p 0 -l 0 -od "$day $month $year $hour:$minute:$second" $map

/usr/local/bin/wxtoimg -c -e MCIR -m $map -t "NOAA $1" -o $output $mcir
/usr/bin/curl -T $mcir ftp://ftp.server/ --user $login

if [ $? -eq 0 ]; then
	convert $mcir -resize 25% $thumbnail
	/usr/bin/curl -T $thumbnail ftp://ftp.server/ --user $login
	
	/usr/local/bin/wxtoimg -c -e NO -m $map -A -t "NOAA $1" -o $output $no
	/usr/bin/curl -T $no ftp://ftp.server/ --user $login
	
	#hourt=`echo $hour | sed 's/^0//'`
	python ./sun.py
	if [ $? -eq 1 ]; then
		/usr/local/bin/wxtoimg -c -e HVC -m $map -A -K -t "NOAA $1" -o $output $hvc
		/usr/bin/curl -T $hvc ftp://ftp.server/ --user $login
		rm $hvc
	fi	
	rm $no $mcir $thumbnail $output
fi

rm $tmp $filename $map

exit 0

Pro určení zda je den či noc používám ne zcela přesný výpočet místního hvězdného času a jeho porovnání vypočteného času východu a západu Slunce, opět hvězdného. To vše kvůli senzoru 4, který se odesílá až hodinu po východu Slunce, zde je ještě prostor pro vylepšení.
sun.py

#!/usr/bin/python

from datetime import datetime
import time
import math

def days_between(d1, d2):
    return (d2 - d1).days + (d2 - d1).seconds/84600.0
    
def echo_time(t):
    print (math.modf(t)[1])+math.ceil((math.modf(t)[0])*60)/100

#QTH
lat = math.radians(48.1505)
lon = math.radians(11.5809)
ele = 300
utcoffset = 1

#Twilight
now = datetime.utcnow()
equinox = datetime.strptime("03-20", "%m-%d")
equinox = equinox.replace(now.year)

if days_between(equinox, now) < 0:
	equinox = equinox.replace(now.year - 1)

d = days_between(equinox, now)
lamb = math.radians(360.0/365.2422 * d)
print "Lambda=", math.degrees(lamb), "d=", d

beta = 0.0
epsilon = math.radians(23.5)

delta = math.asin(math.sin(lamb)*math.sin(epsilon))
alfa = math.acos(math.cos(lamb)/math.sqrt(1-math.pow(math.sin(lamb)*math.sin(epsilon),2)))
if math.sin(lamb)*math.cos(epsilon)/math.sqrt(1-math.pow(math.sin(lamb)*math.sin(epsilon),2)) < 0:
	alfa = 2*math.pi - alfa

ts = math.acos(-math.tan(lat)*math.tan(delta))
if ts>math.pi:
	t1 = ts
	t2 = 2.0*math.pi - ts	
else:
	t1 = 2.0*math.pi - ts
	t2 = ts
	
t1 = math.degrees(t1)
t2 = math.degrees(t2)
while t1 > 360:
    t1 -= 360
while t2 > 360:
    t2 -= 360
print alfa
riselst = t1+math.degrees(alfa)
dawnlst = t2+math.degrees(alfa)
while riselst > 360:
    riselst -= 360
while dawnlst > 360:
    dawnlst -= 360

print "rise in angle",
echo_time(t1/15.0)
print "in LST",
echo_time(riselst/15.0)
print "set in angle",
echo_time(t2/15.0)
print "in LST",
echo_time(dawnlst/15.0)
	
#JD
r=now.year
m=now.month
d=now.day + (now.hour+(now.minute+now.second/60)/60)/24
if m > 2:
    r = r - 1
    m = m + 12
JD = 1720994.5 + 365.25*r + r/400.0 +- r/100.0 + 30.6*(m+1) + d + 2.0
print "JD=",JD

#JD best
unix = time.time()
J = 2440587.5 + unix/86400.0
T = (J - 2451545.0) / 36525.0
TT = 64.184 + 59.0 * T - 51.2 * math.pow(T,2) - 67.1 * math.pow(T,3) - 16.4 * math.pow(T,4)
JTT = J + TT / 86400.0
d = JTT - 2451545.0
print "JD=",JTT
print d

# Sidereal time
J2k = datetime.strptime("2000-01-01 12:00 UTC", "%Y-%m-%d %H:%M %Z")
#d = days_between(J2k, now)
print d

GMST = 280.46061837 + 360.98564736629 * d
while GMST > 360:
    GMST -= 360
LMST = 280.46061837 + 360.98564736629 * d + math.degrees(lon)
while LMST > 360:
    LMST -= 360
print "GMST",
echo_time(GMST/15.0)
print "LMST",
echo_time(LMST/15.0)

if LMST > riselst or LMST < dawnlst:
    print "Day"
    exit(1)
else:
    print "Night"
    exit(2)

t = GMST + math.degrees(lon) + 1.00273790935*(now.hour+now.minute/60 - 1)
print "t=",
echo_time(t/15.0)

h = math.degrees(math.asin(math.sin(delta)*math.sin(lat)+math.cos(delta)*math.cos(math.radians(LMST))*math.cos(lat)))
print "h=", h, "d=",
echo_time(math.degrees(delta))

2 komentáře u „Automatizovaný příjem snímků NOAA

  1. David

    Mám identickou sestavu. Tedy RTL USB s QFH anténou. Daří se mi chytat NOAA družice přes gqrx. Nicméně pokud se pokouším příjmout družici přes utilitu rtl_fm, nebyl jsem nikdy úspěšný. Ve výsledném wav souboru po všech úpravách pomocí sox byl sice někaký náznak zvuku, ale zdaleka nic co by šlo dekódovat.

    Skript jsem nepoužíval a zatím jsem to jen spouštěl ručně:
    rtl_fm -f 137.5M -s 150k > noaa.s16
    sox -r 150k -c 1 -L noaa.s16 noaatmp.wav
    sox noaatmp.wav noaa.wav rate 11025

    Netušíte kde by mohla být chyba?

    1. Honza Autor příspěvku

      Už to dlouho aktivně neprovozuji, takže si jen matně vzpomínám, že z rtl_fm byl horší výstup než z gqrx, které šlo nastavovat ručně během přeletu, jde o to, že se lehce mění frekvence a z rtl_fm pak není tak kvalitní výstup, ale zase je to automatické. Také si vzpomínám, že velkou roli hrálo zesílení, ale protože naše příkazy vypadají stejně, tak si nejsem jistý, kde je chyba.

Napsat komentář: Honza Zrušit odpověď na komentář

Vaše e-mailová adresa nebude zveřejněna. Vyžadované informace jsou označeny *