xmet/lib/xmet/sounding.py

import datetime
import shapely

from typing import Self

from xmet.db     import Database, DatabaseTable, DatabaseOrder
from xmet.coord  import COORD_SYSTEM
from xmet.list   import nearest
from xmet.series import Series, SeriesIntersection
from xmet.thermo import follow_dry_adiabat, \
                        follow_moist_adiabat, \
                        follow_saturated_mixing_ratio, \
                        pressure_height, \
                        kelvin, \
                        celsius

LAPSE_RATE_DRY   = 9.8 # degrees C per 1000m
LAPSE_RATE_MOIST = 4.0

class SoundingSample(DatabaseTable):
    __slots__ = (
        'id', 'sounding_id', 'elapsed', 'pressure', 'pressure_qa',
        'height', 'height_qa', 'temp', 'temp_qa', 'humidity',
        'dewpoint', 'wind_dir', 'wind_speed'
    )

    __table__ = 'xmet_sounding_sample'
    __key__   = 'id'

    __columns__ = (
        'id', 'sounding_id', 'elapsed', 'pressure', 'pressure_qa',
        'height', 'height_qa', 'temp', 'temp_qa', 'humidity',
        'dewpoint', 'wind_dir', 'wind_speed'
    )

    def __init__(self):
        super().__init__()
        self.id: int = None
        self.sounding_id: int   = None
        self.elapsed:     int   = None
        self.pressure:    float = None
        self.pressure_qa: str   = None
        self.height:      float = None
        self.height_qa:   str   = None
        self.temp:        float = None
        self.temp_qa:     str   = None
        self.humidity:    float = None
        self.dewpoint:    float = None
        self.wind_dir:    float = None
        self.wind_speed:  float = None

    def is_saturated(self) -> bool:
        return self.humidity >= 100.0

class Sounding(DatabaseTable):
    __slots__ = (
        'id', 'station', 'timestamp_observed', 'timestamp_released',
        'data_source_pressure', 'data_source_other', 'samples', 'location'
    )

    __table__ = 'xmet_sounding'
    __key__   = 'id'

    __columns__ = (
        'id', 'station', 'timestamp_observed', 'timestamp_released',
        'data_source_pressure', 'data_source_other', 'location'
    )

    __columns_read__ = {
        'location': 'ST_AsText(location) as location'
    }

    __values_read__ = {
        'location': shapely.from_wkt
    }

    __columns_write__ = {
        'location': 'ST_GeomFromText(:location, {crs})'.format(crs=COORD_SYSTEM)
    }

    __values_write__ = {
        'location': lambda v: {'location': shapely.to_wkt(v)}
    }

    id:                   int
    station:              str
    timestamp_observed:   datetime.datetime
    timestamp_released:   datetime.datetime
    data_source_pressure: str
    data_source_other:    str
    location:             shapely.Point
    samples:              list[SoundingSample]

    def __init__(self):
        super().__init__()
        self.id = None

    @staticmethod
    def valid_by_station(db: Database,
                         station: str,
                         timestamp: datetime.datetime=None):
        sql = """
            select
                id, station, timestamp_observed, timestamp_released,
                data_source_pressure, data_source_other,
                ST_AsText(location) as location
            from
                xmet_sounding
            where
                station = :station
                and timestamp_observed <= :timestamp
            order by
                timestamp_observed desc
            limit 1
        """

        if timestamp is None:
            timestamp = datetime.datetime.now(datetime.UTC)

        st = db.query_sql(Sounding, sql, {
            'station':   station,
            'timestamp': timestamp
        })

        sounding = st.fetchone()
        sounding.samples = list(db.query(SoundingSample,
            clauses = [
                'sounding_id = :sounding_id'
            ],
            values = {
                'sounding_id': sounding.id
            },
            order_by = [[
                'pressure', DatabaseOrder.DESC
            ]]).fetchall())

        return sounding

    @staticmethod
    def valid_by_location(db: Database,
                          location: shapely.Point,
                          timestamp: datetime.datetime):
        sql = """
            select
                id, station, timestamp_observed, timestamp_released,
                data_source_pressure, data_source_other,
                ST_AsText(location) as location,
                ST_Distance(location, MakePoint(:lon, :lat, {crs})) as distance
            from
                xmet_sounding
            where
                timestamp_observed <= :timestamp
            order by
                distance asc,
                timestamp_observed desc
            limit 1
        """.format(crs=COORD_SYSTEM)

        if timestamp is None:
            timestamp = datetime.datetime.now(datetime.UTC)

        st = db.query_sql(Sounding, sql, {
            'lon':       location.x,
            'lat':       location.y,
            'timestamp': timestamp
        })

        sounding = st.fetchone()
        sounding.samples = list(db.query(SoundingSample,
            clauses = [
                'sounding_id = :sounding_id'
            ],
            values = {
                'sounding_id': sounding.id
            },
            order_by = [[
                'pressure', DatabaseOrder.DESC
            ]]).fetchall())

        return sounding

    def follow_temp(self) -> Series:
        series = Series()

        for sample in self.samples:
            series[sample.pressure] = sample.temp

        return series

def between(n, a, b):
    return n > a and n < b

class SoundingParameters():
    @staticmethod
    def cape(temp_line: Series,
             moist_adiabat: Series,
             lfc: tuple[float],
             el: tuple[float]) -> float:
        cape = 0.0

        p_lfc, p_el = lfc[1], el[1]

        neighbors = temp_line.neighbors(moist_adiabat)
        gph_last  = None

        for pair in neighbors:
            p_env, p_parcel = pair

            if not between(p_env, p_el, p_lfc):
                continue

            t_env    = kelvin(temp_line[p_env])
            t_parcel = kelvin(moist_adiabat[p_parcel])
            gph      = pressure_height(p_env)

            if gph_last is not None:
                gph_delta = gph - gph_last

                cape += ((t_parcel - t_env) / t_env) * gph_delta

            gph_last = gph

        return 9.8076 * cape

    @staticmethod
    def from_sounding(sounding: Sounding) -> Self:
        temp_line   = sounding.follow_temp()
        dry_adiabat = follow_dry_adiabat(sounding.samples[0].temp,
                                         sounding.samples[0].pressure)

        saturated_mr_line = follow_saturated_mixing_ratio(sounding.samples[0].dewpoint,
                                                          sounding.samples[0].pressure)

        lcl = dry_adiabat.intersect(saturated_mr_line, SeriesIntersection.LESSER)

        moist_adiabat = follow_moist_adiabat(*lcl)

        lfc = moist_adiabat.intersect(temp_line, SeriesIntersection.GREATER)
        el  = moist_adiabat.intersect(temp_line, SeriesIntersection.LESSER, lfc[1])

        params = SoundingParameters()
        params.lcl  = lcl
        params.lfc  = lfc
        params.el   = el
        params.cape = SoundingParameters.cape(temp_line,
                                              moist_adiabat,
                                              lfc,
                                              el)

        return params