Wykład “Badania sejsmiczne Polski w 3D”

Jeszcze przed obroną doktoratu zostałem poproszony o wygłoszenie wykładu o tematyce mojej pracy doktorskiej w ramach studenckiego koła geofizyki na wydziale fizyki uniwersytetu warszawskiego.

Było to dla mnie ważne z dwóch powodów. Po pierwsze był to pierwszy wykład w ramach koła w tym roku akademickim, a po drugie było to moje pierwsze wystąpienie jako doktora. OK, nie czepiajmy się szczegółów, doktorem zostanę po uchwale rady wydziału, ale to tylko kwestia formalna.

Wszystkich zainteresowanych zapraszam do oglądania. Niestety kamera sportowa zawiodła mnie po raz kolejny – tym razem jasność nie pozwala na odczytanie niczego ze slajdów, więc znów dodałem je w montażu. Druga kamera również mnie zawiodła przerywając nagrywanie po 30 minutach. Mam nadzieję, że kiedyś nauczę się jak porządnie nagrywać.

| Komentarze

Jak w poniedziałek zostałem doktorem

W poniedziałek 17 października broniłem na Wydziale Fizyki Uniwersytetu Warszawskiego swoją pracę doktorską zatytułowaną „Lokalizacja ognisk trzęsień ziemi metodą propagacji wstecznej z wykorzystaniem implementacji metody fast marching w modelu 3D obszaru Polski”.

img_2278

Praca doktorska była podsumowaniem 4-letnich studiów doktoranckich, podczas których realizowałem cały temat. Oczywiście najważniejszym celem studiów doktoranckich jest zdobycie stopnia naukowego doktora, ale równie ważne jest wykorzystanie tego czasu na wyrabianie sobie “naukowej marki”. Dowodem wyników pracy w nauce są publikacje w czasopismach naukowych (recenzowanych czasopismach z tak zwanej listy filadelfijskiej). W 4 lata zostałem współautorem 8 takich publikacji.

Drugą połowę czwartego roku studiów doktoranckich poświęciłem na zebranie wszystkich zagadnień, nad którymi wcześniej pracowałem w spójną całość – 180 stron pracy doktorskiej. Moją pracę doktorską możecie pobrać z mojej prywatnej strony: http://marcinpolkowski.com/en/ (po prawej stronie w boksie “Education”).

Należy też zaznaczyć, że jednym z ważnych elementów studiów doktoranckich jest dydaktyka – ja uczyłem głównie programowania w pythonie i C++.

Samo przygotowanie dysertacji (wiki) to dopiero początek drogi do uzyskania stopnia doktora. Trzeba otworzyć przewód doktorski, powołać komisje na dwa (lub trzy jak ktoś nie ma certyfikatu z angielskiego) egzaminy, wyznaczyć dwóch recenzentów i powołać komisję, która przyjmie obronę pracy. Każdy z tych elementów przechodzi przez Radę Wydziału. Gotowa praca trafia najpierw do recenzentów, którzy ją opiniują i dopuszczają do dalszego postępowania – mają na to dwa miesiące (recenzje mojej pracy: recenzja prof. Pietsch, recenzja prof. Guterch). Po otrzymaniu recenzji komisja ds. obrony sprawdza wszystkie dokumenty i wyznacza termin obrony. Moja obrona miała miejsce 17 października 2016.

W celach głównie pamiątkowych postanowiłem dyskretnie zarejestrować całe wydarzenie. Małą kamerę sportową schowałem w piórniku koło komputera, żeby nikogo nie kuła w oczy. Całoś trwała równo dwie godziny. Myślę, że możecie być ciekawi jak taka obrona wygląda (moja prezentacja startuje około 8-9 minuty):

Gdybyście mieli pytania o studia doktorancie, proces uzyskiwania stopnia czy samą moją pracę to chętnie będę odpowiadał – komentujcie poniżej!

Nie zapomnijcie subskrybować mojego kanału na YouTube:

| Komentarze (4)

Czego nie wiemy o liczbach zmiennoprzecinkowych

