Skip to content

Instantly share code, notes, and snippets.

@sickel
Last active Sep 10, 2021
Embed
What would you like to do?
Read lightning data from frost,met.no
#!/usr/bin/python3
"""Reads lightning data from frost.met.no. Inserts all points and auxilliary data into a PostGIS database. After allt
the points are inserted, uncertainty ellipses are calculated for each point and stored in the same record. In QGIS it is
then possible to use a calculated geometry for either the point or the ellipse"""
# Using format strings, so needs python >= 3.6
# If this file is stored in a version control system, rewrite the next two lines to
# read clientid and connection parameters from a file that is not checked in in the VCS
client_id="FROSTCLIENTID"
# See https://frost.met.no/howto.html about how to get the client id.
dbconn="dbname=<database> user=<username> password=<password> port=<port> host=<hostname>"
# any value can be undefined. Default values:
# dbname and user: logged in username
# password: blank (passwordless login)
# port: 5432
# host: localhost
CRS=25833
# The CRS in which the data are to be stored. This must be a meter based CRS.
# For a normal setup, it should not be neccesary to change anything below here
import requests
import sys
import os
import psycopg2
import psycopg2.extras
import argparse
parser = argparse.ArgumentParser(description='Fetches lightning data from frost.met.no, see https://frost.met.no/api.html#!/lightning/getLightning')
parser.add_argument('-r','--reftime',help="The time range to get observations for in either extended ISO-8601 format or the single word 'latest'.",default='latest')
parser.add_argument('-m','--maxage',help="The maximum observation age as an ISO-8601 period, like 'P1D'. Applicable only when reftime=latest. default P1D", default='P1D')
parser.add_argument('-p','--polygon',help="Coordinates of the area of interest")
parser.add_argument('-q','--quiet',help="1, no summary per line, 2 no output at all",default=0,type=int)
args = parser.parse_args()
reftime=args.reftime
maxage=args.maxage
quiet=args.quiet
if args.polygon is None:
polygon='10 61,12 61,12 60,10 60, 10 61'
else:
polygon=args.polygon
endpoint= f'https://frost.met.no/lightning/v0.ualf?referencetime={reftime}&maxage={maxage}&geometry=POLYGON(({polygon}))'
CONN = psycopg2.connect(dbconn,cursor_factory=psycopg2.extras.DictCursor)
CUR = CONN.cursor()
if quiet<1:
print('db connected ok')
SQL = f"""insert into lightning (timestamp, point, peakcurr, multiplicy, sensors, dof, angle, semimaj, semimin, chi2, risetime,
peak2zerotime, maksror, cloud, angleind, signalind, timingind)
values (%s, ST_Transform(ST_setSRID(ST_MakePoint(%s,%s),4326),{CRS}), %s,
%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s)"""
# Reprojecting the points as etrs89/UTMzone33N since I need that projection for my work later on
CREATEELLIPSES=f"""update lightning
set ellipse = st_translate(st_rotate(st_scale(st_buffer(st_setsrid(st_point(0::double precision, 0::double precision), {CRS}),
1000::double precision), lightning.semimin::double precision, lightning.semimaj::double precision),
lightning.angle::double precision * pi() / '-180'::integer::double precision), st_x(lightning.point), st_y(lightning.point))
where ellipse is null"""
# Basically: Make a point at 0,0 in the coordinate system, make a 1 km buffer around the point (since the axis are given in kms),
# scale the buffer as defined by semi min and semimaj, rotate it as definied by angle
# and in the end, translate it to where it is to sit
# =====================================
#All set up, ready to get data:
r=requests.get(endpoint, auth=(client_id,''))
# The api returns data as UALF format text, https://opendata.smhi.se/apidocs/lightning/parameters.html
data = r.text.split('\n')
nread=0
nins=0
ndata=len(data)-1 # Always an empty line at the end, may be the only line received
for strike in data:
if quiet < 1:
print(f"{nread} read / {nins} inserted / {ndata}")
if strike=='':
continue
nread+=1
params=strike.split(" ")
if len(params)<=1:
continue
time=f'{params[1]}-{params[2]}-{params[3]}T{params[4]}:{params[5]}:{int(params[6])+int(params[7])/1E9}+00'
useparams=params[10:]
useparams.insert(0,params[8]) #Need to switch position of lat and lon
useparams.insert(0,params[9])
useparams.insert(0,time)
# Now useparams contains time, lon, lat, and then the data from position 9 and onwards
try:
CUR.execute(TESTSQL,useparams[:3])
n=CUR.fetchall()[0][0]
if n > 0:
continue
CUR.execute(SQL,useparams)
nins+=1
CONN.commit() # Commit in case error on later insert, maybe a bit conservative as existing key is checked above
except psycopg2.errors.UniqueViolation:
if quiet<1:
print("Already inserted")
CONN.commit() # As an error has happened, needs to start a new transaction
continue
except Exception as e:
print(useparams)
print(e)
sys.exit(2)
CONN.commit()
if quiet < 1:
print("Creating ellipses")
CUR.execute(CREATEELLIPSES)
CONN.commit()
if quiet < 2:
print(f"{nread} read / {nins} inserted / {ndata} total")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment