Som mange sikkert har fått med seg gikk Håvard Tveite bort i slutten av mai 2021. Med det mistet Norge en av de beste og mest dedikerte fagpersoner innen GIS, og QGIS mistet en betydelig bidragsyter.
Håvard Tveite var førsteamanuensis på NLH/UMB/NMBU/ hele tiden fra 1992. Han hadde ansvar for utvikling og gjennomføring av en rekke GIS-emner på ulike nivå og bidro sterkt til introduksjonsemnet LAD102 som i dag tas av over 400 studenter fra «alle» studieprogram på NMBU. I tillegg var han en faglig brobygger mot andre miljø på Campus. Håvard var sentral også i oppstarten av GIS-delen av det som i dag er et egen studium i folkehelsevitenskap på LANDSAM.
Håvards betydning for QGIS kommer tydelig frem i forbindelse med QGIS versjon 3.24 som er gitt navn Tisler – den lille øya på Hvaler som var et av Håvards favorittsteder.
I Skaff og bruk lynnedslagsdata fra met.no beskrev jeg hvordan man henter inn lynnedslagsdata og lagrer dem i en tekstfil som kan brukes i QGIS. For å gjøre rutinemessige analyser av lynnedslagsdata, henter vi data fra MET med jevne mellomrom og lagrer dem i en PostGIS-tabell.
create table lightning (
timestamp timestamp with time zone not null,
point geometry(POINT,25833) not null,
peakcurr integer not null,
multiplicy integer not null,
sensors integer not null,
dof integer not null,
angle numeric(7,3) not null,
semimaj numeric(6,3) not null,
semimin numeric(6,3) not null,
chi2 numeric(5,3) not null,
risetime numeric(3,1) not null,
peak2zerotime numeric(6,2) not null,
maksror numeric(6,2) not null,
cloud boolean not null,
angleind boolean not null,
signalind boolean not null,
timingind boolean not null,
ellipse geometry(POLYGON,25833),
primary key (timestamp,point)
);
Scriptet som kjøres for å hente inn data, kan hentes som en gist fra github. Det er en pythonfil som henter data fra Frost, lagrer hvert punkt som en geometri i PostGIS, deretter beregnes usikkerhetsellipsene og lagres i samme tabell. Hver rad i tabellen inneholder da to geometrier. Når vi skal laste inn en slik tabell i QGIS, vil hver geometri komme som en rad i oversikten:
Ligthningtabellen listes to ganger i QGIS, en gang for hver geometritype
Det er da mulig å legge inn to lag i QGIS, et for punktene og et for usikkerhetsellipsene. Det er imidlertid et par problemer med dette. I tillegg til at et datasett vises som to separate lag, er det ikke mulig å samsortere tegning av punkter og usikkerhetsellipser. For eksempel vil man ofte ønske å lage et kart med de nyeste dataene øverst. Med punkter og ellipser i hvert sitt lag, vil da alle punkter komme oppå alle ellipser, eller motsatt.
Det er imidlertid mulig i QGIS å definere et lag hvor en geometri er en bit av dekorasjonen på en annen geometri. Når disse to ligger i samme rad i en PostGIS database, er det da (nesten) så enkelt som å hente ut en feltverdi.
Her blir usikkerhetsellipsen tegnet opp med nedslagspunktet som en dekorasjon
Dersom man ønsker å bruke usikkerhetsellipsen for overlay i QGIS, må dette være datapunktet som plottes i utgangspunktet. Velg lightning raden med polygongeometri. For å legge inn punktet i dette, legg til et lag til i symbolet, sett typen til Geometry generator og skriv inn følgende tekst:
geom_from_wkt( string_to_array(«point»,’;’)[1] )
Punktet vises som en wkt-tekst, men har projeksjonsinfo som et prefix skilt fra wkt-teksten med et semikolon. Derfor er det nødvendig med en string_to_array som deler på «;» og plukker ut det andre elementet. geom_from_wkt lager da punktgeometrien.
Dersom man ønsker å plotte punktet med ellipsen rundt, er det samme fremgangsmåte, men «geometry type» må settes til «Polygon / Multipolygon» og teksten som skrives inn må i stedet være
«point» og «ellipse» i eksemplene over har ikke noe å gjøre med geometritypene, men er henvisninger til kolonnenavnene i PostGIS.
Legg også merke til at dersom feltet er definert som «geography» og ikke «geometry» er ikke projeksjonen med i feltet (i og med at geography håndterer uprojiserte data) så da bruker man geom_from_wkt(«<feltnavn>») direkte uten noe string_to_array.
I «Skaff og bruk lynnedslagsdata fra MET, hentet vi ned lyndata fra Meteorologisk institutt og filtrerte ut så bare registeringene av lynnedslag ble vist. Nå skal vi bruke usikkerhetsellipsene som finnes i datasettet for å se om noe av interesse kan ha blitt truffet av lyn.
NORSAR har noen seismiske stasjoner i Innlandet hvor sensorer over et stort område er koblet sammen med et dedikert kablet nett. Av og til detter en sensor ut, dette kan være forårsaket av et lynnedslag som har truffet kabelene. Da er det interessant å se om dette kan ha skjedd.
Kartet viser et utsnitt av NORSARs kabler sammen med lynnedslag for en valgt periode
Lyndataene har en usikkerhetsellipse, om denne overlapper en kabel, tolkes det som et mulig lynnedslag i kabelen. I QGIS finnes funksjonen «Geometry by expression» under «Vector geometry» i Processing Toolbox. Her kan vi kjøre funksjonen make_ellipse som trenger argumentene «major semi axis», «minor semi axis» og angle som vi kan hente fra lyndatasettet. Det er bare et lite problem, for at dette skal fungere må datasettet være lagret i en projeksjon med meter som enhet. Dataene som er hentet inn fra MET er i geografisk projeksjon med grader som enhet.
Første trinn er å få ønsket projeksjon på datasettet. Åpne Processing – Toolbox og finn «Reproject Layer» under «Vector General». Velg laget med lyndata og en passende projeksjon – gjerne UTM med riktig sone, i mitt tilfelle med data fra Hedmarken velger jeg ETRS89 UTM zone 32N.
«Reprojected» velger hvordan det nye laget skal lagres. Ved å la det stå som [Create temporary layer] lages et midlertidig lag som forsvinner når QGIS avsluttes. Trykk «Run» for å lage det nye midlertidige laget.
Fra lagoversikten i QGIS etter at jeg har kjørt «Reproject Layer» Ikonet i høyre kant ved «Reprojected» viser at det er et midlertidig lag
Da er neste trinn å lage ellipsene. Gå til Processing Toolbox igjen, finn Vector Geometry og Geometry by expression under denne. Velg det reprojiserte laget og skriv inn make_ellipse($geometry,major*1000,minor*1000,angle) i Geometry Expression. (Dette forutsetter at lynlaget er definert som beskrevet i del 1). Aksene må multipliseres med 1000 siden de er oppgitt i km og vi trenger dem i meter. Trykk eventuelt på sigma-knappen ved siden av for å få opp en større formeleditor med mulighet til å velge felter og andre parametere.
Trykk «Run» og vi får et nytt midlertidig lag «Modified geometry»
Lynnedslag med usikkerhetsellipser. De fleste har den ideelle systemusikkerheten på 200 meter, mens et par har registeringer som gir en høyere usikkerhet langs en akse.
Nå er vi en standard overlay unna å se hvilke kabler som kan ha blitt truffet av lynnedslag. Bruk Vector, Geoprocessing Tools, Intersection. Jeg velger kartet med kabler som input layer og kartet med usikkerhetsellipsene som Overlay layer.
Når dette kjøres, får vi et nytt midlertidig lag – «Intersection». Dette viser at i dette tilfellet var en kabel i faresonen, selv om akkurat dette tordenværet ikke gav noen driftsforstyrrelser.
Et kabelstrekk var i faresonen under dette tordenværet
Alt arbeidet som er gjort er lagret i midlertidige lag som forsvinner når QGIS avsluttes. Dersom vi ønsker å ta vare på noen av disse lagene, høyreklikk på det midlertidige laget, velg «Make permanent» og lagre det som et passende geodatasett. Eventuelt for postgis eller andre geodatabasebrukere: Gå til DB Manager, velg aktuell database (og schema for postgis, og muligens for andre) og «Import layer file».
Kommer snart: Hvordan få dataene rett inn i postgis og går det montro an å lage en animasjon ut av dette? – stay tuned.
Jeg har fortsatt planer om å skrive om analyse, lagring og animasjon av de arkiverte lyndataene fra Meteorologisk Institutt, men jeg fikk litt inside informasjon om at det finnes en SSE server på lyn.met.no/events som sender ut i sanntid de samme dataene som vises på lyn.met.no. Denne spytter ut fortløpende forenklede lyndata i et tekstformat som ser ut som dette:
(heartbeat blir sendt ut når det har gått for lang tid uten lyndata)
Dette var litt for kult til å la ligge…
For å få dette inn i QGIS, må vi ha en prosess som går i bakgrunnen og oppdaterer et datasett. Vi kan få til dette ved å bruke prosesseringsgrensesnittet. Det første forsøket var med å kjøre kode direkte i pythonkonsollet. Dette fungerte som et «proof of concept» for å få inn data, men var ellers fullstendig ubrukelig da pythonkonsollet kjører i samme tråd som brukerinterfacet så programmet var fullstendig uresponsivt mens scriptet kjørte…
Først av alt må vi opprette et lag som dataene kan legges inn i. Laget må lagre punkter og ha to felt. Et felt fid, integer og et felt Epoch, tekst, minst 30 tegn. (Noen dataformater forlanger at det finnes en feltid, fid). Når dette kartlaget er satt opp, kan prosesseringsscriptet kjøres. Last ned http://qgis.no/filer/lynprocess.py og åpne dette i prosesseringsverktøykassen.
Scriptet vil bli åpnet i et editor-vindu. Trykk på den grønne startpilen for å starte opp, da kommer prosesseringsdialog, velg kartlaget som ble laget og start prosesseringen. Hvis det er tordenværområder innen Skandinavia, vil lynregistreringene begynne å komme opp på kartet som punkter.
For å få et plott tilsvarende lyn.met.no, må det kontinuerlig beregnes hvor gamle registeringene er. Dette kan gjøres med et kalkulert felt. Gå inn i attributt-tabellen. Trykk ctrl-i eller kulerammeikonet for å åpne feltkalkulatoren. Hak av for å opprette et virtuelt felt, gi det et passende navn, «Age» for eksempel, og velg «Rediger funksjon»-fanen. Fjern alt som står i det øverste tekst-vinduet og skriv inn:
from qgis.core import *
from qgis.gui import *
from datetime import datetime
@qgsfunction(args='auto', group='Custom')
def ageOfEvent(value, feature, parent):
"""
reads value as a timestamp, calculates age
2021-06-20T18:39:48.624772864Z
"""
eventtime=dt_object1 = datetime.strptime(value[:25], "%Y-%m-%dT%H:%M:%S.%f")
age=(datetime.utcnow()-eventtime).total_seconds()
return(age)
(I funksjonen kutter jeg nøyaktigheten på registeringstidspunktet fra nanosekunder til mikrosekunder siden dette er oppløsningen et python timestamp har, det er greit nok for det det skal brukes til her)
Trykk «Lagre og last funksjoner» og gå tilbake til Uttrykk. Skriv så inn
ageOfEvent("Epoch")
Hvis dette har fungert som det skal, skal det stå et tall i forhåndsvisningen nederst. Dette er alderen på en av registeringene i sekunder. Dette tallet kan brukes som grunnlag for en gradert symbologi.
0-46 sekunder er en stor gul ring – grensene er tilfeldige verdier fra datasettet jeg hadde inne da jeg satt det opp
Ved å gå på «Styr objektenes opptegningsrekkefølge» kan vi sørge for at de nyeste registreringene alltid kommer øverst. Sett uttrykk til «Epoch» og Stigende/Synkende til Stigende.
Finn fram popcornet og vent på tordenværet… Dataene i videoen under blinker litt. Det skyldes sannsynligvis at OSM er et litt for tungt bakgrunnskart for denne bruken. Tekstterminalen med gul bakgrunn kjører et python script som viser rådataene som kommer inn, QGIS prosesseringsvinduet som ligger halvveis i bakgrunnen viser tidspunkt for siste registering. Dette er kommentert ut i scriptet slik det ligger for nedlasting.
Dersom jeg skulle hatt dette systemet i drift, ville jeg sannsynligvis valgt å splitte det i to: En prosess som kjører og lytter til lyn.met.no/events og lagrer eventene fortløpende for eksempel i en postgis tabell. Deretter kan det settes opp et prosjekt i QGIS som leser dette datasettet og oppdaterer med passende frekvens. Dette vil sannsynligvis gi noe som er enklere å vedlikeholde, men en en ren QGIS-løsning vil ofte være enklere og raskere å få opp. Det kan sannsynligvis også være en ide å legge til et system i innlesingsprosessen som sletter data som har blitt for gamle, eldre enn en eller to timer for eksempel. Skal man bruke historiske lyndata er jo disse tilgjengelige som beskrevet i Skaff og bruk lynnedslagsdata fra MET, disse dataene inneholder også mer metainformasjon.
Meteorologisk Institutt har en policy med mest mulig åpne data. Da må¨vi kunne hente dem inn og bruke dem i QGIS. Et interessant datasett er lyndata. Sanntidsdata vises på lyn.met.no, men arkiverte data kan i likhet med mange andre meteorologiske datasett hentes fra https://frost.met.no/api.html#!/lightning/getLightning. For å hente data derifra, må man registere seg og opprette en brukerid (beskrevet under «How to use Frost«). Når man har denne id-en, kan man hente ned data fra nettsiden. Det gir en tekstfil i såkalt ualf-format (https://opendata.smhi.se/apidocs/lightning/parameters.html). For at den skal være grei å bruke videre, lim inn denne linjen som den første:
d year month day hour min sec nanos lat lon pI multi nsens dof angle major minor chi2 rt ptz mrr cloud aI sI tI
Eventuelt kjør dette pythonscriptet (testet i python3.6.8, men bør fungere i alle versjoner >= 3):
#!/usr/bin/python3
import requests
client_id='<Skaffes fra frost.met.no>'
reftime='latest'
maxage='P1D'
# Kombinasjonen av reftime og maxage gir lynregistereringer fra siste døgn
# Bruk reftime='start/end' for et tidsintervall,
# fx '2021-05-13/2021-05-14' for lyndata fra 13. mai '21.
# I så fall blir maxage ignorert
polygon='10.3227 61.4215,12.0436 61.4712,12.0112 60.4963,10.3680 60.6211'
# Dette er et område på Hedmarken
endpoint= 'https://frost.met.no/lightning/v0.ualf?referencetime={}&maxage={}&geometry=POLYGON(({}))'.format(reftime,maxage,polygon)
header='d year month day hour min sec nanos lat lon pI multi nsens dof angle major minor chi2 rt ptz mrr cloud aI sI tI'
r=requests.get(endpoint, auth=(client_id,''))
print(header)
print(r.text)
Dette gav en output som begynte med disse linjene da det ble kjørt 10. juni ’21 kl 13:50
Dette kan lagres og leses rett inn i QGIS som en avgrenset tekstfil.
Innlesing av lyndata. Må sette opp som mellomromseparert. QGIS vil gjette at lon er X field og lat er Y field. Pass på at CRS er riktig
Med kartverkets TOPO4 gråtone og CRS satt til ETRS89 / UTM zone 32N, får jeg dette kartet:
Lynkart med alle registreringer fra MET
Dette viser alle lynregisteringer i det aktuelle tidsrommet. Det interessante er ofte lynnedslag. Da må vi filtrere datasettet, Høyreklikk på settet og velg «Filter». I uttrykksboksen, skriv «cloud»=0
Her filterer vi så bare lynnedslag vises. «cloud»=1 vil da tilsvarende vise sky-sky lyn (Jeg bruker en litt eldre fil, så feltnavnene stemmer ikke 100% med det som beskrives i teksten)
Med dette forsvinner de aller fleste registreringene:
Lynkart med bare lynnedslag fra MET
Det går an å bruke flere kolonner, noen av de mest interessante er kanskje angle, major og minor som definerer en usikkerhetsellipse for lynregisteringen. pI kolonnen er maks strøm i lynet og med alt jeg har gjort med elektronikk og elektro i mitt liv er dette første gangen jeg har brukt enheten kiloAmpere.
Polygonargumentet til api’et er et polygon i WKT-format. Skriv inn lon og lat verdier skilt med mellomrom for hjørnene i rekkefølge rundt. Komma skiller mellom hvert koordinatpar, det siste paret må gjenta det første. For eksempel, for å få alt mellom 60 og 61 grader nord og 10 og 11 grader øst: POLYGON(( 10 60, 10 61, 11 61, 11 60, 10 60)) (Doble paranteser stemmer, jeg er ikke sikker på om det er signifikant om det er med eller mot klokka, men med klokka fungerer i hvert fall)
Hvis man ikke ønsker å håndlage WKT, kan man bruke pluginen get WKT med denne er det mulig å hente ut WKT-beskrivelsen av et valgt polygon. En annen mulighet for å se at et område ble som planlagt er å bruke pluginen QuickWKT. Da er det mulig å skrive eller lime inn et WKT-polygon og få det vist på kartet. Bruker man denne, må man bare være oppmerksom på at WKT-koordinatene er forutsatt å være i samme CRS som prosjektet man arbeider med, så enten må prosjektet settes midlertidig til WGS84, eller så må projeksjonen på det WKT-laget settes til WGS84 etter at det er laget. (Layer properties, Source, Assigned Coordinate Reference System (CRS))
Dersom lyndataene brukes i noen sammenheng, husk å kreditere Meteorologisk Institutt.
Kommer snart med oppfølging, hvordan få dataene inn i postgis, plotte usikkerhetsellipsen og se om noe interessant kan ha blitt truffet. Stay tuned!