Przez ostatnie 4 lata prowadziłem na Wydziale Fizyki ćwiczenie do kilku kursów programowania w pythonie i c++. Kilka razy spotkałem się z pytaniem o liczby zmiennoprzecinkowe. Poruszamy na zajęciach kwestię zapisu binarnego liczb całkowitych i wtedy pojawia się pytanie o liczby rzeczywiste. Do tej pory odpowiadałem bardzo skrótowo, że “jest to skomplakowne”, bardzo sprytne, ale nie pozbawione wad. W przypadku liczb całkowitych sprawa jest prosta. W danej ilości bitów możemy przechować dany zakres liczb całkowitych. Z liczbami zmiennoprzecinkowymi problem polega na tym, że w każdym przedziale jest ich nieskończenie wiele, a my mamy tylko 32 albo 64 bity. Więc jak to właściwie jest? O samej formie zapisu można oczywiście przeczytać wszystko w Wikipedii – ja chcę zrobić eksperyment. Na pierwszy ogień liczby typu float, czyli 32 bitowe (1 bit na znak, 8 bitów na wykładnik i 23 na mantysę). 32 bity daje jakieś 4 294 967 296 kombinacji – to na tyle mało, że mogę sprawdzić wszystkie. W ramach testu sprawdzam ile liczb mieści się w 4 przedziałach: 0-1, 1-2, 1000-1001 i 1000000-1000001:

int main()
{
	float f;
	uint i;	 
	std::vector range = {0, 1, 1000, 1000000};
	std::vector count = {0, 0, 0, 0};
	for(i = 0; i=range[a] && f < range[a]+1)
				count[a]++;
		}
	}
	for(int a = 0; a < range.size(); a++)
	{
		std::cout << range[a] << "-" << range[a]+1 << ":\t" << count[a] << std::endl;
	}
}

Wynik mamy po niecałych 2 minutach. Policzone ilości liczb są bardzo ciekawe:

screenshot_33

Od razu widać kilka zależności. W zakresie od 1 do 2 jest dokładnie 2^23 liczb - tak jakby były to wszystkie "mantysy" dla jednego wykładnika. W zakresie od 0 do 1 jest dokładnie 2^23*(2^7-1)-1 - czyli wszystkie mantysy dla połowy wykładników (bez jednego). Im większa liczba tym mniejsza dokładność - mamy do dyspozycji 2^31 a musimy móc pokazać zarówno bardzo małe liczby i te gigantyczne.

Z liczbami w pojedynczej precyzji sprawa była prosta - udało się sprawdzić wszystkie kombinacje w niecałe 2 minuty. Przy podwójnej precyzji mamy 64 bity, więc sprawdzenie wszystkich zajęłoby 16000 lat. Szkoda czasu. Możemy natomiast pobawić się takim kodem:

int main()
{
	union
	{
		uint64_t input;
		uint64_t   output;
	} data;

	uint64_t x = 5;
	data.input = x;
	std::bitset   bits(data.output);

	for(uint64_t a = 0; a<=2048; a++)
	{
		for(uint64_t b = 0; b<=pow(2,52)-1; b+=pow(2,52)-1)
		{
			x = b+(a<<52);
			data.input = x;
			bits = data.output;
			double f;
			memcpy(&f, &x, sizeof(f));
			std::cout << std::setw(4)<< a<< '\t'<< std::setw(16)<<  b << '\t'<< std::setw(12) << f << '\t' << bits  << std::endl;
		}
		std::cout << std::endl;
	}
}

Powyższy kod wypisze zakres liczb dla każdego z 2^11 wykładników. Mantysa ma 52 bity (~4.5*10^15), więc od 1 do 2 mamy właśnie 2^52 liczb, Tyle samo od 2 do 4, 4 do 8 itd:

screenshot_34

Na powyższym przykładzie pierwsza kolumna to wykładnik, druga to mantysa, trzecia to wartość liczby, a czwarta to reprezentacja binarna.

