Skip to content

Commit

Permalink
feat: add a WIP DVN RTP product
Browse files Browse the repository at this point in the history
  • Loading branch information
@akrherz
akrherz committed Nov 17, 2023
1 parent e3a1acb commit 699c7b5
Show file tree
Hide file tree
Showing 2 changed files with 235 additions and 143 deletions.
197 changes: 123 additions & 74 deletions scripts/00z/generate_rtp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,96 +7,128 @@
import subprocess
import tempfile

import pandas as pd
from pyiem.network import Table as NetworkTable
from pyiem.tracker import loadqc
from pyiem.util import get_dbconnc, logger, utc
from pyiem.util import get_sqlalchemy_conn, logger, utc
from sqlalchemy import text

LOG = logger()
NETWORKS = ["IA_ASOS", "ISUSM", "IA_DCP", "IA_RWIS"]
PRECIP_WORKS = ["IA_ASOS", "ISUSM"]


def main():
JOBS = [
{
"wfo": "DMX",
"filename": "awos_rtp_00z.shef",
"networks": ["IA_ASOS", "ISUSM", "IA_DCP", "IA_RWIS"],
"precip_works": ["IA_ASOS", "IL_ASOS", "ISUSM"],
"limiter": "state = :state",
"filter": ["state", "IA"],
"state": "IA",
},
{
"filename": "awos_rtp_00z_dvn.shef",
"networks": "IA_ASOS ISUSM IA_DCP IA_RWIS IL_RWIS IL_ASOS".split(),
"precip_works": ["IA_ASOS", "IL_ASOS", "ISUSM"],
"limiter": "wfo = :wfo",
"filter": ["wfo", "DVN"],
"wfo": "DVN",
},
]


def main(job):
"""Go Main Go."""
qdict = loadqc()
pgconn, icursor = get_dbconnc("iem")

ets = utc().replace(hour=0, minute=0, second=0, microsecond=0)
sts12z = ets + datetime.timedelta(hours=-12)
sts6z = ets + datetime.timedelta(hours=-18)
sts24h = ets + datetime.timedelta(days=-1)
job["ets"] = utc().replace(hour=0, minute=0, second=0, microsecond=0)
job["sts12z"] = job["ets"] + datetime.timedelta(hours=-12)
job["sts6z"] = job["ets"] + datetime.timedelta(hours=-18)
job["sts24h"] = job["ets"] + datetime.timedelta(days=-1)

asosfmt = "%-6s:%-43s: %3s / %3s / %5s\n"
fmt = "%-6s:%-43s: %3s / %3s\n"

# We get 18 hour highs
highs = {}
sql = """SELECT t.id as station,
round(max(tmpf)::numeric,0) as max_tmpf,
job["badtemps"] = [sid for sid in qdict if qdict[sid].get("tmpf")]
job["badprecip"] = [sid for sid in qdict if qdict[sid].get("precip")]
with get_sqlalchemy_conn("iem") as conn:
# We get 18 hour highs
data = pd.read_sql(
text(
f"""
SELECT id, round(max(tmpf)::numeric,0) as max_tmpf,
count(tmpf) as obs FROM current_log c, stations t
WHERE t.iemid = c.iemid and t.network = ANY(%s) and valid > %s
and valid < %s
and tmpf between -70 and 130 GROUP by t.id """
icursor.execute(sql, (NETWORKS, sts6z, ets))
for row in icursor:
if qdict.get(row["station"], {}).get("tmpf"):
continue
highs[row["station"]] = row["max_tmpf"]

# 00 UTC to 00 UTC Preciptation
pcpn = {}
icursor.execute(
"""
WHERE t.iemid = c.iemid and t.network = ANY(:networks) and
valid >= :sts6z and valid < :ets and {job['limiter']}
and tmpf > -99 and not id = ANY(:badtemps) GROUP by id
"""
),
conn,
params=job,
index_col="id",
)
# 18 hour low temperature
lows = pd.read_sql(
text(
f"""
SELECT id, round(min(tmpf)::numeric,0) as min_tmpf,
count(tmpf) as obs FROM
current_log c JOIN stations t on (t.iemid = c.iemid)
WHERE t.network = ANY(:networks) and valid >= :sts6z
and valid < :ets and tmpf > -99 and not
id = ANY(:badtemps) and {job['limiter']} GROUP by id
"""
),
conn,
params=job,
index_col="id",
)
data["min_tmpf"] = lows["min_tmpf"]

# 12z to 12z precip
pcpn = pd.read_sql(
text(
f"""
select id as station, sum(precip) from
(select t.id, extract(hour from valid) as hour,
(select id, extract(hour from valid) as hour,
max(phour) as precip from current_log c, stations t
WHERE t.network = ANY(%s) and t.iemid = c.iemid
and valid >= %s and valid < %s
GROUP by t.id, hour) as foo
GROUP by id""",
(NETWORKS, sts24h, ets),
)
for row in icursor:
if qdict.get(row["station"], {}).get("precip") or row["sum"] is None:
continue
if 0 < row["sum"] < 0.009:
pcpn[row["station"]] = "T"
else:
pcpn[row["station"]] = f"{row['sum']:5.2f}"

lows = {}
icursor.execute(
"""SELECT t.id as station,
round(min(tmpf)::numeric,0) as min_tmpf,
count(tmpf) as obs FROM current_log c, stations t
WHERE t.iemid = c.iemid and t.network = ANY(%s) and valid > %s and
valid < %s and tmpf > -99 GROUP by t,id""",
(NETWORKS, sts6z, ets),
)

for row in icursor:
if qdict.get(row["station"], {}).get("tmpf"):
continue
lows[row["station"]] = row["min_tmpf"]
WHERE t.network = ANY(:networks) and t.iemid = c.iemid
and valid >= :sts24h and valid < :ets
and not id = ANY(:badprecip) and {job['limiter']}
GROUP by id, hour) as foo
GROUP by id
"""
),
conn,
params=job,
index_col="station",
)
pcpn = pcpn[pcpn["sum"].notna()]
data["precip"] = pcpn["sum"]

