Meteorologist: MQL recipes

Recipes

Dew point and cumulus base
Heat index of a heat wave
Wind chill table and frostbite threshold
Reduce a station pressure to sea level
100 year rainfall with the Gumbel law
Downscale a model temperature to a station
Vector mean wind, steadiness and wind rose of a station
Geostrophic and gradient wind from a 500 hPa chart
Static stability of a sounding layer
Verification of two temperature forecast models

executed on MentDB

Dew point and cumulus base

From the afternoon temperature in degF and the relative humidity, compute the dew point with the Magnus formula and estimate the base of the cumulus clouds.

Pipeline

01 · convertdegF to degCunit convert→02 · magnusDew pointmath eval→03 · cloud125 m per degree of spreadmath calc

Script

admin

#1. Afternoon observation: 86 degF and 62 percent relative humidity;
-> "[T]" (unit convert 86 "degF" "degC");
-> "[RH]" 62;
#2. Magnus formula for the dew point, constants 17.62 and 243.12 degC;
-> "[g]" (math eval "ln(RH/100)+a*T/(b+T)" (concat "{\"RH\": " [RH] ", \"T\": " [T] ", \"a\": 17.62, \"b\": 243.12}"));
-> "[Td]" (math calc "243.12*[g]/(17.62-[g])");
#3. Base of cumulus clouds: about 125 m per degree of spread;
-> "[base]" (math calc "125*([T]-[Td])");
concat "T = " (math round [T] 1) " degC, dew point = " (math round [Td] 1) " degC = " (math round (unit convert [Td] "degC" "degF") 1) " degF, cloud base = " (math round [base] 0) " m = " (math round (unit convert [base] "m" "ft") -2) " ft"

Result

mentdb

T = 30 degC, dew point = 21.9 degC = 71.5 degF, cloud base = 1010 m = 3300 ft

Go further

Compute the relative humidity back from T and Td
Get the wet-bulb temperature with num solve

executed on MentDB

Heat index of a heat wave

Apply the NWS Rothfusz regression to a 35 degC, 55 percent humidity afternoon and classify the risk.

Pipeline

01 · convertdegC to degFunit convert→02 · indexRothfusz regressionmath calc→03 · classifyNWS categoryif

Script

admin

#1. Heat wave: 35 degC and 55 percent relative humidity, the NWS formula works in degF;
-> "[T]" (unit convert 35 "degC" "degF");
-> "[RH]" 55;
#2. Rothfusz regression of the heat index;
-> "[HI]" (math calc "-42.379+2.04901523*[T]+10.14333127*[RH]-0.22475541*[T]*[RH]-0.00683783*[T]^2-0.05481717*[RH]^2+0.00122874*[T]^2*[RH]+0.00085282*[T]*[RH]^2-0.00000199*[T]^2*[RH]^2");
-> "[HIc]" (unit convert [HI] "degF" "degC");
#3. NWS category;
-> "[cat]" "Caution";
if (> [HI] 90) {-> "[cat]" "Extreme caution";};
if (> [HI] 103) {-> "[cat]" "Danger";};
if (> [HI] 124) {-> "[cat]" "Extreme danger";};
concat "Heat index = " (math round [HI] 1) " degF = " (math round [HIc] 1) " degC: " [cat]

Result

mentdb

Heat index = 108.9 degF = 42.7 degC: Danger

Go further

Add the low and high humidity adjustments of the NWS
Run it on a day of hourly data

executed on MentDB

Wind chill table and frostbite threshold

At -15 degC, build the wind chill table of Environment Canada for 10 to 50 km/h and find the wind speed where the index reaches -27, the frostbite threshold.

Pipeline

01 · formulaIndex as a function of Vconcat→02 · tableIndex from 10 to 50 km/hnum table→03 · solveIndex = -27num solve

Script

admin

#1. Air temperature -15 degC, wind chill index of Environment Canada, wind V in km/h;
-> "[T]" -15;
-> "[wci]" (concat "13.12+0.6215*" [T] "-11.37*V^0.16+0.3965*" [T] "*V^0.16");
#2. Table for winds of 10 to 50 km/h;
-> "[tab]" (num table (concat "round(10*(" [wci] "))/10") "V" 10 50 10);
#3. Wind speed from which frostbite is possible in 30 min, index -27;
-> "[V27]" (num solve (concat [wci] "+27") "V" 5 100);
concat "wind chill table, V and index: " [tab] ", index -27 reached at " (math round [V27] 1) " km/h"