Nie chcę tutaj wnikać w szczegóły samego formatu. Zależy mi na pokazaniu specyfiki takiego pamiętania liczb. Możecie modyfikować i uruchamiać powyższe kody, a dowiecie się więcej o liczbach zmiennoprzecinkowych. Dodatkowo powyższe przykłady powinny dobrze obrazować kolosalną różnicę pomiędzy pojedynczą i podwójną precyzją.

Musimy pamiętać, że operując na bardzo dużych liczbach, gdzie dokładność zapisu zmiennoprzecinkowego jest mała, możemy wprowadzać błędy numeryczne - w szczególności gdy łączymy wielkie i małe liczby.

Kto dowiedział się czegoś nowego klika łapkę w górę!

| Komentarze

Karta dźwiękowa jako przetwornik A/D dla sejsmografu

Dostałem dziś mailem pytanie dotyczące wcześniejszych wpisów dotyczących budowy prostego czujnika sejsmicznego (Sejsmograf domowej roboty – pierwsza przymiarka oraz Sejsmograf domowej roboty – krok drugi):

Dzień dobry,

Jestem studentem geofizyki na AGH i instruktorem harcerskim w ZHP. Przygotowuję zajęcia dla młodzieży popularyzujące sejsmologię. Chcieliśmy na nich zbudować sejsmograf według Pańskiego pomysłu, jednak jedna rzecz jest dla mnie niejasna – jak dokładnie należy podłączyć cewkę do karty dźwiękowej w komputerze? Niestety, nie znalazłem odpowiedzi na to pytanie w artykule. Będę bardzo wdzięczy za odpowiedź.

Z wyrazami szacunku,
Marek Ziobro

Każdy taki mail jest dla mnie bardzo ważny, bo oznacza, że treść, którą tworzę jest wartościowa. Dziś postanowiłem odpowiedzieć publicznie na blogu – myślę, że odpowiedź może się komuś jeszcze przydać.

Cewka zastosowana w sejsmometrze to zwykły zwój drutu mający dwa końce. W wyniku zmiany pola magnetycznego spowodowanej ruchem magnesu w cewce indukuje się napięcie, czyli różnica potencjału pomiędzy obydwoma końcami drutu. Wartość tego napięcia jest proporcjonalna do tego, jak szybko zmienia się pole magnetyczne – im szybciej rusza się magnes tym większe napięcie. Do pomiaru stałego napięcia można oczywiście użyć multimetru, natomiast w naszym przypadku zależy nam na pomiarze zmian napięcia w czasie. Pomiaru napięcia w funkcji czasu możemy dokonać za pomocą oscyloskopu. W przypadku, gdy nie mamy pod ręką oscyloskopu można ratować się kartą dźwiękową komputera z wejściem mikrofonowym. Sygnał z mikrofonu jest napięciem zmiennym w czasie dokładnie tak jak w przypadku czujnika sejsmicznego. Niektóre karty dźwiękowe w komputerach stacjonarnych mają również wejścia zwane line-in, które działają bardzo podobnie.

Oczywiście karta dźwiękowa nie będzie doskonałym oscyloskopem (i tym bardziej rejestratorem sejsmicznym) ponieważ:

  1. Ma ograniczony zakres napięć wejściowych – to bardzo ważne, bo przekroczenie maksymalnego dopuszczalnego napięcia może kartę uszkodzić. Dokładne wartości należy sprawdzić w dokumentacji danego modelu karty.
  2. Pracuje w domyślnym dla dźwięku próbowaniu (44.1, 48kHz) jest to zdecydowanie za dużo, ale lepiej za dużo niż za mało.
  3. Ma spore szumy, co jest istotne przy rejestracji słabych sygnałów.
  4. Ma zazwyczaj niską rodzielczość.

Do samego podłączenia będziemy potrzebowali wtyk jack 3.5mm (mono lub stereo). Do masy podpinamy jeden koniec drutu cewki, a do plusa (w przypadku wtyku stereo do kanału L lub P) drugi. Zamiany miejscami kabli od cewki będzie miała wpływ na kierunek pomiaru: dodatnie napięcie przy ruch w górę lub w dół – na początek zabawy nie ma to żadnego znaczenia, przy głębszej interpretacji warto to wiedzieć i kontrolować. Czasami łatwiej będzie kupić kabel zakończony złączem jack (np. przedłużacz do słuchawek, lub kabel do podłączania odtwarzacza do radia), przeciąć i połączyć kabelki. W przypadku kabla stereo trzeba uważać, żeby nie połączyć cewki do kanałów L i P, bo wtedy nie zadziała – jeden koniec musi być do masy. Ja stosowałem właśnie taki przerobiony kabel:

Na komputrze można użyć dowolnego oprogramowania do rejestracji dźwięku lub specjalnego programu udającego oscyloskop: Zelscope lub Soundcard Oscilloscope. Drugi z nich widać na moich starych zdjęciach (to zdjęcie jest 2006, dlatego w warsztacie miałem monitor CRT):

Pozdrawiam i jak zawsze czekam na Wasze maile!

Zapraszam w środę o 20 na nowy wpis – w tym tygodniu będzie o liczbach zmiennoprzecinkowych.

| Komentarze

Rekrutacja do dużych firm na stanowiska IT – zadania rekrutacyjne

Za niecałe dwa tygodnie będę bronił doktorat i skończy się moja przygoda ze studiowaniem. Zaproponowano mi pozostanie na uczelni w charakterze pracownika naukowego, ale mimo tego – trochę z ciekawości – postanowiłem wysłać swoje CV do kilku firm z branży IT w Polsce i w USA. Praca w USA to dla mnie pewien cel na przyszłość. Na razie rozmowy o ew. przeprowadzce zostały odsunięte w czasie, żeby choć trochę odchować potomstwo (a druga córka w drodze).

Google

Pierwsza przygoda z poszukiwaniem pracy związana była z polskim oddziałem Google. Wysłałem CV i po ok. tygodniu odezwał się rekruter z Londynu. Pierwsza rozmowa z rekruterem trwała ok. pół godziny. Po kilku dniach zaproponowano mi termin rozmowy technicznej. Pozwolono mi wybrać język programowania (C++, java lub Python, którego wybrałem). Rozmowa odbyła się przez telefon (wyświetlił mi się numer ze Szwajcarii) i trwała 40 minut. Swoją odpowiedź tworzyłem w specjalnie do tego celu przygotowanym dokumencie google, który w momencie zakończenia rozmowy stał się tylko do odczytu i chwilę później zniknął. Pierwsze zadanie to prosta łamigłówka (nie trzeba było nic kodować):

Bawisz się z dziećmi w rzuty monetą (do dyspozycji jest wiele monet). Monety nie są idealne i prawdopodobieństwo reszki wynosi 3/7 a orzełka 4/7. Zaproponuj dzieciom najprostszy sposób uzyskania “mniej więcej” uczciwego wyniku losowania (czyli takiego, który wypada z prawdopodobieństwem ok. 0.5).

Po powyższym dostałem zadanie z właściwego programowania i muszę przyznać, że zajęło mi chwile zrozumienie o co chodziło autorowi. Teraz jak patrze z perspektywy czasu treść jest oczywista a rozwiązanie banalne. Zadanie brzmiało:

Twój kolega prowadzi nieskończony eksperyment fizyczny, gdzie mierzona jest wartość liczbowa. Napisz dla niego program w pythonie, który umożliwi ciągłe obliczanie średniej kroczącej z tych pomiarów. Zaproponuj interfejs pomiędzy eksperymentem a twoim programem.

Po około tygodniu otrzymałem telefon, że proces rekrutacji nie będzie kontynuowany.

Tesla Motors

Druga próba była związana ze stanowiskiem inżyniera oprogramowania map i nawigacji w Tesla Motors. Tutaj cały kontakt był przez e-mail a czas pomiędzy wiadomościami liczony był w tygodniach. Po 3 mailach zostałem poproszony o rozwiązanie zadania programistycznego (Tesla Coding Challenge). Zadanie zostało mi udostępnione do pobrania. Miałem tydzień na zabranianie się za rozwiązywanie, ale od pobrania do odesłania zadania czas miał nie przekroczyć 4h (rekomendowane 3h). Rozwiązanie miało być w C++. Zadanie brzmiało:

Napisz program działający w wierszu poleceń, którego zadaniem będzie rozwiązanie układu równań zapisanego w pliku tekstowym. Uproszczenia polegały na tym, że układ ma mieć jedno rozwiązanie, po lewej stronie każdego równania jest tylko jedna zmienna, a po prawej stronie każdego równania mogą być inne zmienne i liczby całkowite dodatnie. Jedyną operacją po prawej stronie może być dodawanie. Zadaniem programu jest rozwiązanie układu i wypisanie zmiennych z wartościami w kolejności alfabetycznej.

Tak dział mój program:

screenshot_31

Po trzech tygodniach dostałem maila, że doszli do wniosku, że potrzebują kogoś od zaraz, a mnie trzeba by “sprowadzić”.

HSBC

Wysłałem CV również na stanowisko specjalisty ds. oceny ryzyka w banku – było to jedno z niewielu stanowisk, gdzie oczekiwano znajomości matlaba. Praca niestety w Krakowie, ale zostałem zaproszony na rozmowę przez telekonferencję do siedziby w Warszawie. Zostałem poproszony o próbki mojego kodu w matlabie, a że nie mam nic, co nie byłoby tajemnicą mojego dotychczasowego pracodawcy poproszono mnie o rozwiązanie zadania:

1. Zapytanie użytkownika o numer telefonu.
2. Na klawiaturze telefonu kazda cyfra, ma odpowiadające jej kilka
literek (1: ABC, 2:DEF, 3:GHI, …, 9:XYZ, 0:spacja). <- to sobie sam zmieniłem na układ zgodny z faktyczną klawiaturą telefonu
3. Utworzyć w Matlab skrypt, który będzie zamieniał cyfry, na
odpowiadające im litery. (Każdy numer wiele wariantów utworzonego słowa).
4. Z utworzonych wariantów słów wybrać 3 ciekawe – np. poprzez
podłączenie do google i wybranie 3, które zwracają najwięcej wyników
wyszukiwania (lub podłączenie do jakiegokolwiek innego źródła danych –
słownika czy tym podobne i wybranie słów istniejących).

Nie wiem czy dostali moje rozwiązanie, bo wrzuciłem w załącznik duży plik ze słownikiem sjp.pl, ale już się nie odezwali.

Nie publikuję na razie rozwiązań – pobawcie się sami, kiedyś pokażę swoje. Tak na prawdę najbardziej szkoda mi było pracy w google. Ani do USA, ani (tym bardziej) do Krakowa w obecnej sytuacji nie zdecydowałbym się jechać całą rodziną. Google był na miejscu. Niestety była to moja pierwsza rekrutacja – nie miałem pojęcia czego się spodziewać i poległem na banalnym zadaniu. Zachęcali mnie, żeby spróbować w przyszłości i pewnie się jeszcze skuszę. Na razie mam nową pracę.

| Komentarze (3)

Ostatni dzień września z drona: nowy budynek IGF

Udało mi się zrobić dziś kilka ładnych zdjęć z drona. Najpierw widok na Warszawę z okolic Mostu Północnego:


003_sm004_sm

Następnie widok na nowy budynek Instytutu Geofizyki Wydziału Fizyki UW:


001_sm002_sm

Wokół powyższego budynku nagrałem również kilka ujęć wideo. Moje ujęcia nie są niestety gładkie jaku u @CaseyNeistat, ale ćwiczę. Zapraszam, to tylko 54 sekundy:

| Komentarze

27 wyjazdów na pomorze, 31000 km zapisane w GPS

O eksperymencie sejsmicznym “13 BB Start” już wspominałem (obiecuję napisać więcej). W ramach eksperymentu przez 3 lata w 13 miejscach na pomorzu umieszczone były szerokopasmowe stacje sejsmiczne. Ponieważ nasza siedziba jest w Warszawie trzeba było trochę na pomorze jeździć: zaplanować, zainstalować, serwisować i zdemontować. Ostatni wyjazd był w zeszłym tygodniu. Ponieważ od początku jeździłem prywatnym samochodem potrzebowałem rejestrować za pomocą GPS przebytą trasę, żeby mieć podkładkę do rozliczenia kosztów. Wykorzystałem do tego dwa bliźnicze odbiorniki: starszy Garmin GPSmap 60CSx i młodszy GPSmap 62sc. Podczas wyjazdów odbiornik zawsze był w samochodzie, uruchamiał się po uruchomianiu silnika i zapisywał pozycję do pliku co 1 sekundę.