with tempfile.NamedTemporaryFile("w", delete=False) as fh:
fh.write("\n\n\n")
for netid in NETWORKS:
networkname = "AWOS" if netid == "IA_ASOS" else netid
precip_on = netid in PRECIP_WORKS
for netid in job["networks"]:
nt = NetworkTable(netid)
ids = list(nt.sts.keys())
ids.sort()
for sid, meta in nt.sts.items():
if meta[job["filter"][0]] != job["filter"][1]:
ids.remove(sid)
if not ids:
continue
networkname = netid.replace("ASOS", "AWOS")
precip_on = netid in job["precip_works"]
fh.write(
f".BR DMX {ets:%m%d} Z "
f".BR {job['wfo']} {job['ets']:%m%d} Z "
"DH00/TAIRVS/TAIRVI"
f"{'/PPDRVZ' if precip_on else ''}\n"
f": IOWA {networkname} RTP FIRST GUESS PROCESSED BY THE IEM\n"
f": 06Z TO 00Z HIGH TEMPERATURE FOR {sts12z:%d %b %Y}\n"
f": 06Z TO 00Z LOW TEMPERATURE FOR {sts12z:%d %b %Y}\n"
f": {networkname} RTP FIRST GUESS PROCESSED BY THE IEM\n"
": 06Z TO 00Z HIGH TEMPERATURE FOR "
f"{job['sts12z']:%d %b %Y}\n"
": 06Z TO 00Z LOW TEMPERATURE FOR "
f"{job['sts12z']:%d %b %Y}\n"
": 00Z YESTERDAY TO 00Z TODAY RAINFALL\n"
": ...BASED ON REPORTED OBS...\n"
)
nt = NetworkTable(netid)
ids = list(nt.sts.keys())
ids.sort()
for sid in ids:
if (
netid == "IA_ASOS"
Expand All @@ -105,11 +137,28 @@ def main():
continue
if netid == "IA_DCP" and nt.sts[sid]["name"].find("RAWS") < 0:
continue
myP = pcpn.get(sid, "M")
myH = highs.get(sid, "M")
myL = lows.get(sid, "M")
my_high = "M"
my_low = "M"
my_precip = "M"
if sid in data.index:
row = data.loc[sid]
if not pd.isna(row["max_tmpf"]):
my_high = int(row["max_tmpf"].round(0))
if not pd.isna(row["min_tmpf"]):
my_low = int(row["min_tmpf"].round(0))
if not pd.isna(row["precip"]):
if 0 < row["precip"] < 0.009:
my_precip = "T"
else:
my_precip = f"{row['precip']:.02f}"
_fmt = asosfmt if precip_on else fmt
args = [sid, nt.sts[sid]["name"], myH, myL, myP]
args = [
sid,
nt.sts[sid]["name"],
my_high,
my_low,
my_precip,
]
if not precip_on:
args = args[:-1]
fh.write(_fmt % tuple(args))
Expand All @@ -119,14 +168,14 @@ def main():
cmd = [
"pqinsert",
"-p",
f"plot ac {ets:%Y%m%d}0000 awos_rtp_00z.shef awos_rtp_00z.shef shef",
f"plot ac {job['ets']:%Y%m%d}0000 {job['filename']} "
f"{job['filename']} shef",
fh.name,
]
LOG.info(" ".join(cmd))
subprocess.call(cmd)
os.unlink(fh.name)
pgconn.close()


if __name__ == "__main__":
main()
[main(job) for job in JOBS]
Loading

0 comments on commit 699c7b5

Please sign in to comment.