Result

mentdb

wind chill table, V and index: [
  [10,-21.2],
  [20,-24.2],
  [30,-26],
  [40,-27.4],
  [50,-28.6]
], index -27 reached at 36.5 km/h

Go further

Add the temperature dimension with a loop
Convert the winds from knots with unit convert

executed on MentDB

Reduce a station pressure to sea level

A mountain station at 1609 m reads 836.2 hPa: reduce it to sea level with the standard lapse rate and give the value in inches of mercury.

Pipeline

01 · convertdegC to Kunit convert→02 · reduceBarometric formulamath calc→03 · reporthPa to inHgunit convert

Script

admin

#1. Mountain station at 1609 m: station pressure 836.2 hPa, temperature 12 degC;
-> "[h]" 1609;
-> "[p]" 836.2;
-> "[T]" (unit convert 12 "degC" "K");
#2. Reduction to sea level with the standard lapse rate 0.0065 K/m;
-> "[qff]" (math calc "[p]*(1-0.0065*[h]/([T]+0.0065*[h]))^-5.257");
#3. Scale height and inches of mercury for the US report;
-> "[H]" (math calc "287.05*[T]/_g0");
concat "Sea level pressure = " (math round [qff] 1) " hPa = " (math round (unit convert [qff] "hPa" "inHg") 2) " inHg, scale height = " (math round (unit convert [H] "m" "km") 2) " km"

Result

mentdb

Sea level pressure = 1010.5 hPa = 29.84 inHg, scale height = 8.35 km

Go further

Get the altitude from the pressure: solve the same formula with num solve
Compute the QNH used by pilots

executed on MentDB

100 year rainfall with the Gumbel law

From 20 years of annual maximum daily rainfall, fit a Gumbel law by the method of moments, get the 100 year rainfall and the return period of a 120 mm day.

Pipeline

01 · statsMean and standard deviationstat mean→02 · gumbelScale and locationmath calc→03 · quantile90th percentilestat percentile

Script

admin

#1. Annual maximum of daily rainfall in mm over 20 years;
-> "[x]" "[62, 48, 75, 91, 55, 68, 103, 59, 71, 84, 47, 66, 79, 58, 112, 63, 70, 88, 52, 74]";
-> "[m]" (stat mean [x]);
-> "[s]" (stat stddev [x]);
#2. Gumbel law by the method of moments;
-> "[beta]" (math calc "[s]*sqrt(6)/pi");
-> "[mu]" (math calc "[m]-0.5772156649*[beta]");
#3. 100 year return level, and return period of a 120 mm day;
-> "[x100]" (math calc "[mu]-[beta]*ln(-ln(1-1/100))");
-> "[T120]" (math calc "1/(1-exp(-exp(-(120-[mu])/[beta])))");
concat "mean = " (math round [m] 1) " mm, s = " (math round [s] 2) " mm, 90th percentile = " (math round (stat percentile [x] 90) 1) " mm, 100 yr rain = " (math round [x100] 1) " mm, 120 mm returns every " (math round [T120] 0) " years"

Result

mentdb

mean = 71.3 mm, s = 17.57 mm, 90th percentile = 92.2 mm, 100 yr rain = 126.4 mm, 120 mm returns every 63 years

Go further

Test the fit with stat ks_test
Fit a GEV law with stat fit_model on the empirical frequencies

executed on MentDB

Downscale a model temperature to a station

A forecast model gives the 2 m temperature on a 0.25 degree grid with its own terrain height. The script interpolates the temperature and the model terrain at a valley station (bilinear and bicubic), then corrects the altitude difference with the standard lapse rate.

Pipeline

01 · linearBilinear temperaturenum interpolate2→02 · cubicBicubic temperaturenum interpolate2→03 · terrainModel height at the stationnum interpolate2→04 · lapseLapse rate correctionmath calc

Script

admin

