Geologisch onderzoek met de Pi en Mathematica
Geologisch onderzoek met de Raspberry Pi en Mathematica
Mathematica kan niet alleen vergelijkingen doorrekenen en oplossen, maar biedt ook toegang tot Wolframs wetenschappelijke database. Je mag Mathematica op een Raspberry Pi gratis gebruiken – bijvoorbeeld om aardbevingen inzichtelijk te maken.
Wolframs wiskundige software Mathematica is gratis beschikbaar op de Pi-minicomputer voor privé- en educatief gebruik. Bij Mathematica hoort de krachtige wetenschappelijke online database Wolfram Knowledgebase. Die stelt continu bijgewerkte data beschikbaar uit allerie onderzoeksgebieden. Met behulp van onderwerpspecifieke query-functies krijg je uitsluitsel over ruimte-objecten van satellietbrokstukken tot melkwegclusters, over moleculen van H2O tot DNA, en er zijn allerlei details over plaatsen van Druten tot en met Los Angeles. De database levert bovendien tienduizenden actuele records over geregistreerde aardbevingen. Elk van die datarecords bevat onder meer het tijdstip, de geografische coördinaten en de intensiteit van een aardbeving.
Aan de hand van die aardbevingsdata leer je hier hoe je met de Wolfram Language een query aanmaakt, die naar de Knowledgebase stuurt en de resultaten vervolgens grafisch weergeeft.
Zo kun je het voorkomen van aardbevingen in bepaalde regio's en de spreiding qua tijd en intensiteit goed weergeven, bijvoorbeeld voor een lezing over de laatste ontwikkelingen op het Indonesische vakantie-eiland Lombok. Als je wat programmeerervaring hebt, krijg je vanzelf een indruk hoe het voorbeeldprogramma in dit artikel werkt. Voor nieuwelingen is het wellicht raadzaam de tutorials van Wolfram eens door te nemen. Die staan, net als het voorbeeldprogramma, vermeld bij de link aan het eind van dit artikel.
Raspbian heeft Mathematica als grafische X-toepassing beschikbaar voor de standaardgebruiker pi. Voor andere scenario's kun je het met het pakketmanagement van Raspbian installeren met sudo apt install wolfram-engine.
Universeel notebook
Na het starten van Mathematica verschijnt er op het beeldscherm een leeg notebook – zo heten projectbestanden in Mathematica. Daar typ je de code in, voer je hem uit en krijg je de resultaten te zien. Dat kan tekstoutput zijn, maar ook grafieken, animaties en zelfs geluiden. Bij het hier beschreven zoeken naar aardbevingen opent de code bij het uitvoeren een eigen notebook voor de grafische weergave van de resultaten.
Notebooks zijn opgebouwd uit cellen, die ook genest mogen zijn. Mathematica markeert het bereik waarover een cel zich uitstrekt met een rechthoekige haak aan de rechter vensterrand. Om een cel te selecteren, klik je op die rechthoekige haak of in het celbereik. Met Shift+Enter evalueert Mathematica de inhoud van de geselecteerde cel en laat het resultaat zien. Met Alt+. breek je een lopende berekening af.
Omdat het notebookformaat bedoeld is voor presentaties en desktop-publishing, kunnen cellen behalve code en resultaten ook titels, commentaren en dergelijke bevatten. Een selectie van de alle beschikbare celtypen krijg je via de rechter muisknop.
Anders dan de meeste programmeertalen accepteert Mathematica zonder te protesteren ook onbekende namen voor variabelen en functies. Elke expressie wordt net zo ver uitgevoerd als mogelijk is: als x niet bekend is, wordt 3+4+x dan 7+x. Dat is niet zo als een formule syntaxisfouten heeft zoals niet afgesloten haakjes. Dan is de expressie niet compleet en krijg je een foutmelding.
Als je aardbevingen.nb downloadt en opent met Mathematica, verschijnt een notebook met één enkele cel. Die bevat de hele code van deze applicatie. Voer deze cel uit met Shift+Enter, dan springt de weergave naar het einde van de cel en krijg je invoervelden te zien waar al standaardwaarden in zijn ingevuld. Daar kun je de gegevens invullen voor het centrum en de straal van het gewenste zoekgebied, evenals de begin-en eindtijdstippen van de gewenste periode en een minimale kracht (intensiteit) van de aardbevingen die meegenomen moeten worden. Bovendien kun je een naam invullen waarmee de grafische weergave in notebookformaat moet worden opgeslagen.
Een klik op 'Start zoeken' voert dan zowel de zoektocht in de Knowledgebase uit als de aansluitende verdere verwerking van de data. Bij de invoervelden zie je bij het tekstveld 'Status' wat de voortgang van de lopende berekening is. Tot slot verschijnt een
ander notebook met de grafische weergave van de resultaten. De duur van de query hangt sterk af van de gewenste periode waarbinnen de aardbevingen gezocht moeten worden en hoe zwaar de Knowledgebase momenteel belast wordt. In ongunstige gevallen kan het een paar minuten duren.
Rondleiding door de code
In de Knowledgebase wordt naar posities op het aardoppervlak verwezen op basis van hun geografische coördinaten. Om bijvoorbeeld de lengte- en breedtegraad van Tokyo op te vragen, kun je de ingebouwde functie Interpreter gebruiken. Die zet gewone namen om in een voor Mathematica begrijpelijk formaat. Om precies te zijn: een GeoPosition-object. Als je die functie aanroept met het argument "Location" (voor plaats) en bijvoorbeeld een plaatsnaam als Tokyo, dan krijg je de bijbehorende coördinaten in de uitvoer:
Interpreter["Location"]["Tokyo"]
Dan verschijnt de uitvoer:
GeoPosition[{35.67,139.77}
Ook invoer als "Eiffel Tower" en "Machu Picchu" is mogelijk.
De voorbeeldcode definieert de functie getPositionFromString, waarmee het aanroepen van Interpreter["Location"] met een meer beschrijvende functienaam mogelijk is:
getPositionFromString[input_] :=
Interpreter["Location"][input]
Daarbij is input de placeholder voor de plaatsaanduiding. Aan de linkerkant van die aanduiding moet een underscore (_) volgen. Die geeft aan dat het om een zogeheten pattern gaat. Meer daarover staat in de hulpfunctie van Mathematica bij 'Defining Functions' en 'Transformation Rules for Functions'. Voor het vervolg is het alleen van belang om te weten dat op de parameternamen een underscore moet volgen, maar niet op de identifier rechts van de toewijzingsoperatie :=. Daardoor levert bijvoorbeeld
getPositionFromString["Eiffel Tower"]
de waarde GeoPosition[{48.8583,2.29444}] op. Om de periode in te perken waarbinnen naar de aardbevingsdata gezocht wordt, is het aangeven van kalenderdatums in een eigen formaat van Wolfram Language in plaats van als string nodig. Ook dat kan de Interpreter doen, en wel door het aangeven van "Date" in plaats van "Location". De voorbeeldcode definieert voor een betere leesbaarheid een functie met de aanroep:
getDateFromString[input_] :=
Interpreter["Date"][input];
Door het intypen van
getDateFromString["1.1.2018"]
verschijnt in de output een klein kalenderpictogram en daarnaast 'Day' en 'Mon 1 Jan 2018'. Die output komt van Mathematica. Als je wilt weten welke Wolfram Language-code daarachter zit, klik je daaronder op 'date formats/text'. Dan verschijnt
DateString[DateObject[{2018,1,1}, "Day",
"Gregorian", 2.]]
Wat de betekenis is van de andere elementen achter jaar, maand en dag, kun je bij het Wolfram Language Documentation Center opvragen door erop te klikken.
Data ophalen
Dan moet je de aardbevingsdata opvragen met de Mathematica-functie EarthquakeData. Als doelgebied kun je daar een GeoDisk aan meegeven, die in de huidige Mathematicaversie tot het onderzoeken van een rechthoekig gebied leidt. Als andere filtercriteria kun je aan EarthquakeData de gewenste minimale kracht en de gewenste periode meegeven. De functie getEarthquakeData geeft de parameters door aan EarthquakeData:
getEarthquakeData[center_, radius_, t0_,
t1_, magnitude_] :=
Module[{data, area}, area = GeoDisk[center, Quantity[radius,
"Kilometers"]]; data = EarthquakeData[area,
magnitude, {t0, t1}]; Map[{#["Period"],
#["Magnitude"],
#["Position"]} &,
Values[data]]
]
De Module kun je net als bij andere talen als functie opvatten. Hij zet de tussenberekeningen in area en data, zodat die niet onnodig als symbolen uit de functie naar buiten hoeven te komen. Dergelijke lokale variabelen zet je als eerste parameters tussen accolades.
GeoDisk verwacht als eerste parameter het middelpunt van een cirkel en als tweede een waarde voor de straal. Om de eenheid vrij te kunnen kiezen wordt bijvoorbeeld niet naar een getal in kilometers gevraagd, maar naar een Quantity, waar je een willekeurige eenheid aan mee kunt geven, Zinvol zijn daarbij echter alleen lengtematen als "Kilometers" en "Miles".
EarthquakeData levert een lijst terug van sleutel-waardeparen (Association). Van alle mogelijke sleutels zijn voor de voorbeeldtoepassing alleen de duur ("Period"), de kracht ("Magnitude") en het middelpunt ("Position") van elk resultaat interessant. Om ervoor te zorgen dat getEarthquakeData alleen die waarden per lijstelement terug levert, loopt Map alle waarden van data af en maakt daar een nieuwe lijst van. De # refereert naar de eerste parameter van een pure functie (in andere talen heet dat een lambda- of anonieme functie). Parameters zijn in dit geval de na elkaar doorgegeven elementen van data. Daarbij refereert #2 naar de tweede parameter, #3 naar de derde enzovoorts. Met de & wordt deze constructie afgesloten. Meer daarover bij het onderdeel Pure Functions in de documentatie van Mathematica.
Met de volgende aanroep kun je dan een overzicht krijgen van de laatste aardbevingen op het Indonesische vakantieeiland Lombok:
getEarthquakeData[ getPositionFromString["Lombok"], 100, getDateFromString["4.8.2018"], getDateFromString["6.8.2018"], 4]
In het tweede element van de uitvoer herken je meteen de flinke aardbeving met kracht 6,9 die het noordelijk deel van het eiland om 11:46 MET teisterde:
{{DateObject[{2018, 8, 4, 5, 44}], 4.4,
GeoPosition[{-8.3939, 116.611}]}, {DateObject[{2018, 8, 5, 11, 46}], 6.9,
GeoPosition[{-8.2871, 116.451}]}, {DateObject[{2018, 8, 5, 12, 42},], 4.9,
GeoPosition[{-8.3373, 116.185}]}, ...}
Bij deze output hebben we het weergeven van de waarden in dit geval even beperkt tot de hier relevante waarden.
Dat het epicentrum inderdaad in het noorden van het eiland lag, kun je bijvoorbeeld zien door de coördinaten (in GeoPosition[{-8.2871, 116.451}) in Google Maps in te typen. Dat is echter nogal omslachtig.
Op het scherm tonen
Het voorbeeld-notebook maakt het makkelijker: alle resultaten worden daarbij op een landkaart weergegeven. Dat gebeurt met de functie renderEarthquakeMap:
renderEarthquakeMap[data_] := GeoGraphics[
{Apply[{RGBColor[1, 0, 1, 0.4], PointSize[0.001 (#2^2)], Point[#3]} &, data, {1}]},
GeoRange -> Automatic, GeoScaleBar -> "Kilometers", ImageSize -> 700] De functie GeoGraphics genereert de kaart uit de in de eerste parameter meegegeven geografische kenmerken. In dit geval zijn dat kleine cirkels. RGBColor[1, 0, 1, 0.4] zorgt ervoor dat ze violet worden en een dekking van 40 procent krijgen. De grootte van de cirkels (PointSize) neemt toe met de kracht (#2) tot een grootte van 2,5 keer. Point zorgt er tenslotte voor dat de cirkel op de met #3 geextraheerde geografische coördinaten getekend wordt. De laatste parameter {1} bepaalt de laag waaruit de #2 en #3 de waarden uit data moeten extraheren. Als je die zou weglaten of zou vervangen door {0}, dan hadden ze betrekking op de nulde laag, oftewel het tweede en derde element in data. We willen hier echter de kracht (#2) en de geografische coördinaten (#3) van elk element in data hebben, dus moeten we een laag dieper.
Op de kaart zie je dan vervolgens wel de locatie van de aardbevingen en de kracht, maar niet hoe vaak er aardbevingen van een bepaalde kracht binnen een tijdsbestek optreden. Voor het visualiseren van de frequenties biedt Mathematica het Histogram:
renderMagnitudeHistogram[data_] := Histogram[Map[#[[2]] &, data], {0.1},
AxesLabel -> {"Kracht",
"Frequentie"},
ImagePadding -> 60, ChartElementFunction ->
"FadingRectangle",
ChartStyle -> Blue
]
In de eerste parameter verwacht Histogram de te verwerken waarden, in de tweede de intervalgrootte. Een intervalgrootte van 0,1 geeft aan dat de frequenties in intervallen met stappen van 0,1 geteld moeten worden. De andere optionele parameters bepalen dingen als de aslabels en het uiterlijk van het histogram.
En actie!
Het invoerformulier ontstaat uit verschillende invoervelden (InputField) met beschrijvingen van de velden (Text) die het lay-outelement Grid dan in tabelvorm weergeeft.
Een klik op de met Button gedeclareerde startknop brengt de boel aan het lopen. De tweede parameter daarvan is quasi het hoofdprogramma. Dat wordt bij het klikken op de knop uitgevoerd. De daar uitgevoerde code haalt de waarden uit de invoervelden en geeft die door aan getEarthquakeData om de gewenste data uit de Wolfram Knowledgebase te halen. Het resultaat daarvan wordt op zijn beurt dan weer doorgegeven aan de functies voor het tekenen van de kaart en de diagrammen.
Epiloog
Ondanks de bescheiden rekencapaciteit van de Raspberry Pi is de gratis Mathematicaversie voor deze minicomputer een serieus te nemen tool, die zich leent voor de meest uiteenlopende taken. Een groot voordeel van Mathematica is de toegang tot de Wolfram Knowledgebase – en een in dertig jaar tijd opgebouwd assortiment van duizenden ingebouwde functies. Een ander pluspunt is het feit dat je in de Mathematica-notebooks berekeningen, teksten, animaties en dergelijke in een presentatiewaardige vorm kunt verzamelen – die notebooks kunnen dan ook de basis zijn voor interactieve presentaties en publicaties.
Als je uit de objectgeoriënteerde wereld komt, zul je bij het werken met Mathematica even moeten omdenken – op zich een uitdaging, maar wel een die volgens ons de moeite loont. Als je Lisp, Clojure of Haskell kent of bekend bent met de werking van functionele programmeertalen, zul je je snel thuis voelen.
Als je het voorbeeldprogramma wilt aanpassen, dan nog even een belangrijke tip: voer na een modificatie altijd het commando ClearAll["Global`*"] uit, waarmee alle declaraties en definities verwijderd worden. Anders kunnen overschreven functiedefinities werkzaam blijven, ook al staan ze allang niet meer in de code.