W sumie wyjazdów było 27, a całkowita odległość przekroczyła 31 tysięcy kilometrów. Zapis tras obejmuje ponad półtora miliona punktów, co daje ponad 400 godzin jazdy (średnio 16 godzin na wyjazd). Jest to piękny zbiór danych, który wymaga wizualizacji na mapie. Oczywiście nie byłbym sobą, gdybym najpierw nie wczytał wszystkich danych do Matlaba (zawsze po wyjeździe wczytywałem zapis do Matlaba i skryptem generowałem raport wszystkich odcinków tras: skąd – dokąd i od kiedy – do kiedy). Goła wizualizacja danych jest bardzo ładna, ale dopiero naniesienie tras na mapę daje pełny obraz sytuacji. Wykorzystałem do tego oczywiście mapy google. Z Matlaba zapisałem moje ścieżki do pliku KML (pamiętając o wpisaniu NaN na końcu każdego odcinka, żeby podczas rysowania była przerwa pomiędzy odcinkami) i ten wgrałem na mapę Google (po zalogowaniu idziemy w moje miejsca i moje mapy, gdzie można wgrać plik KML). Plik ze wszystkimi moimi danymi miał rozmiar ponad 40 mega, więc żeby zaspokoić ograniczenia google (max 5mb) podczas zapisu brałem co 10 punkt. Efekty wyglądają tak:

screencapture-google-maps-d-u-2-edit-1474983786479

Zbliżenie na obszar startowy i docelowy:

screencapture-google-maps-d-u-2-edit-1474983297187

screencapture-google-maps-d-u-2-edit-1474983254545

Trzeba zaznaczyć, że podczas pierwszych naszych wyjazdów (pierwszy był w 2012), autostrada A1 była jeszcze w budowie na odcinku z Łodzi do Torunia. W tym miejscu przyszedł mi do głowy pewien pomysł. Ponieważ ponad połowa moich wyjazdów była właśnie przez autostradę A2 i A1 zrobiłem prostą analizę. Dla każdego wyjazdu tą drogą brałem czas przejechania przez punkt wspólny (zaraz za połączeniem S8 i A2 pod Warszawą, pomarańczowy znacznik na mapie) i zaznaczyłem położenie po 30 (zielone), 60 (brązowe), 120 (żółte), 180 (fioletowe) i 240 minutach (jasno zielone). Oczywiście jeździłem o różnych porach, czasem z tankowaniem, w różnych warunkach, itp., natomiast ładnie widać, że im dalej tym większy rozrzut:

screencapture-google-maps-d-u-2-edit-1475009898800

Zupełnie przy okazji pisania tego wpisu zakochałem się ponownie w mapach google. Dodatkowym, ważnym atutem takiej mapy google, jest możliwość wczytania mapy z nałożonymi trasami w aplikacji w telefonie, co znacząco ułatwia nawigację np. powtórzenie trasy po lesie. Postaram się kiedyś więcej na ten temat napisać.

Bonus! Jeden z przejazdów nagrałem w całości kamerą sportową (takie GoPro, ale od Sony):

| Komentarze

Wyjazd pomorze z lotu ptaka [4k]

Nie zapomnijcie subskrybować mojego kanału na YouTube:

| Komentarze

Wpływ sezonowych zmian w jonosferze na pomiary GPS