#1. Forecast model grid 0.25 deg: latitudes, longitudes, 2 m temperature in degC and model terrain height in m;
-> "[lat]" "[45.0, 45.25, 45.5, 45.75]"; -> "[lon]" "[5.0, 5.25, 5.5, 5.75]";
-> "[T]" "[[8.2, 7.1, 4.3, 2.0], [7.6, 6.0, 3.1, 1.2], [7.9, 5.4, 2.2, 0.4], [8.4, 6.2, 3.0, 1.1]]";
-> "[Z]" "[[310, 480, 920, 1350], [350, 610, 1120, 1480], [330, 700, 1290, 1610], [290, 560, 1100, 1450]]";
#2. Station at 45.32 N, 5.41 E, altitude 214 m: bilinear and bicubic temperature;
-> "[Tl]" (num interpolate2 [lat] [lon] [T] 45.32 5.41 "linear");
-> "[Tc]" (num interpolate2 [lat] [lon] [T] 45.32 5.41 "cubic");
#3. Model terrain height at the station;
-> "[zg]" (num interpolate2 [lat] [lon] [Z] 45.32 5.41 "cubic");
#4. Altitude correction with the standard lapse rate 6.5 K/km;
concat "bilinear " (math round [Tl] 2) " degC, bicubic " (math round [Tc] 2) " degC, model height " (math round [zg] 0) " m, corrected temperature at 214 m = " (math round (math calc "[Tc]+0.0065*([zg]-214)") 1) " degC"

Result

mentdb

bilinear 3.92 degC, bicubic 3.82 degC, model height 991 m, corrected temperature at 214 m = 8.9 degC

Go further

Use the mode nearest to see the error of a nearest grid point forecast.
Several stations at once: pass JSON arrays of latitudes and longitudes.

executed on MentDB

Vector mean wind, steadiness and wind rose of a station

From twelve hourly wind reports, the script splits each wind into u and v components, gets the vector mean wind and its direction with atan2, the steadiness of the wind, the prevailing 45 deg sector with the mode and the frequency table of the sectors (value 0 = N, 1 = NE ... 7 = NW).

Pipeline

01 · splitReports and sectorsatom get→02 · meanVector mean windmath hypot→03 · directionDirection it comes frommath atan2→04 · prevailingPrevailing sectorstat mode→05 · roseFrequency of the sectorsstat frequency

Script

admin

#1. Hourly wind reports of a coastal station, direction in deg and speed in kt;
-> "[obs]" "250:12,260:15,240:10,270:18,250:14,300:8,250:11,320:6,260:16,250:13,90:4,240:9";
-> "[names]" "N,NE,E,SE,S,SW,W,NW";
-> "[su]" 0; -> "[sv]" 0; -> "[ss]" 0; -> "[sect]" "";
#2. Wind components u and v, the wind blows from the direction, and the 45 deg sector of each report;
for (-> "[i]" 1) (<= [i] (atom size [obs] ",")) (++ "[i]") {
	-> "[o]" (atom get [obs] [i] ",");
	-> "[d]" (math deg_to_rad (atom get [o] 1 ":")); -> "[s]" (atom get [o] 2 ":");
	-> "[su]" (math calc "[su]-[s]*sin([d])"); -> "[sv]" (math calc "[sv]-[s]*cos([d])"); -> "[ss]" (math calc "[ss]+[s]");
	-> "[sect]" (concat [sect] (if (> [i] 1) {","} {""}) (math calc (concat "mod(" (math round (math calc (concat (atom get [o] 1 ":") "/45")) 0) ", 8)")));
};
#3. Vector mean wind, direction it comes from, steadiness = vector mean / scalar mean;
-> "[n]" (atom size [obs] ",");
-> "[vm]" (math calc (concat (math hypot [su] [sv]) "/[n]"));
-> "[dir]" (math calc (concat "mod(" (math rad_to_deg (math atan2 (math calc "-[su]") (math calc "-[sv]"))) "+360, 360)"));
#4. Prevailing sector with the mode, frequency table of the sectors;
-> "[mode]" (string regexp_extract (stat mode (concat "[" [sect] "]")) "[0-9]+" 0);
concat "vector mean wind " (math round [dir] 0) " deg / " (math round [vm] 1) " kt, steadiness " (math round (math calc (concat [vm] "*[n]/[ss]")) 3) ", prevailing sector " (atom get [names] (math calc "[mode]+1") ",") [_n_] (stat frequency (concat "[" [sect] "]"))

Result

mentdb

vector mean wind 259 deg / 10.1 kt, steadiness 0.89, prevailing sector W
[
  {"value": 2, "count": 1, "frequency": 0.08333333333333333, "cumulative_count": 1, "cumulative_frequency": 0.08333333333333333},
  {"value": 5, "count": 2, "frequency": 0.16666666666666666, "cumulative_count": 3, "cumulative_frequency": 0.25},
  {"value": 6, "count": 7, "frequency": 0.5833333333333334, "cumulative_count": 10, "cumulative_frequency": 0.8333333333333334},
  {"value": 7, "count": 2, "frequency": 0.16666666666666666, "cumulative_count": 12, "cumulative_frequency": 1}
]

