Read lightning data from frost,met.no

readfrost.py

#!/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")