Od czterech lat biorę udział w pasywnym eksperymencie sejsmicznym “13BB star”. Będę na ten temat pisał jeszcze wiele razy, ale muszę tutaj zrobić dwa zdania wstępu. Celem eksperymentu jest rejestracja naturalnych zjawisk sejsmicznych, co pozwala na określanie głębokich struktur w Ziemi. W ramach eksperymentu na pomorzu pracowało 13 samodzielnych stacji poukrywanych głęboko w lasach. Każda ze stacji wyposażona była w odbiornik GPS, których celem było dostarczanie bardzo precyzyjnych informacji o czasie. Gdzieś na boku i przy okazji raz na godzinę każda stacja zapisywała pozycję według gps: długość i szerokość geograficzną oraz wysokość n.p.m.

Któregoś dnia postanowiłem wczytać do matlaba te dane dla jednej ze stacji. Gdy narysowałem wykres zależności mierzonej wysokości n.p.m. od czasu zainteresowała mnie mało widoczna zmiana sezonowa. Zupełnie nie dziwi rozrzut pomiarów o +- 10 metrów, bo pomiary “cywilnym” (jako przeciwieństwo do geodezyjnego) GPS są obarczone mniej więcej takim błędem. Zaciekawiło mnie sezonowe przesunięcie pomiarów, które jest bardzo dobrze widoczne po narysowaniu tygodniowej średniej kroczącej:

out_13bb_gps_001

Sprawdziłem od razu czy podobny efekt widać na pozostałych stacjach:

out_13bb_gps_002

Jeszcze ciekawej sprawa wygląda gdy porówna się średnie kroczące dla wszystkich stacji poprzez odjęcie średniej z całego okresu. Poniżej średnie 7-dniowe:

out_13bb_gps_003_7

30-dniowe:

out_13bb_gps_003_30

i 60-dniowe:

out_13bb_gps_003_60

Na powyższych wykresach widać wyraźną sezonowość zmian. Poprosiłem o komentarz specjalistę od pomiarów referencyjnych GPS: “jonosfera, standard, nuda”. Natomiast dla mnie ciekawe jest, że skoro efekt jest tak dobrze znany i tak powtarzalny, czemu nie została wprowadzona na to poprawka. To właśnie dynamika zmian jonosfery powoduje, że cywilne GPS mają ograniczoną dokładność. Sprzęt geodezyjny wykorzystuje stacje permanentne, dzięki którym można określi wpływ jonosfery na pomiar w danym momencie (stacja permanentna się nie rusza, więc odchylenie pomiaru wynika właśnie ze zmian w jonosferze). Takich szybkich lokalnych zmian nie da się przewidzieć, ale zmiany między latem a zimą wydają się całkowicie powtarzalne i oczywiste.

Gdyby ktoś był zainteresowany kodem użytym do wygenerowania 3 dolnych rysunków:

%% 7 dni
clf; hold on
files = dir('*.gps');

colors = lines(length(files));

for fi = 1:length(files)
	X = load(files(fi).name);
	D = X(:,1) - circshift(X(:,1),[1 0]); 
	idx = [1; find(D>1); length(X(:,1))]; % Szukanie przerw w danych większych niż 1 dzień
	X(X(:,4)>500,4)=NaN; 
	X(X(:,4)<50,4)=NaN; % Usuwanie pomiarów "od czapy"
	for a = 1:length(idx)-1
		i = idx(a):idx(a+1)-1;
		i(isnan(X(i,4)))=[];
		f = nanmean(X(:,4)); % Średnia dla całego okresu
		if length(i) > 7*24 % Badamy tylko fragmenty o długości > 7 dni
			% Obliczenie i narysowanie średniej kroczącej przesuniętej o średnią globalną
			plot(X(i,1), tsmovavg(X(i,4),'s',7*24,1)-f, 'Color', colors(fi,:), 'LineWidth',3)
		end
	end
end


% Rysowanie i podział rysunku na pory roku

axis([datenum(2013,09,01) datenum(2016,07,31), -6 4])
pory = [datenum(2013,09,22) datenum(2013,12,21) datenum(2014,03,20) datenum(2014,06,21) datenum(2014,09,23)  datenum(2014,12,21) datenum(2015,03,20) datenum(2015,06,21) datenum(2015,09,23) datenum(2015,12,22)  datenum(2016,03,20) datenum(2016,06,20) datenum(2016,09,22)];