Go further

Draw the wind rose with plot bar on the counts of the frequency table.
Weight the sectors by the speed to get the energy rose of a wind farm site.

executed on MentDB

Geostrophic and gradient wind from a 500 hPa chart

The forecaster reads the height gradients on the 500 hPa chart. The script gets the Coriolis parameter, the geostrophic wind components, speed and direction, the Rossby number of the trough and the slower gradient wind around a low.

Pipeline

01 · coriolisf = 2 Omega sin phimath calc→02 · geostrophicSpeed of the geostrophic windmath hypot→03 · directionMeteorological directionmath atan2→04 · gradientGradient wind, cyclonic caseunit convert

Script

admin

#1. 500 hPa chart at 45 N: height falls 40 m over 300 km to the north and falls 60 m over 300 km to the east;
-> "[phi]" (math deg_to_rad 45); -> "[g]" 9.80665;
-> "[dzdy]" (math calc "-40/300e3"); -> "[dzdx]" (math calc "-60/300e3");
#2. Coriolis parameter f = 2 Omega sin phi;
-> "[f]" (math calc "2*7.2921e-5*sin([phi])");
#3. Geostrophic wind: ug = -g/f dZ/dy, vg = g/f dZ/dx, speed and meteorological direction;
-> "[ug]" (math calc "-[g]/[f]*[dzdy]"); -> "[vg]" (math calc "[g]/[f]*[dzdx]");
-> "[Vg]" (math hypot [ug] [vg]);
-> "[dir]" (math calc (concat "mod(" (math rad_to_deg (math atan2 (math calc "-[ug]") (math calc "-[vg]"))) "+360, 360)"));
#4. Rossby number for a 1000 km trough and gradient wind around a low of radius 800 km;
-> "[Ro]" (math calc "[Vg]/([f]*1000e3)");
-> "[R]" 800e3;
-> "[Vgr]" (math calc "-[f]*[R]/2+sqrt(([f]*[R])^2/4+[f]*[R]*[Vg])");
concat "f = " (math round_sig [f] 4) " s-1, ug = " (math round [ug] 2) " m/s, vg = " (math round [vg] 2) " m/s, geostrophic wind " (math round [dir] 0) " deg / " (math round (unit convert [Vg] "m/s" "kn") 1) " kt, Rossby = " (math round [Ro] 3) ", gradient wind around the low = " (math round (unit convert [Vgr] "m/s" "kn") 1) " kt"

Result

mentdb

f = 0.0001031 s-1, ug = 12.68 m/s, vg = -19.02 m/s, geostrophic wind 326 deg / 44.4 kt, Rossby = 0.222, gradient wind around the low = 36.2 kt

Go further

Do the anticyclonic case: the root exists only if the pressure gradient is weak enough.
Compute the thermal wind between 1000 and 500 hPa from the thickness gradient.

executed on MentDB

Static stability of a sounding layer

Between the 850 and 700 hPa levels of a radiosonde, the script gets the potential temperatures, the lapse rate, the Brunt-Vaisala frequency and the buoyancy period, and classifies the layer with the signs of N^2 and of the lapse rate minus the moist adiabatic rate.

Pipeline

01 · kelvinCelsius to kelvinunit convert→02 · thetatheta = T (1000/p)^0.2857math power→03 · stabilitySign of N^2math sign→04 · moistAgainst the moist ratemath signum

Script

admin

#1. Radiosonde: 850 hPa at 1457 m, 14.2 degC, 700 hPa at 3012 m, 2.6 degC;
-> "[T1]" (unit convert 14.2 "degC" "K"); -> "[T2]" (unit convert 2.6 "degC" "K");
-> "[z1]" 1457; -> "[z2]" 3012;
#2. Potential temperature theta = T 1000/p ^ R/cp;
-> "[th1]" (math calc (concat [T1] "*" (math power (math calc "1000/850") 0.2857)));
-> "[th2]" (math calc (concat [T2] "*" (math power (math calc "1000/700") 0.2857)));
#3. Lapse rate and Brunt-Vaisala frequency N^2 = g/theta dtheta/dz;
-> "[gamma]" (math calc "-([T2]-[T1])/([z2]-[z1])*1000");
-> "[N2]" (math calc "9.80665/(([th1]+[th2])/2)*([th2]-[th1])/([z2]-[z1])");
#4. Stability: sign of N^2 for dry air, then the lapse rate against the moist rate 6 K/km;
-> "[dry]" (math sign [N2]);
-> "[moist]" (math signum (math calc "[gamma]-6"));
concat "theta850 = " (math round [th1] 2) " K, theta700 = " (math round [th2] 2) " K, lapse rate = " (math round [gamma] 2) " K/km, N = " (math round (math calc "sqrt(abs([N2]))") 5) " s-1, buoyancy period = " (math round (math calc "2*pi/sqrt(abs([N2]))/60") 1) " min, layer: " (if (> [dry] 0) {"stable for dry air"} {"dry unstable"}) (if (> [moist] 0) {", conditionally unstable"} {", absolutely stable"})