names = {'Wiosna', 'Lato', 'Jesien', 'Zima'};
for a = 1:length(pory)-1
	text(mean(pory(a:a+1)),3.5,names{rem(a+1,4)+1},'HorizontalAlignment', 'Center', 'FontSize', 15);
	plot([pory(a) pory(a)], [-10 10],'k--')
end

set(gca,'xTick', pory)

datetick('x','YYYY-mm-dd','keepticks')

hgexport(gcf, '13BB_GPS_003_7.png', hgexport('factorystyle'), 'Format', 'png', 'resolution', 300, 'Width', 40, 'Height', 10);

| Komentarze

Burze w Polsce, edycja 2014

Pod koniec 2013 roku wpadłem na pomysł, aby regularnie zapisywać sobie obrazek (mapę) przedstawiający aktualne lokalizacje wyładowań atmosferycznych. Korzystałem wtedy ze strony blide.de. Oczywiście nie miałem na myśli robienia tego ręcznie, więc uruchomiłem na serwerze skrypt, który równo co 15 minut zapisywał aktualny obrazek. Skrypt działa do dziś i cały czas zapisuje aktualne mapy. Przypomniałem sobie o nim, gdy na serwerze zaczęło brakować miejsca na dane. Pobrałem obrazki na swój komputer i zacząłem zabawę. Na pierwszy ogień poszedł pierwszy pełny rok – 2014.

Teoretycznie w ciągu 365 dni powinienem zapisać 365*24*4 = 35040 obrazków. Faktycznie zapisałem 34999, czyli całkiem nieźle. Pierwsza, oczywista myśl to animacja. Połączyłem wszystkie obrazki w filmik:

W kolejnym korku chciałem podejść do tematu bardziej kompleksowo. Jak oszacować ilość wyładowań dla każdej z map. Oczywiście ilości się nie da, ale da się oszacować obszar objęty burzą. Na mapach kolorem czerwonym oznaczone są miejsca, gdzie wystąpiło co najmniej jedno wyładowanie w przeciągu ostatnich 15 minut. Niestety pliki są zapisane w mocno skompresowanym jpg, ale wyraźnie widać co czerwone a co nie. Napisałem na kolanie szybki skrypt w Matlabie:

clear;clc
files=dir('2014/blids/*.png');
S = zeros(length(files), 2);
for a = 1:length(files)
	name = ['2014/blids/' files(a).name];
	C = strsplit(name, '/');
	C2 = strsplit(C{end},'.');
	date = C2{1};
	time = C2{2};
	I = imread(name);
	idx = (I(:,:,1)>100 & I(:,:,2)<100 & I(:,:,3)<100);
	S(a,1) = sum(sum(idx));
	S(a,2) = datenum(str2double(date(1:4)), str2double(date(5:6)), str2double(date(7:8)), str2double(time(1:2)), str2double(time(3:4)), 0);
	
	fprintf('%d\t%d\n',a,S(a,1));
end

Dla każdej mapy skrypt obliczył i zapisał ilość pikseli, w których miały miejsca wyładowania. Teraz wystarczy to ładnie nanieść na wykres: np. sumarycznie dla każdego dnia.

clf
days = datenum(2014,01,01):datenum(2014,12,31);
storm = zeros(size(days));

for d = 1:length(days)
	storm(d) = sum(S(S(:,2)>=days(d)&S(:,2)0)=log10(storm(storm>0));
stairs(days-datenum(2014,01,01),storm, 'LineWidth',2)
axis([0 365 0 2.5e5])
hold on
for month = 1:12
	plot(ones(1,2)*datenum(2014,month,1)-datenum(2014,01,01),[0 2.5e5], 'k--')
end

A oto efekt:

001

Pionowe linie oddzielają kolejne miesiące. Widać wzmożoną aktywność burzową od kwietnia do września. Maksimum przypada na początku sierpnia, co widać na filmie.

Jak wykorzystalibyście takie dane? W wolnej chwili zrobię animację i wykres dla 2015 roku.

| Komentarze

« Nowsze wpisy - Starsze wpisy »