Result

mentdb

theta850 = 301.01 K, theta700 = 305.33 K, lapse rate = 7.46 K/km, N = 0.00948 s-1, buoyancy period = 11 min, layer: stable for dry air, conditionally unstable

Go further

Load the full sounding as JSON and loop over the layers to find the inversions.
Add the humidity to get the equivalent potential temperature and the convective instability.

executed on MentDB

Verification of two temperature forecast models

Ten days of maximum temperature forecasts from two models are checked against the observations: bias, MSE split into bias squared plus error variance, RMSE, mean absolute error, and a Mann-Whitney test to see whether the absolute errors of model A are really smaller.

Pipeline

01 · errorsForecast - observationvector sub→02 · decomposeMSE = bias^2 + variancestat pvariance→03 · maeMean absolute errormath avg→04 · testNon-parametric teststat mann_whitney

Script

admin

#1. Daily maximum temperature: observations and two forecast models, degC, 10 days;
-> "[obs]" "[21.4, 23.0, 24.8, 22.1, 19.5, 18.7, 20.9, 25.3, 26.1, 24.0]";
-> "[A]" "[22.0, 23.9, 24.1, 23.4, 20.8, 18.1, 21.7, 24.2, 27.4, 24.6]";
-> "[B]" "[23.5, 23.3, 26.9, 22.5, 22.4, 16.2, 23.1, 27.8, 23.6, 26.7]";
#2. Errors forecast - observation, bias and error variance;
-> "[eA]" (vector sub [A] [obs]); -> "[eB]" (vector sub [B] [obs]);
#3. MSE = bias^2 + population variance of the error, RMSE;
-> "[mseA]" (math calc (concat (stat mean [eA]) "^2+" (stat pvariance [eA])));
-> "[mseB]" (math calc (concat (stat mean [eB]) "^2+" (stat pvariance [eB])));
#4. Absolute errors of each model, element by element from JSON;
json load "rx_p4e_err" (concat "{\"A\": " [eA] ", \"B\": " [eB] "}");
-> "[aA]" ""; -> "[aB]" "";
for (-> "[i]" 0) (< [i] 10) (++ "[i]") {
	-> "[sep]" (if (> [i] 0) {","} {""});
	-> "[aA]" (concat [aA] [sep] (math abs (json select "rx_p4e_err" (concat "/A/[" [i] "]"))));
	-> "[aB]" (concat [aB] [sep] (math abs (json select "rx_p4e_err" (concat "/B/[" [i] "]"))));
};
json unload "rx_p4e_err";
#5. Mann-Whitney test on the absolute errors: is model A really better?;
-> "[mw]" (stat mann_whitney (concat "[" [aA] "]") (concat "[" [aB] "]"));
concat "model A: bias " (math round (stat mean [eA]) 2) ", RMSE " (math round (math sqrt [mseA]) 3) ", MAE " (math round (math avg (concat "[" [aA] "]")) 3) " / model B: bias " (math round (stat mean [eB]) 2) ", RMSE " (math round (math sqrt [mseB]) 3) ", MAE " (math round (math avg (concat "[" [aB] "]")) 3) [_n_] "Mann-Whitney on the absolute errors: " [mw]

Result

mentdb

model A: bias 0.44, RMSE 0.964, MAE 0.92 / model B: bias 1.02, RMSE 2.199, MAE 2.02
Mann-Whitney on the absolute errors: {"u": 20, "p_value": 0.025414885134874533, "method": "asymptotic"}

Go further

Add the skill score against persistence: the forecast of today is the observation of yesterday.
Use stat wilcoxon instead: the two models forecast the same days, the samples are paired.